diff --git a/analysis/data/badchannels2021.dat b/analysis/data/badchannels2021.dat new file mode 100644 index 000000000..aaf20c467 --- /dev/null +++ b/analysis/data/badchannels2021.dat @@ -0,0 +1,873 @@ +svt_channel_id +470 +802 +1145 +1320 +2371 +2448 +3443 +3444 +3966 +4737 +4739 +4745 +4747 +4749 +4751 +4753 +4755 +5364 +5374 +6624 +7819 +7966 +8001 +8574 +9203 +11776 +11777 +11778 +11779 +11780 +11781 +11782 +11783 +11784 +11785 +11786 +11787 +11788 +11789 +11790 +11791 +11792 +11793 +11794 +11795 +11796 +11797 +11798 +11799 +11800 +11801 +11802 +11803 +11804 +11805 +11806 +11807 +11808 +11809 +11810 +11811 +11812 +11813 +11814 +11815 +11816 +11817 +11818 +11819 +11820 +11821 +11822 +11823 +11824 +11825 +11826 +11827 +11828 +11829 +11830 +11831 +11832 +11833 +11834 +11835 +11836 +11837 +11838 +11839 +11840 +11841 +11842 +11843 +11844 +11845 +11846 +11847 +11848 +11849 +11850 +11851 +11852 +11853 +11854 +11855 +11856 +11857 +11858 +11859 +11860 +11861 +11862 +11863 +11864 +11865 +11866 +11867 +11868 +11869 +11870 +11871 +11872 +11873 +11874 +11875 +11876 +11877 +11878 +11879 +11880 +11881 +11882 +11883 +11884 +11885 +11886 +11887 +11888 +11889 +11890 +11891 +11892 +11893 +11894 +11895 +11896 +11897 +11898 +11899 +11900 +11901 +11902 +11903 +11904 +11905 +11906 +11907 +11908 +11909 +11910 +11911 +11912 +11913 +11914 +11915 +11916 +11917 +11918 +11919 +11920 +11921 +11922 +11923 +11924 +11925 +11926 +11927 +11928 +11929 +11930 +11931 +11932 +11933 +11934 +11935 +11936 +11937 +11938 +11939 +11940 +11941 +11942 +11943 +11944 +11945 +11946 +11947 +11948 +11949 +11950 +11951 +11952 +11953 +11954 +11955 +11956 +11957 +11958 +11959 +11960 +11961 +11962 +11963 +11964 +11965 +11966 +11967 +11968 +11969 +11970 +11971 +11972 +11973 +11974 +11975 +11976 +11977 +11978 +11979 +11980 +11981 +11982 +11983 +11984 +11985 +11986 +11987 +11988 +11989 +11990 +11991 +11992 +11993 +11994 +11995 +11996 +11997 +11998 +11999 +12000 +12001 +12002 +12003 +12004 +12005 +12006 +12007 +12008 +12009 +12010 +12011 +12012 +12013 +12014 +12015 +12016 +12017 +12018 +12019 +12020 +12021 +12022 +12023 +12024 +12025 +12026 +12027 +12028 +12029 +12030 +12031 +12032 +12033 +12034 +12035 +12036 +12037 +12038 +12039 +12040 +12041 +12042 +12043 +12044 +12045 +12046 +12047 +12048 +12049 +12050 +12051 +12052 +12053 +12054 +12055 +12056 +12057 +12058 +12059 +12060 +12061 +12062 +12063 +12064 +12065 +12066 +12067 +12068 +12069 +12070 +12071 +12072 +12073 +12074 +12075 +12076 +12077 +12078 +12079 +12080 +12081 +12082 +12083 +12084 +12085 +12086 +12087 +12088 +12089 +12090 +12091 +12092 +12093 +12094 +12095 +12096 +12097 +12098 +12099 +12100 +12101 +12102 +12103 +12104 +12105 +12106 +12107 +12108 +12109 +12110 +12111 +12112 +12113 +12114 +12115 +12116 +12117 +12118 +12119 +12120 +12121 +12122 +12123 +12124 +12125 +12126 +12127 +12128 +12129 +12130 +12131 +12132 +12133 +12134 +12135 +12136 +12137 +12138 +12139 +12140 +12141 +12142 +12143 +12144 +12145 +12146 +12147 +12148 +12149 +12150 +12151 +12152 +12153 +12154 +12155 +12156 +12157 +12158 +12159 +12160 +12161 +12162 +12163 +12164 +12165 +12166 +12167 +12168 +12169 +12170 +12171 +12172 +12173 +12174 +12175 +12176 +12177 +12178 +12179 +12180 +12181 +12182 +12183 +12184 +12185 +12186 +12187 +12188 +12189 +12190 +12191 +12192 +12193 +12194 +12195 +12196 +12197 +12198 +12199 +12200 +12201 +12202 +12203 +12204 +12205 +12206 +12207 +12208 +12209 +12210 +12211 +12212 +12213 +12214 +12215 +12216 +12217 +12218 +12219 +12220 +12221 +12222 +12223 +12224 +12225 +12226 +12227 +12228 +12229 +12230 +12231 +12232 +12233 +12234 +12235 +12236 +12237 +12238 +12239 +12240 +12241 +12242 +12243 +12244 +12245 +12246 +12247 +12248 +12249 +12250 +12251 +12252 +12253 +12254 +12255 +12256 +12257 +12258 +12259 +12260 +12261 +12262 +12263 +12264 +12265 +12266 +12267 +12268 +12269 +12270 +12271 +12272 +12273 +12274 +12275 +12276 +12277 +12278 +12279 +12280 +12281 +12282 +12283 +12284 +12285 +12286 +12287 +12288 +12289 +12290 +12291 +12292 +12293 +12294 +12295 +12296 +12297 +12298 +12299 +12300 +12301 +12302 +12303 +12304 +12305 +12306 +12307 +12308 +12309 +12310 +12311 +12312 +12313 +12314 +12315 +12316 +12317 +12318 +12319 +12320 +12321 +12322 +12323 +12324 +12325 +12326 +12327 +12328 +12329 +12330 +12331 +12332 +12333 +12334 +12335 +12336 +12337 +12338 +12339 +12340 +12341 +12342 +12343 +12344 +12345 +12346 +12347 +12348 +12349 +12350 +12351 +12352 +12353 +12354 +12355 +12356 +12357 +12358 +12359 +12360 +12361 +12362 +12363 +12364 +12365 +12366 +12367 +12368 +12369 +12370 +12371 +12372 +12373 +12374 +12375 +12376 +12377 +12378 +12379 +12380 +12381 +12382 +12383 +12384 +12385 +12386 +12387 +12388 +12389 +12390 +12391 +12392 +12393 +12394 +12395 +12396 +12397 +12398 +12399 +12400 +12401 +12402 +12403 +12404 +12405 +12406 +12407 +12408 +12409 +12410 +12411 +12412 +12413 +12414 +12415 +12461 +12526 +12630 +13311 +13889 +14143 +15290 +16063 +16127 +16767 +18868 +18869 +19968 +19969 +19970 +19971 +19972 +19973 +19974 +19975 +19976 +19977 +19978 +19979 +19980 +19981 +19982 +19983 +19984 +19985 +19986 +19987 +19988 +19989 +19990 +19991 +19992 +19993 +19994 +19995 +19996 +19997 +19998 +19999 +20000 +20001 +20002 +20003 +20004 +20005 +20006 +20007 +20008 +20009 +20010 +20011 +20012 +20013 +20014 +20015 +20016 +20017 +20018 +20019 +20020 +20021 +20022 +20023 +20024 +20025 +20026 +20027 +20028 +20029 +20030 +20031 +20032 +20033 +20034 +20035 +20036 +20037 +20038 +20039 +20040 +20041 +20042 +20043 +20044 +20045 +20046 +20047 +20048 +20049 +20050 +20051 +20052 +20053 +20054 +20055 +20056 +20057 +20058 +20059 +20060 +20061 +20062 +20063 +20064 +20065 +20066 +20067 +20068 +20069 +20070 +20071 +20072 +20073 +20074 +20075 +20076 +20077 +20078 +20079 +20080 +20081 +20082 +20083 +20084 +20085 +20086 +20087 +20088 +20089 +20090 +20091 +20092 +20093 +20094 +20095 +20863 +20991 +21117 +21119 +21121 +21183 +21247 +21311 +21313 +21315 +21439 +21501 +21503 +21567 +21631 +21697 +21759 +21761 +21821 +21823 +21887 +21889 +21949 +21951 +21953 +21977 +21978 +21979 +21980 +21981 +21982 +21983 +21984 +21985 +21986 +21987 +21988 +21989 +21990 +21991 +21992 +21993 +21994 +21995 +21996 +21997 +21998 +21999 +22000 +22001 +22002 +22003 +22004 +22005 +22006 +22007 +22008 +22009 +22010 +22011 +22012 +22013 +22014 +22015 +22271 +23872 +23874 diff --git a/analysis/include/RawSvtHitHistos.h b/analysis/include/RawSvtHitHistos.h index c6a83082a..632d8d33f 100644 --- a/analysis/include/RawSvtHitHistos.h +++ b/analysis/include/RawSvtHitHistos.h @@ -26,11 +26,12 @@ class RawSvtHitHistos : public HistoManager{ ~RawSvtHitHistos(); void DefineHistos(); - void FillHistograms(RawSvtHit* rawSvtHit,float weight = 1.,int Ireg=0,unsigned int nhit = 0,Float_t TimeDiff = -42069.0,Float_t AmpDiff = -42069.0); - void saveHistosSVT(TFile* outF,std::string folder); + void FillHistograms(RawSvtHit* rawSvtHit, float weight = 1., int Ireg = 0, unsigned int nhit = 0, Float_t TimeDiff = -42069.0, Float_t AmpDiff = -42069.0, int strip = -10000, int hitc = 0, int hitl = 0, float otherTime = 69420.0); + void saveHistosSVT(TFile* outF, std::string folder); + private: - int Event_number=0; + int Event_number = 0; int debug_ = 1; int adcs_[6]; diff --git a/analysis/src/RawSvtHitHistos.cxx b/analysis/src/RawSvtHitHistos.cxx index ef3d6f0bb..85a2d7574 100644 --- a/analysis/src/RawSvtHitHistos.cxx +++ b/analysis/src/RawSvtHitHistos.cxx @@ -12,153 +12,145 @@ RawSvtHitHistos::~RawSvtHitHistos() { } void RawSvtHitHistos::DefineHistos(){ - //Define vector of hybrid names using ModuleMapper - //Use this list to define multiple copies of histograms, one for each hybrid, from json file + // Define vector of hybrid names using ModuleMapper + // Use this list to define multiple copies of histograms, one for each hybrid, from json file std::vector hybridNames; mmapper_->getStrings(hybridNames); std::string makeMultiplesTag = "SvtHybrids"; - //std::cout<<"hello1"< hybridStrings={}; +void RawSvtHitHistos::FillHistograms(RawSvtHit* rawSvtHit, float weight, int i, unsigned int i2, Float_t TimeDiff, Float_t AmpDiff, int strip, int hitc, int hitl, float otherTime) { + std::vector hybridStrings = {}; std::string histokey; //std::cout<=10000){return;} //if(Event_number==11) std::cout<getModule()); auto lay = std::to_string(rawSvtHit->getLayer()); - swTag= mmapper_->getStringFromSw("ly"+lay+"_m"+mod); - std::string helper = mmapper_->getHwFromSw("ly"+lay+"_m"+mod); + swTag = mmapper_->getStringFromSw("ly" + lay + "_m" + mod); + std::string helper = mmapper_->getHwFromSw("ly" + lay + "_m" + mod); char char_array[helper.length()+1]; - std::strcpy(char_array,helper.c_str()); + std::strcpy(char_array, helper.c_str()); int feb = (int)char_array[1]-48; int hyb = (int)char_array[3]-48; histokey = swTag + "_SvtHybrids_getFitN_h"; - //std::cout<<"hello3"<getFitN(),weight); + + Fill1DHisto(histokey, rawSvtHit->getFitN(), weight); histokey = swTag +"_SvtHybrids_T0_h"; - //std::cout<getT0(i)<getAmp(i)<getT0err(i)<getAmpErr(i)<getT0(i),weight); - //}//else{ - // Fill1DHisto(histokey, rawSvtHit->getT0(i)-27.0,weight); - //} - //std::cout<<"hello6"<getT0(i), weight); + histokey = swTag + "_SvtHybrids_Am_h"; - Fill1DHisto(histokey, rawSvtHit->getAmp(i),weight); + Fill1DHisto(histokey, rawSvtHit->getAmp(i), weight); histokey = swTag + "_SvtHybrids_Chi_Sqr_h"; - Fill1DHisto(histokey, rawSvtHit->getChiSq(i),weight); - //std::cout<getStrip()<getChiSq(i), weight); + + histokey = swTag + "_SvtHybrids_HitsPerCluster_h"; + Fill1DHisto(histokey, hitc, weight); + + histokey = swTag + "_SvtHybrids_HitsOnALayer_h"; + Fill1DHisto(histokey, hitl, weight); + histokey = swTag + "_SvtHybrids_ADCcount_hh"; - int * adcs=rawSvtHit->getADCs(); + int * adcs = rawSvtHit->getADCs(); int maxx = 0; - for(unsigned int K=0; K<6; K++){ - if(maxxgetT0(i)),((Float_t)(adcs[K]))/(rawSvtHit->getAmp(i)),weight); } - for(unsigned int K=1; K<6; K++){ - if(feb<=1){ - Fill2DHisto(histokey,24.0*K-(rawSvtHit->getT0(i)),((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K]))/(rawSvtHit->getAmp(i)),weight); - }else{ - Fill2DHisto(histokey,24.0*K-(rawSvtHit->getT0(i)),((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K]))/(rawSvtHit->getAmp(i)),weight); + for (unsigned int K=1; K<6; K++){ + if (feb <= 1){ + Fill2DHisto(histokey, 24.0*K-(rawSvtHit->getT0(i)), ((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K]))/(rawSvtHit->getAmp(i)), weight); + } else { + Fill2DHisto(histokey, 24.0*K-(rawSvtHit->getT0(i)), ((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K]))/(rawSvtHit->getAmp(i)), weight); } //((Float_t)maxx),weight); } - //NOW THE CODE FOR THE TGRAPH SHITE - // - - - - + // The following code is for the TGraph histokey = swTag + "_SvtHybrids_ADCcountdeshift_hh"; - for(unsigned int K=1; K<6; K++){ - if(feb<=1){ - if(std::abs(rawSvtHit->getT0(i)+60)<25){ - Fill2DHisto(histokey,K,((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K])));//(rawSvtHit->getAmp(i)),weight); - }else{ - Fill2DHisto(histokey,K,((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K])));//(rawSvtHit->getAmp(i)),weight); + for (unsigned int K=1; K<6; K++){ + if (feb <= 1){ + if (std::abs(rawSvtHit->getT0(i)+60) < 25){ + Fill2DHisto(histokey, K, ((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K]))); //(rawSvtHit->getAmp(i)),weight); + } else { + Fill2DHisto(histokey, K, ((Float_t)(adcs[K])-Float_t(baseErr1_[feb][hyb][(int)(rawSvtHit->getStrip())][K]))); //(rawSvtHit->getAmp(i)),weight); } - }else{ - if(std::abs(rawSvtHit->getT0(i)+60)<25){ - Fill2DHisto(histokey,K,((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K])));//(rawSvtHit->getAmp(i)),weight); - }else{ - Fill2DHisto(histokey,K,((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K])));//(rawSvtHit->getAmp(i)),weight); + } else { + if (std::abs(rawSvtHit->getT0(i)+60) < 25){ + Fill2DHisto(histokey, K, ((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K]))); //(rawSvtHit->getAmp(i)),weight); + } else { + Fill2DHisto(histokey, K, ((Float_t)(adcs[K])-Float_t(baseErr2_[feb-2][hyb][(int)(rawSvtHit->getStrip())][K]))); //(rawSvtHit->getAmp(i)),weight); } } //((Float_t)maxx),weight); } //adcs_=rawSvtHit->getADCs(i); //Fill1DHisto(histokey, -(rawSvthit->getT0(i)),weight); + + histokey = swTag + "_SvtHybrids_T0Strip_hh"; + Fill2DHisto(histokey, rawSvtHit->getT0(i), rawSvtHit->getStrip()-strip, weight); + + histokey = swTag + "_SvtHybrids_2HitCluDiff_hh"; + Fill2DHisto(histokey, rawSvtHit->getT0(i), otherTime, weight); + + histokey = swTag + "_SvtHybrids_AmChi_Sqr_hh"; + Fill2DHisto(histokey, rawSvtHit->getChiSq(i), rawSvtHit->getAmp(i), weight); - //std::cout<<"hello7"<getChiSq(i), rawSvtHit->getT0(i), weight); + histokey = swTag + "_SvtHybrids_T0Err_hh"; - Fill2DHisto(histokey, rawSvtHit->getT0(i), rawSvtHit->getT0err(i),weight); + Fill2DHisto(histokey, rawSvtHit->getT0(i), rawSvtHit->getT0err(i), weight); + histokey = swTag + "_SvtHybrids_AmErr_hh"; - //std::cout<<"hello8"<getAmp(i), rawSvtHit->getAmpErr(i),weight); + Fill2DHisto(histokey, rawSvtHit->getAmp(i), rawSvtHit->getAmpErr(i), weight); + histokey = swTag + "_SvtHybrids_AmT0_hh"; - //std::cout<<"hello9"<getT0(i), rawSvtHit->getAmp(i),weight); + Fill2DHisto(histokey, rawSvtHit->getT0(i), rawSvtHit->getAmp(i), weight); + histokey = swTag + "_SvtHybrids_AmErrT0Err_hh"; - //std::cout<<"hello10"<getT0err(i), rawSvtHit->getAmpErr(i),weight); + Fill2DHisto(histokey, rawSvtHit->getT0err(i), rawSvtHit->getAmpErr(i), weight); histokey = swTag + "_SvtHybrids_AmT0Err_hh"; - //std::cout<<"hello10"<getT0err(i), rawSvtHit->getAmp(i),weight); + Fill2DHisto(histokey, rawSvtHit->getT0err(i), rawSvtHit->getAmp(i), weight); histokey = swTag + "_SvtHybrids_AmErrT0_hh"; - //std::cout<<"hello10"<getT0(i), rawSvtHit->getAmpErr(i),weight); + Fill2DHisto(histokey, rawSvtHit->getT0(i), rawSvtHit->getAmpErr(i), weight); - if(i==1){ + if (i==1) { histokey = swTag + "_SvtHybrids_PT1PT2_hh"; - Fill2DHisto(histokey, rawSvtHit->getT0(1),rawSvtHit->getT0(0)); - }else{ + Fill2DHisto(histokey, rawSvtHit->getT0(1), rawSvtHit->getT0(0)); + } else { histokey = swTag + "_SvtHybrids_PT1PT2_hh"; - Fill2DHisto(histokey, rawSvtHit->getT0(0),rawSvtHit->getT0(1)); + Fill2DHisto(histokey, rawSvtHit->getT0(0), rawSvtHit->getT0(1)); } - if(TimeDiff==-42069){return;} - else{ + if (TimeDiff==-42069) { return; } + else { histokey = swTag + "_SvtHybrids_TD_h"; - Fill1DHisto(histokey, TimeDiff,weight); + Fill1DHisto(histokey, TimeDiff, weight); histokey = swTag + "_SvtHybrids_T0TD_hh"; - Fill2DHisto(histokey, rawSvtHit->getT0(i),TimeDiff,weight); + Fill2DHisto(histokey, rawSvtHit->getT0(i), TimeDiff, weight); histokey = swTag + "_SvtHybrids_AmErrTD_hh"; - Fill2DHisto(histokey, rawSvtHit->getAmpErr(i),TimeDiff,weight); + Fill2DHisto(histokey, rawSvtHit->getAmpErr(i), TimeDiff, weight); histokey = swTag + "_SvtHybrids_AmpTD_hh"; - Fill2DHisto(histokey, rawSvtHit->getAmp(i),TimeDiff,weight); + Fill2DHisto(histokey, rawSvtHit->getAmp(i), TimeDiff, weight); histokey = swTag + "_SvtHybrids_Amp12_hh"; - Fill2DHisto(histokey, rawSvtHit->getAmp(0),rawSvtHit->getAmp(1),weight); + Fill2DHisto(histokey, rawSvtHit->getAmp(0), rawSvtHit->getAmp(1), weight); histokey = swTag + "_SvtHybrids_ADTD_hh"; - Fill2DHisto(histokey, AmpDiff,TimeDiff,weight); + Fill2DHisto(histokey, AmpDiff, TimeDiff, weight); } - //} - //std::cout<<"hello11"<(regname, regionSelections_[i_reg]); - // } - //} + return; } diff --git a/processors/config/trackHitClusterAna.py b/processors/config/trackHitClusterAna.py new file mode 100644 index 000000000..40b50520d --- /dev/null +++ b/processors/config/trackHitClusterAna.py @@ -0,0 +1,57 @@ +import HpstrConf +import sys +import baseConfig as base +from baseConfig import bfield + +base.parser.add_argument("-L", "--layer", type=int, dest="layer", + help="Layer Under Investigation", metavar="layer", default=-1) +base.parser.add_argument("-M", "--module", type=int, dest="module", + help="Module Under Investigation", metavar="module", default=-1) +base.parser.add_argument("-MC", "--MC", type=int, dest="isMC", + help="Is the file used generated from Monte Carlo", metavar="module", default=0) +base.parser.add_argument("-doT", "--doT", type=int, dest="doTrack", + help="We plot the tracking based cluster performance metrics", metavar="doTrack", default=0) +base.parser.add_argument("-cut", "--cut", type=float, dest="cut", + help="Momentum Cut for NShared Profile", metavar="cut", default=-1.0) + + +options = base.parser.parse_args() + +# Use the input file to set the output file name +root_file = options.inFilename +ana_file = options.outFilename + +#print('LCIO file: %s' % root_file) +#print('Root file: %s' % ana_file) + +p = HpstrConf.Process() + +#p.max_events = 1000 +p.run_mode = 1 +p.skip_events = options.skip_events +p.max_events = options.nevents + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### +clua = HpstrConf.Processor('clua', 'SVTClusterAnaProcessor') +clua.parameters["debug"] = 0 +clua.parameters["layer"] = options.layer +clua.parameters["module"] = options.module +clua.parameters["isMC"] = options.isMC +clua.parameters["doTrack"] = options.doTrack +clua.parameters["cut"] = options.cut +clua.parameters["badchannels"] = os.environ['HPSTR_BASE']+"/analysis/data/badchannels2021.dat" + +sequence = [clua] + +p.sequence = sequence + +print("processors::clusterAna: The input file is: " + str(root_file)) +p.input_files = root_file +p.output_files = [ana_file] + +p.printProcess() diff --git a/processors/config/trackHitCompareAna.py b/processors/config/trackHitCompareAna.py new file mode 100644 index 000000000..b9f4e058c --- /dev/null +++ b/processors/config/trackHitCompareAna.py @@ -0,0 +1,57 @@ +import HpstrConf +import sys +import baseConfig as base +from baseConfig import bfield + +base.parser.add_argument("-L", "--layer", type=int, dest="layer", + help="Layer Under Investigation", metavar="layer", default=-1) +base.parser.add_argument("-M", "--module", type=int, dest="module", + help="Module Under Investigation", metavar="module", default=-1) +base.parser.add_argument("-MC", "--MC", type=int, dest="isMC", + help="Is the file used generated from Monte Carlo", metavar="module", default=0) +base.parser.add_argument("-doT", "--doT", type=int, dest="doTrack", + help="we plot tracking related cluster performance metrics", metavar="doTrack", default=0) +base.parser.add_argument("-cut", "--cut", type=float, dest="cut", + help="Momentum Cut for NShared Profile", metavar="cut", default=-1.0) + + +options = base.parser.parse_args() + +# Use the input file to set the output file name +root_file = options.inFilename +ana_file = options.outFilename + +#print('LCIO file: %s' % root_file) +#print('Root file: %s' % ana_file) + +p = HpstrConf.Process() + +#p.max_events = 1000 +p.run_mode = 1 +p.skip_events = options.skip_events +p.max_events = options.nevents + +# Library containing processors +p.add_library("libprocessors") + +############################### +# Processors # +############################### +cclua = HpstrConf.Processor('cclua', 'TrackHitCompareAnaProcessor') +cclua.parameters["debug"] = 0 +cclua.parameters["layer"] = options.layer +cclua.parameters["module"] = options.module +cclua.parameters["isMC"] = options.isMC +cclua.parameters["doTrack"] = options.doTrack +cclua.parameters["cut"] = options.cut +cclua.parameters["badchannels"] = os.environ['HPSTR_BASE']+"/analysis/data/badchannels2021.dat" + +sequence = [cclua] + +p.sequence = sequence + +print("processors::clusterAna: The input file is: " + str(root_file)) +p.input_files = root_file +p.output_files = [ana_file] + +p.printProcess() diff --git a/processors/include/SVTClusterAnaProcessor.h b/processors/include/SVTClusterAnaProcessor.h new file mode 100644 index 000000000..0d0ecda30 --- /dev/null +++ b/processors/include/SVTClusterAnaProcessor.h @@ -0,0 +1,175 @@ +#ifndef __SVTCLUSTER_ANAPROCESSOR_H__ +#define __SVTCLUSTER_ANAPROCESSOR_H__ + +// HPSTR +#include "HpsEvent.h" +#include "RawSvtHit.h" +#include "TrackerHit.h" +#include "RawSvtHitHistos.h" +#include "AnaHelpers.h" +#include "Event.h" +#include "BaseSelector.h" +#include "RawSvtHitHistos.h" +#include "EventHeader.h" +#include "VTPData.h" +#include "TSData.h" +#include "CalCluster.h" +#include "Track.h" +#include "TrackerHit.h" +#include "Collections.h" + +// ROOT +#include "Processor.h" +#include "TClonesArray.h" +#include "TBranch.h" +#include "TTree.h" +#include "TFile.h" +#include "TF1.h" +#include "TGraphErrors.h" +#include "TAxis.h" +#include "TROOT.h" +#include "TPad.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TProfile.h" + +class TTree; + + +class SVTClusterAnaProcessor : public Processor { + + public: + + SVTClusterAnaProcessor(const std::string& name, Process& process); + + ~SVTClusterAnaProcessor(); + + /** + * + * This method fills several histograms for each event directed at evaluating clustering performance. For our tracking + * variables we plot z0 vs the number of shared hits with and without a cluster cut (we expect a good reconstruction to have + * have a low number of shared clusters, especially around the z0 vertex and that this would further improve with high p tracks) + * + * To further characterize the clusters we capture basic information: cluster amplitude, time, charge, and strip cluster position + * and distances. We do this for shared and unshared clusters. We perform this for regular or NTD (next to dead) clusters. + * + **/ + virtual bool process(IEvent* ievent); + + virtual void initialize(TTree* tree); + + /** + * + * The remaining methods (up to TrackPlot) all plot standard histograms described in the description to the process method. + * They include standard function calls you would expect in a ROOT macro and afford some more control than a HistoManager and + * RegionSelector combo. + * Feasible for a limited number of collection cuts. + * + **/ + + virtual void PlotClusterLayers(); + + virtual void PlotClusterLayersNTD(); + + virtual void PlotClusterCharges(); + + virtual void PlotClusterChargesNTD(); + + virtual void PlotClusterPositions(); + + virtual void PlotClusterTimes(); + + virtual void PlotClusterTimesNTD(); + + virtual void TrackPlot(); + + /** + * + * This method fills the collection of dead channels given an input file name. Required for NTD plots, but largely disused in this processor as of yet. + * + **/ + virtual void fillDeads(); + + /** + * + * This method gets the strip count out of roughly 25000 given the FEB and hybrid IDs and the strip number with respect to these IDs. + * + **/ + virtual int GetStrip(int feb, int hyb,int strip); + + /** + * + * The finalize method calls all the plotting macros above. They are created into PNGs into the repository the processor is called in. + * + **/ + virtual void finalize(); + + virtual void configure(const ParameterSet& parameters); + + private: + + //Containers to hold histogrammer info + //RawSvtHitHistos* histos{nullptr}; + ModuleMapper * mmapper_; + + TTree* tree_; + TBranch* bClusters_{nullptr}; + TBranch* bClustersKF_{nullptr}; + TBranch* bsvtraw_{nullptr}; + TBranch* btracks_{nullptr}; + + int layer_{-1}; + int module_{-1}; + + TH1F* layers_; + TH1F* layersOnTrk_; + TH1F* layersOffTrk_; + TH1F* charges_; + TH1F* chargesOnTrk_; + TH1F* chargesOffTrk_; + + TH1F* layersNTD_; + TH1F* layersOnTrkNTD_; + TH1F* layersOffTrkNTD_; + TH1F* chargesNTD_; + TH1F* chargesOnTrkNTD_; + TH1F* chargesOffTrkNTD_; + + TH1F* positions_; + TH1F* positionsOnTrk_; + TH1F* ClusDistances_; + TH1F* ClusDistancesNTD_; + + TH1F* times_; + TH1F* timesOnTrk_; + TH1F* timesOffTrk_; + TH1F* timesNTD_; + TH1F* timesOnTrkNTD_; + TH1F* timesOffTrkNTD_; + + // TRACKING RELATED VARIABLES + + TH2F* Z0VNShare2Hist_; + TH2F* Z0VNShare2HistCut_; + TH1F* SharedAmplitudes_; + TH1F* UnSharedAmplitudes_; + TH1F* SharedTimes_; + TH1F* UnSharedTimes_; + + float pcut_{-1.0}; + + bool doingTracks_{false}; + + float Deads_[24576]; + + std::vector * Clusters_{}; + std::vector * ClustersKF_{}; + std::vector * svtraw_{}; + std::vector * tracks_{}; + + int debug_{0}; + int isMC_{0}; + std::string badchann_{""}; +}; + +#endif diff --git a/processors/include/SvtRawDataAnaProcessor.h b/processors/include/SvtRawDataAnaProcessor.h index 00616e717..47520175b 100644 --- a/processors/include/SvtRawDataAnaProcessor.h +++ b/processors/include/SvtRawDataAnaProcessor.h @@ -1,7 +1,7 @@ #ifndef __RAWSVTHIT_ANAPROCESSOR_H__ #define __RAWSVTHIT_ANAPROCESSOR_H__ -//HPSTR +// HPSTR #include "HpsEvent.h" #include "RawSvtHit.h" #include "RawSvtHitHistos.h" @@ -16,9 +16,7 @@ #include "Track.h" #include "TrackerHit.h" -//#include " -//ROOT - +// ROOT #include "Processor.h" #include "TClonesArray.h" #include "TBranch.h" @@ -43,27 +41,54 @@ class SvtRawDataAnaProcessor : public Processor { ~SvtRawDataAnaProcessor(); + /** + * + * Runs over the region selectors and checks if an event passes a selection json and fills a rawsvthithisto + * if do sample is on, it runs sampling. + * + **/ virtual bool process(IEvent* ievent); + /** + * + * Process initializer. Reads in the offline baselines into local baseline files, reads in the pulse shapes, and finally + * establishes regions which are used along with the region selector class and cuts in analysis/selection/svt to select on + * events for which histograms in rawsvthisto is filled. + * + **/ virtual void initialize(TTree* tree); virtual void sample(RawSvtHit* thisHit, std::string word, IEvent* ievent, long t,int i); + /** + * + * Four pole pulse function and the sum of two of them with baselines borrowed from Alic + * + **/ virtual TF1* fourPoleFitFunction(std::string word, int caser); + /** + * + * This method is implemented because C++ std:of method which converts strings + * to floats is not working. We need this to read in offline baselines and characteristic times. + * + **/ virtual float str_to_float(std::string word); float reverseEngineerTime(float ti, long t); - //virtual int maximum(int arr[]); - + /** + * + * Fills in histograms + * + **/ virtual void finalize(); virtual void configure(const ParameterSet& parameters); private: - //Containers to hold histogrammer info + // Containers to hold histogrammer info RawSvtHitHistos* histos{nullptr}; std::string histCfgFilename_; Float_t TimeRef_; @@ -86,7 +111,7 @@ class SvtRawDataAnaProcessor : public Processor { TBranch* brecoClu_{nullptr}; TBranch* bPart_{nullptr}; TBranch* bTrk_{nullptr}; - + TBranch* bClusters_{nullptr}; TBranch* bevH_; std::vector * svtHits_{}; @@ -95,7 +120,8 @@ class SvtRawDataAnaProcessor : public Processor { std::vector* recoClu_{}; std::vector* Trk_{}; std::vector* Part_{}; - //std::vector Trk_{}; + std::vector* Clusters_{}; + EventHeader * evH_; std::string anaName_{"rawSvtHitAna"}; @@ -109,7 +135,7 @@ class SvtRawDataAnaProcessor : public Processor { std::string timeProfiles_; int tphase_{6}; - //Debug Level + // Debug Level int debug_{0}; }; diff --git a/processors/include/TrackHitCompareAnaProcessor.h b/processors/include/TrackHitCompareAnaProcessor.h new file mode 100644 index 000000000..5caaa4db6 --- /dev/null +++ b/processors/include/TrackHitCompareAnaProcessor.h @@ -0,0 +1,223 @@ +#ifndef __TRACKHITCOMPARE_ANAPROCESSOR_H__ +#define __TRACKHITCOMPARE_ANAPROCESSOR_H__ + +// HPSTR +#include "HpsEvent.h" +#include "RawSvtHit.h" +#include "TrackerHit.h" +#include "RawSvtHitHistos.h" +#include "AnaHelpers.h" +#include "Event.h" +#include "BaseSelector.h" +#include "RawSvtHitHistos.h" +#include "EventHeader.h" +#include "VTPData.h" +#include "TSData.h" +#include "CalCluster.h" +#include "Track.h" +#include "TrackerHit.h" +#include "Collections.h" + +// ROOT +#include "Processor.h" +#include "TClonesArray.h" +#include "TBranch.h" +#include "TTree.h" +#include "TFile.h" +#include "TF1.h" +#include "TGraphErrors.h" +#include "TAxis.h" +#include "TROOT.h" +#include "TPad.h" +#include "TCanvas.h" +#include "TLegend.h" +#include "TProfile.h" + +class TTree; + + +class TrackHitCompareAnaProcessor : public Processor { + + public: + + TrackHitCompareAnaProcessor(const std::string& name, Process& process); + + ~TrackHitCompareAnaProcessor(); + + /** + * + * This method fills several histograms for each event directed at evaluating clustering performance. For our tracking + * variables we plot z0 vs the number of shared hits with and without a cluster cut (we expect a good reconstruction to have + * have a low number of shared clusters, especially around the z0 vertex and that this would further improve with high p tracks) + * We also directly compare the regular and transverse momentum distributions for the two methods. + * + * To further characterize the clusters we capture basic information: cluster amplitude, time, charge, and strip cluster position + * and distances. We do this for shared and unshared clusters. We perform this for regular or NTD (next to dead) clusters. + * + * This processor is undeniably extremely similar to ClusterAnaProcessor, having been its offspring. The main distinction is that + * this processor implements an identity tracker (ident_) which, atm, is implemented to detect whether an event is in type 1 or type 2 + * along with a short macro which associates to differently reconstructed root files a identification monicer, this allows for quick + * direct comparison of two means of reconstruction. + * + **/ + virtual bool process(IEvent* ievent); + + virtual void initialize(TTree* tree); + + /** + * + * The remaining methods (up to TrackPlot) all plot standard histograms described in the description to the process method. + * They include standard function calls you would expect in a root macro and afford some more control than a histomanager and + * region selector combo. + * + **/ + + virtual void PlotClusterLayers(); + + virtual void PlotClusterCharges(); + + virtual void PlotClusterTimes(); + + virtual void PlotClusterPositions(); + + virtual void TrackMomenta(); + + virtual void TrackTransverseMomenta(); + + /** + * + * This method fills the collection of dead channel ids given an input filename. Required for NTD plots, but largely disused in this processor as of yet. + * + **/ + virtual void fillDeads(); + + /** + * + * This method gets the strip count out of roughly 25000 given the feb and hybrid ids and the strip no wrt these IDs + * + **/ + virtual int GetStrip(int feb, int hyb,int strip); + + /** + * + * The finalize method calls all the plotting macros above. They are created into pngs into the repository the processor is called in. + * + **/ + virtual void finalize(); + + virtual void configure(const ParameterSet& parameters); + + private: + + // Containers to hold histogrammer info + ModuleMapper * mmapper_; + + TTree* tree_; + TBranch* bClusters_{nullptr}; + TBranch* bClustersKF_{nullptr}; + TBranch* bsvtraw_{nullptr}; + TBranch* btracks_{nullptr}; + TBranch* bident_{nullptr}; + + int layer_{-1}; + int module_{-1}; + + // FOR THE FIRST FILE + TH1F* layers1_; + TH1F* layersOnTrk1_; + TH1F* layersOffTrk1_; + TH1F* charges1_; + TH1F* chargesOnTrk1_; + TH1F* chargesOffTrk1_; + + TH1F* layersNTD1_; + TH1F* layersOnTrkNTD1_; + TH1F* layersOffTrkNTD1_; + TH1F* chargesNTD1_; + TH1F* chargesOnTrkNTD1_; + TH1F* chargesOffTrkNTD1_; + + TH1F* positions1_; + TH1F* positionsOnTrk1_; + TH1F* ClusDistances1_; + TH1F* ClusDistancesNTD1_; + + TH1F* times1_; + TH1F* timesOnTrk1_; + TH1F* timesOffTrk1_; + TH1F* timesNTD1_; + TH1F* timesOnTrkNTD1_; + TH1F* timesOffTrkNTD1_; + + // FOR THE SECOND FILE + TH1F* layers2_; + TH1F* layersOnTrk2_; + TH1F* layersOffTrk2_; + TH1F* charges2_; + TH1F* chargesOnTrk2_; + TH1F* chargesOffTrk2_; + + TH1F* layersNTD2_; + TH1F* layersOnTrkNTD2_; + TH1F* layersOffTrkNTD2_; + TH1F* chargesNTD2_; + TH1F* chargesOnTrkNTD2_; + TH1F* chargesOffTrkNTD2_; + + TH1F* positions2_; + TH1F* positionsOnTrk2_; + TH1F* ClusDistances2_; + TH1F* ClusDistancesNTD2_; + + TH1F* times2_; + TH1F* timesOnTrk2_; + TH1F* timesOffTrk2_; + TH1F* timesNTD2_; + TH1F* timesOnTrkNTD2_; + TH1F* timesOffTrkNTD2_; + + // TRACKING RELATED VARIABLES + TH2F* Z0VNShare2Hist1_; + TH2F* Z0VNShare2HistCut1_; + TH1F* SharedAmplitudes1_; + TH1F* UnSharedAmplitudes1_; + TH1F* SharedTimes1_; + TH1F* UnSharedTimes1_; + TH1F* TrackMomentumInTime1_; + TH1F* TrackMomentumOutTime1_; + TH1F* TrackMomentumAllTime1_; + TH1F* TrackMomentumTInTime1_; + TH1F* TrackMomentumTOutTime1_; + TH1F* TrackMomentumTAllTime1_; + + TH2F* Z0VNShare2Hist2_; + TH2F* Z0VNShare2HistCut2_; + TH1F* SharedAmplitudes2_; + TH1F* UnSharedAmplitudes2_; + TH1F* SharedTimes2_; + TH1F* UnSharedTimes2_; + TH1F* TrackMomentumInTime2_; + TH1F* TrackMomentumOutTime2_; + TH1F* TrackMomentumAllTime2_; + TH1F* TrackMomentumTInTime2_; + TH1F* TrackMomentumTOutTime2_; + TH1F* TrackMomentumTAllTime2_; + + float pcut_{-1.0}; + + bool doingTracks_{false}; + + float Deads_[24576]; + + std::vector * Clusters_{}; + std::vector * ClustersKF_{}; + std::vector * svtraw_{}; + std::vector * tracks_{}; + + int debug_{0}; + int isMC_{0}; + float ident_{1.0}; + std::string badchann_{""}; +}; + +#endif diff --git a/processors/src/SVTClusterAnaProcessor.cxx b/processors/src/SVTClusterAnaProcessor.cxx new file mode 100644 index 000000000..1e6db1a96 --- /dev/null +++ b/processors/src/SVTClusterAnaProcessor.cxx @@ -0,0 +1,629 @@ +/** + * @file SVTClusterAnaProcessor.cxx + * @brief AnaProcessor used fill histograms to study cluster reconstruction algorithms and dead channels. + * Does not feature a region selector or histomanager (i.e. no configurable json files), rather + * for the limited featured plots allows for more specific manipulation and control. + * @author Rory O'Dwyer and Cameron Bravo, SLAC National Accelerator Laboratory + */ +#include "SVTClusterAnaProcessor.h" +#include + +SVTClusterAnaProcessor::SVTClusterAnaProcessor(const std::string& name, Process& process) : Processor(name,process){ + mmapper_ = new ModuleMapper(2021); +} +SVTClusterAnaProcessor::~SVTClusterAnaProcessor(){} + +void SVTClusterAnaProcessor::configure(const ParameterSet& parameters) { + std::cout << "Configuring SVTClusterAnaProcessor" << std::endl; + try + { + debug_ = parameters.getInteger("debug"); + layer_ = parameters.getInteger("layer"); + module_ = parameters.getInteger("module"); + isMC_ = parameters.getInteger("isMC"); + doingTracks_ = (parameters.getInteger("doTrack")==1); + pcut_ = (float)parameters.getDouble("cut"); + badchann_ = parameters.getString("badchannels"); + } + catch (std::runtime_error& error) + { + std::cout << error.what() << std::endl; + } + +} + +void SVTClusterAnaProcessor::initialize(TTree* tree) { + fillDeads(); + tree_ = tree; + + if (isMC_==1){ + layers_ = new TH1F("layers", "MC Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrk_ = new TH1F("layersOnTrk", "MC Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrk_ = new TH1F("layersOffTrk", "MC Strip Width for Clusters off Track", 12, 0.0, 12.0); + charges_ = new TH1F("charges", "MC Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrk_ = new TH1F("chargesOnTrk", "MC Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrk_ = new TH1F("chargesOffTrk", "MC Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + layersNTD_ = new TH1F("layersNTD", "MC Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrkNTD_ = new TH1F("layersOnTrkNTD", "MC Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrkNTD_ = new TH1F("layersOffTrkNTD", "MC Strip Width for Clusters off Track", 12, 0.0, 12.0); + chargesNTD_ = new TH1F("chargesNTD", "MC Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrkNTD_ = new TH1F("chargesOnTrkNTD", "MC Charge Distribution for On Track", 1000, 0.0, .000016); + chargesOffTrkNTD_ = new TH1F("chargesOffTrkNTD", "MC Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + positions_ = new TH1F("Positions", "MC Location of Cluster Hit;Layer;Hits", 14, 0.0, 14.0); + positionsOnTrk_ = new TH1F("PositionsOnTrk", "MC Location of Cluster Hit for On Track", 14, 0.0, 14.0); + ClusDistances_ = new TH1F("Minimum Cluster Difference", "MC Minimum Distance Between Clusters", 14, 0.0, 14.0); + ClusDistancesNTD_ = new TH1F("Minimum Cluster Difference", "MC Minimum Distance Between Clusters", 14, 0.0, 14.0); + + times_ = new TH1F("Times", "MC Time of Cluster Hit", 1000, -60.0, 60.0); + timesOnTrk_ = new TH1F("TimesOnTrk", "MC Time of On Track Cluster Hit", 1000, -60.0, 60.0); + timesOffTrk_ = new TH1F("TimesOffTrk", "MC Time of Off Cluster Hit", 1000, -60.0, 60.0); + timesNTD_ = new TH1F("TimesNTD", "MC Time of Cluster Hit NTD", 1000, -60.0, 60.0); + timesOnTrkNTD_ = new TH1F("TimesOnTrkNTD", "MC Time of On Track Cluster Hit NTD", 1000, -60.0, 60.0); + timesOffTrkNTD_ = new TH1F("TimesOffTrkNTD", "MC Time of Off Cluster Hit NTD", 1000, -60.0, 60.0); + + tree_->SetBranchAddress("SiClusters", &Clusters_, &bClusters_); + tree_->SetBranchAddress("SiClustersOnTrack_KF", &ClustersKF_, &bClustersKF_); + tree_->SetBranchAddress("SVTRawTrackerHits", &svtraw_, &bsvtraw_); + return; + } + layers_ = new TH1F("layers", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrk_ = new TH1F("layersOnTrk", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrk_ = new TH1F("layersOffTrk", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + charges_ = new TH1F("charges", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrk_ = new TH1F("chargesOnTrk", "Charge Distribution for On Track", 1000, 0.0, .000016); + chargesOffTrk_ = new TH1F("chargesOffTrk", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + layersNTD_ = new TH1F("layersNTD", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrkNTD_ = new TH1F("layersOnTrkNTD", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrkNTD_ = new TH1F("layersOffTrkNTD", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + chargesNTD_ = new TH1F("chargesNTD", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrkNTD_ = new TH1F("chargesOnTrkNTD", "Charge Distribution for On Track", 1000, 0.0, .000016); + chargesOffTrkNTD_ = new TH1F("chargesOffTrkNTD", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + positions_ = new TH1F("Positions", "Location of Cluster Hit;Layer;Hits", 14, 0.0, 14.0); + positionsOnTrk_ = new TH1F("PositionsOnTrk", "Location of Cluster Hit for On Track", 14, 0.0, 14.0); + ClusDistances_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + ClusDistancesNTD_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + + times_ = new TH1F("Times", "Time of Cluster Hit", 1000, -60.0, 60.0); + timesOnTrk_ = new TH1F("TimesOnTrk", "Time of On Track Cluster Hit", 1000, -60.0, 60.0); + timesOffTrk_ = new TH1F("TimesOffTrk", "Time of Off Cluster Hit", 1000, -60.0, 60.0); + timesNTD_ = new TH1F("TimesNTD", "Time of Cluster Hit NTD", 1000, -60.0, 60.0); + timesOnTrkNTD_ = new TH1F("TimesOnTrkNTD", "Time of On Track Cluster Hit NTD", 1000, -60.0, 60.0); + timesOffTrkNTD_ = new TH1F("TimesOffTrkNTD", "Time of Off Cluster Hit NTD", 1000, -60.0, 60.0); + if (doingTracks_){ + Z0VNShare2Hist_ = new TH2F("Z0VNShare2Hist", "Z0 versus Number of Shared Hits No Cut", 100, 0, 3, 8, 0, 8); + Z0VNShare2HistCut_ = new TH2F("Z0VNShare2HistCut", "Z0 versus Number of Shared Hits Momentum Cut", 100, 0, 3, 8, 0, 8); + SharedAmplitudes_ = new TH1F("SharedAmplitudes", "The Amplitudes of Clusters Shared Between Tracks", 1000, 0.0, 0.000016); + UnSharedAmplitudes_ = new TH1F("UnSharedAmplitudes", "The Amplitudes of Clusters Not Shared Between Tracks", 1000, 0.0, 0.000016); + SharedTimes_ = new TH1F("SharedTimes", "The Times of Clusters Shared Between Tracks", 1000, -60.0, 60.0); + UnSharedTimes_ = new TH1F("UnSharedTimes", "The Times of Clusters Not Shared Between Tracks", 1000, -60.0, 60.0); + } + tree_->SetBranchAddress("SiClusters", &Clusters_, &bClusters_); + tree_->SetBranchAddress("SiClustersOnTrack_KF", &ClustersKF_, &bClustersKF_); + tree_->SetBranchAddress("SVTRawTrackerHits", &svtraw_, &bsvtraw_); + if (doingTracks_){ + tree_->SetBranchAddress("KalmanFullTracks", &tracks_, &btracks_); + } +} + +bool SVTClusterAnaProcessor::process(IEvent* ievent) { + if (doingTracks_){ + for (int i = 0; isize(); i++){ + Track* track = tracks_->at(i); + if (track->getTrackTime()*track->getTrackTime() < 100.0){ + Z0VNShare2Hist_->Fill(track->getZ0Err(), track->getNShared()); + if (track->getP() < pcut_){ + Z0VNShare2HistCut_->Fill(track->getZ0Err(), track->getNShared()); + } + } + } + } + for(int i = 0; i < Clusters_->size(); i++){ + TrackerHit * clu = Clusters_->at(i); + Int_t layc = -1; + Int_t modc = -1; + + RawSvtHit * seed = (RawSvtHit*)(clu->getRawHits().At(0)); + + layc = clu->getLayer(); + modc = seed->getModule(); + + if (doingTracks_){ + bool isShared = false; + int increment = 0; + for (int i = 0; isize(); i++){ + Track* track = tracks_->at(i); + for (int j = 0; jgetSvtHits().GetEntries(); j++){ + TrackerHit * test = (TrackerHit *)(track->getSvtHits().At(j)); + if (clu->getTime()==test->getTime()){ + increment += 1; + } + } + } + if (increment > 1){ + isShared = true; + } + if (increment > 0){ + bool general = ((layer_==-1)||(module_==-1)); + if (((layc==layer_)&&(modc==module_))||(general)){ + if (isShared){ + SharedAmplitudes_->Fill(clu->getCharge()); + SharedTimes_->Fill(clu->getTime()); + } else { + UnSharedAmplitudes_->Fill(clu->getCharge()); + UnSharedTimes_->Fill(clu->getTime()); + } + } + } + } + + float seedStrip = (float)(seed->getStrip()); + float nLayers = (float)(clu->getRawHits().GetEntries()); + float ncharges = (float)(clu->getCharge()); + float ntimes = (float)(clu->getTime()); + bool onTrk = false; + bool NTD = false; + for (unsigned int j = 0; j < ClustersKF_->size(); j++){ + if (clu->getID()==(ClustersKF_->at(j)->getID())){ + onTrk = true; + } + } + + std::string input = "ly" + std::to_string(layc+1) + "_m" + std::to_string(modc); + std::string helper = mmapper_->getHwFromSw(input); + + int feb = std::stoi(helper.substr(1, 1)); + int hyb = std::stoi(helper.substr(3, 1)); + + int channelL = seedStrip-1; + int channelR = seedStrip+1; + if (channelL >= 0){ + NTD = (NTD)||(Deads_[GetStrip(feb, hyb, channelL)]==1); + } + if (((feb<=1)&&(channelR<=511))||((feb>1)&&(channelR<=639))){ + NTD = (NTD)||(Deads_[GetStrip(feb, hyb, channelR)]==1); + } + bool general = ((layer_==-1)||(module_==-1)); + if (((layc==layer_)&&(modc==module_))||(general)){ + // Now is the part where I fill the cluster distance histogram. + float dist = 69420; + for(int p = 0; p < Clusters_->size(); p++){ + if (p==i) { continue; } + TrackerHit * clu2 = Clusters_->at(p); + RawSvtHit * seed2 = (RawSvtHit*)(clu2->getRawHits().At(0)); + float layc2 = clu->getLayer(); + float modc2 = seed->getModule(); + if ((!(layc2==layc))||(!(modc2==modc))) { continue; } + float new_dist = ((float)(seed2->getStrip())) - seedStrip; + if (new_dist < 0) { new_dist *= -1.0; } + if (new_dist < dist) { dist = new_dist; } + } + if (dist < 69420){ + ClusDistances_->Fill(dist); + if (NTD) {ClusDistancesNTD_->Fill(dist);} + } + layers_->Fill(nLayers); + charges_->Fill(ncharges); + positions_->Fill(clu->getLayer()); + times_->Fill(ntimes); + if (onTrk) { + layersOnTrk_->Fill(nLayers); + chargesOnTrk_->Fill(ncharges); + timesOnTrk_->Fill(ntimes); + } else { + layersOffTrk_->Fill(nLayers); + chargesOffTrk_->Fill(ncharges); + timesOffTrk_->Fill(ntimes); + } + if (NTD) { + layersNTD_->Fill(nLayers); + chargesNTD_->Fill(ncharges); + positionsOnTrk_->Fill(clu->getLayer()); + timesNTD_->Fill(ntimes); + if (onTrk) { + layersOnTrkNTD_->Fill(nLayers); + chargesOnTrkNTD_->Fill(ncharges); + timesOnTrkNTD_->Fill(ntimes); + } else { + layersOffTrkNTD_->Fill(nLayers); + chargesOffTrkNTD_->Fill(ncharges); + timesOffTrkNTD_->Fill(ntimes); + } + } + } + } + return true; +} + +void SVTClusterAnaProcessor::fillDeads(){ + for(int i = 0; i<24576; i++){ + Deads_[i] = 0.0; + } + std::string FILENAME=badchann_; + std::ifstream file(FILENAME.c_str()); + std::string line; + std::getline(file, line); + while (std::getline(file, line)) { + int value = std::atoi(line.c_str()); + Deads_[value] = 1.0; + } + file.close(); + return; +} + +int SVTClusterAnaProcessor::GetStrip(int feb, int hyb, int strip){ + int count = 0; + if (feb <= 1){ + count += feb*2048 + hyb*512 + strip; + } else { + count += 4096; + count += (feb-2)*2560 + hyb*640 + strip; + } + return count; +} + +void SVTClusterAnaProcessor::PlotClusterLayers(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + layers_->GetXaxis()->SetTitle("Width (strip)"); + layers_->GetYaxis()->SetTitle("Hits"); + layersOnTrk_->GetXaxis()->SetTitle("Width (strip)"); + layersOnTrk_->GetYaxis()->SetTitle("Hits"); + layersOffTrk_->GetXaxis()->SetTitle("Width (strip)"); + layersOffTrk_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Layers for All Clusters"); + layers_->Draw("e"); + c1->SaveAs("allClusters.png"); + c1->Clear(); + + c1->SetTitle("Layers for On Track Clusters"); + layersOnTrk_->Draw("e"); + c1->SaveAs("onClusters.png"); + c1->Clear(); + + c1->SetTitle("Layers for Off Track Clusters"); + layersOffTrk_->Draw("e"); + c1->SaveAs("offClusters.png"); + c1->Clear(); + + layers_->SetTitle("Cluster Strip Width for all Cluster Cuts"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + layers_->SetLineColor(kCyan); + layersOnTrk_->SetLineColor(kGreen); + legend->AddEntry(layers_, "Layers"); + legend->AddEntry(layersOnTrk_, "Layers On Track"); + legend->AddEntry(layersOffTrk_, "Layers Off Track"); + layers_->Draw("e"); + layersOnTrk_->Draw("same e"); + layersOffTrk_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogether.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::PlotClusterCharges(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + + charges_->GetXaxis()->SetTitle("Charge"); + charges_->GetYaxis()->SetTitle("Hits"); + chargesOffTrk_->GetXaxis()->SetTitle("Charge"); + chargesOffTrk_->GetYaxis()->SetTitle("Hits"); + chargesOnTrk_->GetXaxis()->SetTitle("Charge"); + chargesOnTrk_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Charges for All Clusters"); + charges_->Draw("e"); + c1->SaveAs("allClustersCharge.png"); + c1->Clear(); + + c1->SetTitle("Charges for On Track Clusters"); + chargesOnTrk_->Draw("e"); + c1->SaveAs("onClustersCharge.png"); + c1->Clear(); + + c1->SetTitle("Charges for Off Track Clusters"); + chargesOffTrk_->Draw("e"); + c1->SaveAs("offClustersCharge.png"); + c1->Clear(); + + layers_->SetTitle("Putting all Charges Together"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + charges_->SetLineColor(kCyan); + chargesOnTrk_->SetLineColor(kGreen); + legend->AddEntry(charges_, "Charges"); + legend->AddEntry(chargesOnTrk_, "Charges On Track"); + legend->AddEntry(chargesOffTrk_, "Charges Off Track"); + charges_->Draw("e"); + chargesOnTrk_->Draw("same e"); + chargesOffTrk_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogetherCharges.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::PlotClusterLayersNTD(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + layersNTD_->GetXaxis()->SetTitle("Width (strip)"); + layersNTD_->GetYaxis()->SetTitle("Hits"); + layersOnTrkNTD_->GetXaxis()->SetTitle("Width (strip)"); + layersOnTrkNTD_->GetYaxis()->SetTitle("Hits"); + layersOffTrkNTD_->GetXaxis()->SetTitle("Width (strip)"); + layersOffTrkNTD_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("NTD Layers for All Clusters"); + layersNTD_->Draw("e"); + c1->SaveAs("allClustersNTD.png"); + c1->Clear(); + + c1->SetTitle("NTD Layers for On Track Clusters"); + layersOnTrkNTD_->Draw("e"); + c1->SaveAs("onClustersNTD.png"); + c1->Clear(); + + c1->SetTitle("NTD Layers for Off Track Clusters"); + layersOffTrkNTD_->Draw("e"); + c1->SaveAs("offClustersNTD.png"); + c1->Clear(); + + layersNTD_->SetTitle("NTD Cluster Strip Width for all Cluster Cuts"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + layersNTD_->SetLineColor(kCyan); + layersOnTrkNTD_->SetLineColor(kGreen); + legend->AddEntry(layersNTD_, "Layers"); + legend->AddEntry(layersOnTrkNTD_, "Layers On Track"); + legend->AddEntry(layersOffTrkNTD_, "Layers Off Track"); + layersNTD_->Draw("e"); + layersOnTrkNTD_->Draw("same e"); + layersOffTrkNTD_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogetherNTD.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::PlotClusterChargesNTD(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + + chargesNTD_->GetXaxis()->SetTitle("Charge"); + chargesNTD_->GetYaxis()->SetTitle("Hits"); + chargesOffTrkNTD_->GetXaxis()->SetTitle("Charge"); + chargesOffTrkNTD_->GetYaxis()->SetTitle("Hits"); + chargesOnTrkNTD_->GetXaxis()->SetTitle("Charge"); + chargesOnTrkNTD_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("NTD Charges for All Clusters"); + chargesNTD_->Draw("e"); + c1->SaveAs("allClustersChargeNTD.png"); + c1->Clear(); + + c1->SetTitle("NTD Charges for On Track Clusters"); + chargesOnTrkNTD_->Draw("e"); + c1->SaveAs("onClustersChargeNTD.png"); + c1->Clear(); + + c1->SetTitle("NTD Charges for Off Track Clusters"); + chargesOffTrkNTD_->Draw("e"); + c1->SaveAs("offClustersChargeNTD.png"); + c1->Clear(); + + chargesNTD_->SetTitle("NTD Charge Distribution for all Charges"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + chargesNTD_->SetLineColor(kCyan); + chargesOnTrkNTD_->SetLineColor(kGreen); + legend->AddEntry(chargesNTD_, "Charges"); + legend->AddEntry(chargesOnTrkNTD_, "Charges On Track"); + legend->AddEntry(chargesOffTrkNTD_, "Charges Off Track"); + chargesNTD_->Draw("e"); + chargesOnTrkNTD_->Draw("same e"); + chargesOffTrkNTD_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogetherChargesNTD.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::PlotClusterPositions(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + positions_->GetXaxis()->SetTitle("Layer"); + positions_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Cluster Position for all Clusters"); + positions_->Draw("e"); + c1->SaveAs("positions.png"); + c1->Clear(); + + positionsOnTrk_->GetXaxis()->SetTitle("Layer"); + positionsOnTrk_->GetYaxis()->SetTitle("Hits"); + + c1->SetTitle("Cluster Position for NTD Clusters"); + positionsOnTrk_->Draw("e"); + c1->SaveAs("positionsNTD.png"); + c1->Clear(); + + ClusDistances_->GetXaxis()->SetTitle("Separation"); + ClusDistances_->GetYaxis()->SetTitle("Pairs"); + + c1->SetTitle("Cluster Separation on Same Sensors"); + ClusDistances_->Draw("e"); + c1->SaveAs("clusDistances.png"); + c1->Clear(); + + ClusDistancesNTD_->GetXaxis()->SetTitle("Separation"); + ClusDistancesNTD_->GetYaxis()->SetTitle("Pairs"); + + c1->SetTitle("Cluster Separation on Same Sensors for NTD"); + ClusDistancesNTD_->Draw("e"); + c1->SaveAs("clusDistancesNTD.png"); + c1->Clear(); +} + +void SVTClusterAnaProcessor::PlotClusterTimes(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + times_->GetXaxis()->SetTitle("Times (ns)"); + times_->GetYaxis()->SetTitle("Hits"); + timesOnTrk_->GetXaxis()->SetTitle("Time (ns)"); + timesOnTrk_->GetYaxis()->SetTitle("Hits"); + timesOffTrk_->GetXaxis()->SetTitle("Time (ns)"); + timesOffTrk_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Times for All Clusters"); + times_->Draw(); + c1->SaveAs("alltimes.png"); + c1->Clear(); + + c1->SetTitle("Times for On Track Clusters"); + timesOnTrk_->Draw("e"); + c1->SaveAs("ontimes.png"); + c1->Clear(); + + c1->SetTitle("Times for Off Track Clusters"); + timesOffTrk_->Draw("e"); + c1->SaveAs("offtimes.png"); + c1->Clear(); + + times_->SetTitle("Cluster Times for all Cluster Cuts"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + times_->SetLineColor(kCyan); + timesOnTrk_->SetLineColor(kGreen); + legend->AddEntry(times_, "Times"); + legend->AddEntry(timesOnTrk_, "Times On Track"); + legend->AddEntry(timesOffTrk_, "Times Off Track"); + times_->Draw("e"); + timesOnTrk_->Draw("same e"); + timesOffTrk_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogetherTimes.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::PlotClusterTimesNTD(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + timesNTD_->GetXaxis()->SetTitle("Times (ns)"); + timesNTD_->GetYaxis()->SetTitle("Hits"); + timesOnTrkNTD_->GetXaxis()->SetTitle("Time (ns)"); + timesOnTrkNTD_->GetYaxis()->SetTitle("Hits"); + timesOffTrkNTD_->GetXaxis()->SetTitle("Time (ns)"); + timesOffTrkNTD_->GetYaxis()->SetTitle("Hits"); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Times for All Clusters"); + timesNTD_->Draw("e"); + c1->SaveAs("alltimesNTD.png"); + c1->Clear(); + + c1->SetTitle("Times for On Track Clusters"); + timesOnTrkNTD_->Draw("e"); + c1->SaveAs("ontimesNTD.png"); + c1->Clear(); + + c1->SetTitle("Times for Off Track Clusters"); + timesOffTrkNTD_->Draw("e"); + c1->SaveAs("offtimesNTD.png"); + c1->Clear(); + + timesNTD_->SetTitle("Cluster Times for all Cluster Cuts"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + timesNTD_->SetLineColor(kCyan); + timesOnTrkNTD_->SetLineColor(kGreen); + legend->AddEntry(timesNTD_, "Times"); + legend->AddEntry(timesOnTrkNTD_, "Times On Track"); + legend->AddEntry(timesOffTrkNTD_, "Times Off Track"); + timesNTD_->Draw("e"); + timesOnTrkNTD_->Draw("same e"); + timesOffTrkNTD_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("PutOnTogetherTimesNTD.png"); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::TrackPlot(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + + // First, I do the profiles of the shared hits plots + TProfile* p1 = Z0VNShare2Hist_->ProfileX("p1"); + TProfile* p2 = Z0VNShare2HistCut_->ProfileX("p2"); + auto legend = new TLegend(0.3, 0.8, .68, .9); + p2->SetLineColor(kRed); + p2->SetLineWidth(2.0); + p1->SetLineWidth(2.0); + legend->AddEntry(p1, "Shared Profile Not Cut"); + legend->AddEntry(p2, "Shared Profile Cut"); + p1->Draw("e"); + p2->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("Z0VSharedProfile.png"); + c1->Clear(); + + // Now I do the charge distribution for shared and unshared hits overlayed. + legend = new TLegend(0.3, 0.8, .68, .9); + SharedAmplitudes_->SetLineColor(kRed); + SharedAmplitudes_->SetLineWidth(2.0); + UnSharedAmplitudes_->SetLineWidth(2.0); + legend->AddEntry(SharedAmplitudes_, "Charge Distribution for Shared Clusters"); + legend->AddEntry(UnSharedAmplitudes_, "Charge Distribution for UnShared Clusters"); + UnSharedAmplitudes_->Draw("e"); + SharedAmplitudes_->Draw("same e"); + legend->Draw("same e"); + gPad->SetLogy(true); + c1->SaveAs("SharedVUnSharedChargeDist.png"); + c1->Clear(); + + // Now I do the time distribution for shared and unshared hits overlayed. + SharedTimes_->SetTitle("The Time of Clusters Shared Between Tracks"); + legend = new TLegend(0.3, 0.8, .68, .9); + SharedTimes_->SetLineColor(kRed); + SharedTimes_->SetLineWidth(2.0); + UnSharedTimes_->SetLineWidth(2.0); + legend->AddEntry(SharedTimes_, "Time Distribution for Shared Clusters"); + legend->AddEntry(UnSharedTimes_, "Time Distribution for UnShared Clusters"); + UnSharedTimes_->Draw("e"); + SharedTimes_->Draw("same e"); + legend->Draw("same e"); + c1->SaveAs("SharedVUnSharedTimeDist.png"); + gPad->SetLogy(false); + c1->Clear(); + return; +} + +void SVTClusterAnaProcessor::finalize() { + PlotClusterLayers(); + PlotClusterLayersNTD(); + PlotClusterCharges(); + PlotClusterChargesNTD(); + PlotClusterPositions(); + PlotClusterTimes(); + PlotClusterTimesNTD(); + if(doingTracks_){ + TrackPlot(); + } + return; +} +DECLARE_PROCESSOR(SVTClusterAnaProcessor); diff --git a/processors/src/SvtRawDataAnaProcessor.cxx b/processors/src/SvtRawDataAnaProcessor.cxx index 5d4f9b640..b42e400bf 100644 --- a/processors/src/SvtRawDataAnaProcessor.cxx +++ b/processors/src/SvtRawDataAnaProcessor.cxx @@ -1,7 +1,7 @@ /** * @file SvtRawDataAnaProcessor.cxx * @brief AnaProcessor used fill histograms to study svt hit fitting - * @author Cameron Bravo, SLAC National Accelerator Laboratory + * @author Rory O'Dwyer and Cameron Bravo, SLAC National Accelerator Laboratory */ #include "SvtRawDataAnaProcessor.h" @@ -10,7 +10,6 @@ SvtRawDataAnaProcessor::SvtRawDataAnaProcessor(const std::string& name, Process& process) : Processor(name,process){ mmapper_ = new ModuleMapper(); } -//TODO CHECK THIS DESTRUCTOR SvtRawDataAnaProcessor::~SvtRawDataAnaProcessor(){} @@ -18,9 +17,9 @@ void SvtRawDataAnaProcessor::configure(const ParameterSet& parameters) { std::cout << "Configuring SvtRawDataAnaProcessor" << std::endl; try { - debug_ = parameters.getInteger("debug"); - anaName_ = parameters.getString("anaName"); - svtHitColl_ = parameters.getString("trkrHitColl"); + debug_ = parameters.getInteger("debug"); + anaName_ = parameters.getString("anaName"); + svtHitColl_ = parameters.getString("trkrHitColl"); histCfgFilename_ = parameters.getString("histCfg"); regionSelections_ = parameters.getVString("regionDefinitions"); TimeRef_ = parameters.getDouble("timeref"); @@ -38,35 +37,21 @@ void SvtRawDataAnaProcessor::configure(const ParameterSet& parameters) { } -/** - * - *THIS METHOD IS IMPLEMENTED BECAUSE C++ std:of METHOD WHICH CONVERTS STRINGS - TO FLOATS IS NOT WORKING. WE NEED THIS TO READ IN OFFLINE BASELINES AND CHARACTERISTIC TIMES. - * - * - * */ - - float SvtRawDataAnaProcessor::str_to_float(std::string token){ - std::string top1=token.substr(0,token.find(".")); - const char *top=top1.c_str(); - std::string bot1=token.substr(token.find(".")+1); - const char *bottom=bot1.c_str(); + std::string top1 = token.substr(0, token.find(".")); + const char *top = top1.c_str(); + std::string bot1 = token.substr(token.find(".") + 1); + const char *bottom = bot1.c_str(); float base = 0.0; - for(int J=0;J baselines; for(int i=0; i<24576; i++){ - std::getline(myfile,s); - std::getline(myfile2,s2); - int feb=0; - int hyb=0; - int ch=0; - if(i>=4096){ - feb=((i-4096)/2560); - hyb=(i-4096-feb*2560)/640; - ch=i-4096-feb*2560-hyb*640; - }else{ - feb=i/2048; - hyb=(i-feb*2048)/512; - ch=i-feb*2048-hyb*512; + std::getline(myfile, s); + std::getline(myfile2, s2); + int feb = 0; + int hyb = 0; + int ch = 0; + if (i>=4096){ + feb = ((i-4096)/2560); + hyb = (i-4096-feb*2560)/640; + ch = i-4096-feb*2560-hyb*640; + } else { + feb = i/2048; + hyb = (i-feb*2048)/512; + ch = i-feb*2048-hyb*512; } - for(int I=0;I<5;I++){ - std::string token=s2.substr(0,s2.find(",")); - s2=s2.substr(s2.find(",")+1); - if(I>=2){ - if(i<=4096){ - times1_[feb][hyb][ch][I-2]=str_to_float(token); - }else{ - times2_[feb][hyb][ch][I-2]=str_to_float(token); + for (int j=0; j<5; j++){ + std::string token = s2.substr(0,s2.find(",")); + s2 = s2.substr(s2.find(",")+1); + if (j>=2){ + if (i<=4096){ + times1_[feb][hyb][ch][j-2] = str_to_float(token); + } else { + times2_[feb][hyb][ch][j-2] = str_to_float(token); } } } - //if(i<2048){ - //std::cout<0){ - if(i<=4096){ - std::string token=s.substr(0,s.find(" ")); - if(debug_){ + for(int j=0; j<13; j++){ + if (j>0){ + if (i<=4096){ + std::string token = s.substr(0, s.find(" ")); + if (debug_){ std::cout<SetBranchAddress(svtHitColl_.c_str(), &svtHits_, &bsvtHits_); + tree_->SetBranchAddress("FinalStateParticles_KF", &Part_, &bPart_); + tree_->SetBranchAddress("SiClusters", &Clusters_, &bClusters_); + tree_->SetBranchAddress("EventHeader", &evH_, &bevH_); - tree_= tree; - // init histos - //histos = new RawSvtHitHistos(anaName_.c_str(), mmapper_); - //histos->loadHistoConfig(histCfgFilename_); - //histos->DefineHistos(); - //std::cout<size()<SetBranchAddress(svtHitColl_.c_str() , &svtHits_ , &bsvtHits_ ); - //tree_->SetBranchAddress("VTPBank", &vtpBank_ , &bvtpBank_ ); - //tree_->SetBranchAddress("TSBank", &tsBank_ , &btsBank_ ); - //tree_->SetBranchAddress("RecoEcalClusters",&recoClu_,&brecoClu_ ); - //tree_->SetBranchAddress("KalmanFullTracks",&Trk_,&bTrk_); - tree_->SetBranchAddress("FinalStateParticles_KF",&Part_,&bPart_); - tree_->SetBranchAddress("EventHeader",&evH_,&bevH_); - for (unsigned int i_reg = 0; i_reg < regionSelections_.size(); i_reg++) { std::string regname = AnaHelpers::getFileName(regionSelections_[i_reg],false); @@ -234,175 +169,193 @@ void SvtRawDataAnaProcessor::initialize(TTree* tree) { reg_histos_[regname] = std::make_shared(regname,mmapper_); reg_histos_[regname]->loadHistoConfig(histCfgFilename_); - //reg_histos_[regname]->doTrackComparisonPlots(false); reg_histos_[regname]->DefineHistos(); regions_.push_back(regname); } } -/* - * - *RUNS OVER THE REGION SELECTORS AND CHECKS IF AN EVENT PASSES A SELECTION JSON AND FILLS A RAWSVTHITHISTO - IF DO SAMPLE IS ON, IT RUNS SAMPLING. - * - * - */ - - bool SvtRawDataAnaProcessor::process(IEvent* ievent) { - Float_t TimeRef=-0.0; - Float_t AmpRef=1000.0; - double weight = 1.;int count1=0;int count2=0; + Float_t TimeRef = -0.0; + Float_t AmpRef = 1000.0; + double weight = 1.; + int count1 = 0; + int count2 = 0; long eventTime = evH_->getEventTime(); - - //I AM DOING CLUSTER MATCHING HERE :) - //std::cout<<"Here is the eventTime"<T<prescaled.Single_3_Top==1)or(tsBank_->prescaled.Single_3_Bot==1)))){return true;} - + int stripID = -10000; + int hitc = 0; + int hitl = 0; + float otherTime = 69420.0; //ONLY POSITRONS, MAY USE FEE's //ONCE I DETERMINE A CLUSTER WHICH IS IN LINE WITH TRIG, I CAN USE ANY CLUSTERS CLOSE IN TIME. - //std::cout<<"Trigger Time: "<singletrigs.at(0).T<size(); i++){ - RawSvtHit * thisHit = svtHits_->at(i); - int getNum = thisHit->getFitN();//std::cout<<"I got here 10"<size();i++){ - if(Part_->at(i)->getPDG()==22){continue;} - if(Part_->at(i)->getCluster().getEnergy()<0){continue;} - if(not((Part_->at(i)->getCluster().getTime()<=40)and(Part_->at(i)->getCluster().getTime()>=36))){continue;} - //std::cout<<"For each Tracker Hit I now print out Raw Hit Info: "<at(i)->getTrack().getSvtHits().GetEntries();j++){ + int trigPhase = (int)((eventTime%24)/4); + if ((trigPhase!=tphase_)&&(tphase_!=6)) { return true; } + + for (unsigned int i = 0; isize(); i++){ + Clusters_->at(i)->getLayer(); + } + + for (unsigned int I = 0; I < svtHits_->size(); I++){ + RawSvtHit * thisHit = svtHits_->at(I); + int getNum = thisHit->getFitN(); + if (doClMatch){ + bool continue_flag = true; + for (int i = 0; isize(); i++){ + if (Part_->at(i)->getPDG()==22) { continue; } + if (Part_->at(i)->getCluster().getEnergy() < 0) { continue; } + if (!((Part_->at(i)->getCluster().getTime()<=40)&&(Part_->at(i)->getCluster().getTime()>=36))) { continue; } + for (int j = 0; jat(i)->getTrack().getSvtHits().GetEntries(); j++){ TrackerHit * tHit = (TrackerHit*)(Part_->at(i)->getTrack().getSvtHits().At(j)); - //std::cout<getTime()<getRawHits().GetEntries();k++){ + double TrackTime = Part_->at(i)->getTrack().getTrackTime(); + for (int k = 0; kgetRawHits().GetEntries(); k++){ RawSvtHit * rHit = (RawSvtHit*)(tHit->getRawHits().At(k)); - if(rHit->getT0(0)==thisHit->getT0(0)){Continue=false;} - //std::cout<<"Raw Hit T0: "<getT0(0)<getT0(0)==thisHit->getT0(0))&&(mode==0)){//or(mode==2))){ + // This is the HIT ON TRACK Modes + bool inCluster = false; + int layc = 0; // THE PURPOSE OF layc IS TO COUNT THE NUMBER OF HITS PER LAYER + for (int cl = 0; cl < Clusters_->size(); cl++){ + for (int clh = 0; clh < Clusters_->at(cl)->getRawHits().GetEntries(); clh++){ + RawSvtHit * cluHit = (RawSvtHit*)(Clusters_->at(cl)->getRawHits().At(clh)); + if ((cluHit->getLayer()==thisHit->getLayer())&&(cluHit->getModule()==thisHit->getModule())){ + layc++; + } + if ((cluHit->getT0(0)==thisHit->getT0(0))){ // and(not(Clusters_->at(cl)->getID()==tHit->getID()))){ + inCluster = true; + hitc = Clusters_->at(cl)->getRawHits().GetEntries(); + hitl = layc; + if (Clusters_->at(cl)->getRawHits().GetEntries()==2){ + RawSvtHit * otherHit = (RawSvtHit*)(Clusters_->at(cl)->getRawHits().At((clh+1)%2)); + otherTime = otherHit->getT0(0); + } + } + } + } + if (inCluster){ + continue_flag = false; + } + } + + if ((rHit->getT0(0)<=-30)&&(not(rHit->getT0(0)==thisHit->getT0(0)))&&(mode==1)){ + // This is the HIT OFF TRACK Mode + // You are looking at the really early time on track hits, and specifically at other hits in the same layer and module + // to see if you have evidence of a misplaced hit. + + if ((rHit->getLayer()==thisHit->getLayer())&&(rHit->getModule()==thisHit->getModule())){ + stripID = rHit->getStrip(); + // This is conditioned on it being in clusters. + bool inCluster = false; + int layc = 0; // THE PURPOSE OF layc IS TO COUNT THE NUMBER OF HITS PER LAYER + for( int cl = 0; cl < Clusters_->size(); cl++){ + for (int clh = 0; clh < Clusters_->at(cl)->getRawHits().GetEntries(); clh++){ + RawSvtHit * cluHit = (RawSvtHit*)(Clusters_->at(cl)->getRawHits().At(clh)); + if ((cluHit->getLayer()==thisHit->getLayer())&&(cluHit->getModule()==thisHit->getModule())){ + layc++; + } + if ((cluHit->getT0(0)==thisHit->getT0(0))&&(!(Clusters_->at(cl)->getID()==tHit->getID()))){ + inCluster = true; + hitc = Clusters_->at(cl)->getRawHits().GetEntries(); + hitl = layc; + } + } + } + if(inCluster){ + continue_flag = false; + } + } + } } - //std::cout<<" T0: "<at(i)->getTrack().getSvtHits().At(j)->getRawHits().At(0).getT0()<passCutEq("getN_et",getNum,weight))){continue;} - - if(getNum==2){ - TimeDiff=(thisHit->getT0(J))-(thisHit->getT0((J+1)%2)); - AmpDiff=(thisHit->getT0(J))-(thisHit->getT0((J+1)%2)); - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("TimeDiff_lt",TimeDiff*TimeDiff,weight))){continue;} - } - if(!(reg_selectors_[regions_[i_reg]]->passCutEq("getId_lt",J,weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutEq("getId_gt",J,weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("chi_lt",thisHit->getChiSq(J),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutGt("chi_gt",thisHit->getChiSq(J),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutEq("getN_et", getNum,weight))) { continue; } + + if (getNum==2){ + TimeDiff = (thisHit->getT0(j))-(thisHit->getT0((j+1)%2)); + AmpDiff = (thisHit->getT0(j))-(thisHit->getT0((j+1)%2)); + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("TimeDiff_lt", TimeDiff*TimeDiff, weight))) { continue; } + } + + if (!(reg_selectors_[regions_[i_reg]]->passCutEq("getId_lt", j, weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutEq("getId_gt", j, weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("chi_lt", thisHit->getChiSq(j), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutGt("chi_gt", thisHit->getChiSq(j), weight))) { continue; } - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ft",(((thisHit->getT0(J))-TimeRef)*((thisHit->getT0(J))-TimeRef)<((thisHit->getT0((J+1)%getNum)-TimeRef)*(thisHit->getT0((J+1)%getNum)-TimeRef)+.00001)),weight))){continue;} - if(i_regpassCutLt("doing_ct",(((thisHit->getT0(J))-TimeRef)*((thisHit->getT0(J))-TimeRef)>((thisHit->getT0((J+1)%getNum)-TimeRef)*(thisHit->getT0((J+1)%getNum)-TimeRef)+.00001)),weight))){continue;} - }else{ - if(getNum==2){ - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ct",(((thisHit->getT0(J))-TimeRef)*((thisHit->getT0(J))-TimeRef)>((thisHit->getT0((J+1)%getNum)-TimeRef)*(thisHit->getT0((J+1)%getNum)-TimeRef)+.00001)),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ft", (((thisHit->getT0(j))-TimeRef)*((thisHit->getT0(j))-TimeRef) < ((thisHit->getT0((j+1)%getNum)-TimeRef)*(thisHit->getT0((j+1)%getNum)-TimeRef)+.00001)), weight))) { continue; } + if (i_reg < regionSelections_.size()-1){ + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ct", (((thisHit->getT0(j))-TimeRef)*((thisHit->getT0(j))-TimeRef) > ((thisHit->getT0((j+1)%getNum)-TimeRef)*(thisHit->getT0((j+1)%getNum)-TimeRef)+.00001)), weight))) { continue; } + } else { + if (getNum==2){ + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ct", (((thisHit->getT0(j))-TimeRef)*((thisHit->getT0(j))-TimeRef) > ((thisHit->getT0((j+1)%getNum)-TimeRef)*(thisHit->getT0((j+1)%getNum)-TimeRef)+.00001)), weight))) { continue; } } } - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ca",(((thisHit->getAmp(J))-AmpRef)*((thisHit->getAmp(J))-AmpRef)<((thisHit->getAmp((J+1)%getNum)-AmpRef)*(thisHit->getAmp((J+1)%getNum)-AmpRef)+.00001)),weight))){continue;} - - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_fterr",(((thisHit->getT0err(J))-TimeRef)*((thisHit->getT0err(J))-TimeRef)<((thisHit->getT0err((J+1)%getNum)-TimeRef)*(thisHit->getT0err((J+1)%getNum)-TimeRef)+.00001)),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_cterr",(((thisHit->getT0err(J))-0.0)*((thisHit->getT0err(J))-0.0)>((thisHit->getT0err((J+1)%getNum)-0.0)*(thisHit->getT0err((J+1)%getNum)-0.0)+.00001)),weight))){continue;} - - //if(!(std::abs((thisHit->getT0(J))-TimeRef)getT0((J+1)%2)-TimeRef))){continue;} - //if(!(reg_selectors_[regions_[i_reg]]->passCutEq("doing_ca",1.0,weight))){continue;}else{ - //if(!(std::abs((thisHit->getT0(J))-TimeRef)getT0((J+1)%2)-TimeRef))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("amp_lt",thisHit->getAmp(0),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutGt("amp_gt",thisHit->getAmp(0),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_ca", (((thisHit->getAmp(j))-AmpRef)*((thisHit->getAmp(j))-AmpRef) < ((thisHit->getAmp((j+1)%getNum)-AmpRef)*(thisHit->getAmp((j+1)%getNum)-AmpRef)+.00001)), weight))) { continue; } - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("time_lt",thisHit->getT0(0),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutGt("time_gt",thisHit->getT0(0),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_fterr", (((thisHit->getT0err(j))-TimeRef)*((thisHit->getT0err(j))-TimeRef) < ((thisHit->getT0err((j+1)%getNum)-TimeRef)*(thisHit->getT0err((j+1)%getNum)-TimeRef)+.00001)), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("doing_cterr", (((thisHit->getT0err(j))-0.0)*((thisHit->getT0err(j))-0.0) > ((thisHit->getT0err((j+1)%getNum)-0.0)*(thisHit->getT0err((j+1)%getNum)-0.0)+.00001)),weight))) { continue; } + + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("amp_lt", thisHit->getAmp(0), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutGt("amp_gt", thisHit->getAmp(0), weight))) { continue; } - //std::cout<<(float)(thisHit->getT0((J+1)%getNum))<passCutLt("Otime_lt",(float)(thisHit->getT0((J+1)%getNum)),weight))<passCutGt("Otime_gt",(float)(thisHit->getT0((J+1)%getNum)),weight))<passCutLt("Otime_lt",(float)(thisHit->getT0((J+1)%getNum)),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutGt("Otime_gt",(float)(thisHit->getT0((J+1)%getNum)),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("time_lt", thisHit->getT0(0), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutGt("time_gt", thisHit->getT0(0), weight))) { continue; } - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("Stime_lt",(float)(thisHit->getT0((J)%getNum)),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutGt("Stime_gt",(float)(thisHit->getT0((J)%getNum)),weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("Otime_lt", (float)(thisHit->getT0((j+1)%getNum)), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutGt("Otime_gt", (float)(thisHit->getT0((j+1)%getNum)), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("Stime_lt", (float)(thisHit->getT0((j)%getNum)), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutGt("Stime_gt", (float)(thisHit->getT0((j)%getNum)), weight))) { continue; } - if(!(reg_selectors_[regions_[i_reg]]->passCutLt("amp2_lt",thisHit->getAmp(0),weight))){continue;} - if(!(reg_selectors_[regions_[i_reg]]->passCutEq("channel", (thisHit->getStrip()),weight))){continue;} - int * adcs=thisHit->getADCs(); + if (!(reg_selectors_[regions_[i_reg]]->passCutLt("amp2_lt", thisHit->getAmp(0), weight))) { continue; } + if (!(reg_selectors_[regions_[i_reg]]->passCutEq("channel", (thisHit->getStrip()), weight))) { continue; } + int * adcs = thisHit->getADCs(); int maxx = 0; - for(unsigned int K=0; K<6; K++){ - if(maxxpassCutEq("first_max",adcs[0]-maxx,weight))){continue;} - - //if(!(reg_selectors_[regions_[i_reg]]->passCutEq("time_phase",trigPhase,weight))){continue;} + if (!(reg_selectors_[regions_[i_reg]]->passCutEq("first_max", adcs[0]-maxx, weight))) { continue; } bool helper = false; - if(doSample_==1){ + if (doSample_==1){ int len=*(&readout+1)-readout; - for(int KK=0;KKgetEventNumber(); - //std::cout<=2)){continue;} - if((regions_[i_reg]=="CTFit")and((thisHit->getT0(J)<26.0)or(thisHit->getT0(J)>30.0))){continue;} - sample(thisHit,regions_[i_reg],ievent,eventTime,N); + if ((regions_[i_reg]=="CTFit")&&((thisHit->getT0(j)<26.0)||(thisHit->getT0(j)>30.0))) { continue; } + sample(thisHit,regions_[i_reg], ievent, eventTime, N); } - - reg_histos_[regions_[i_reg]]->FillHistograms(thisHit,weight,J,i,TimeDiff,AmpDiff); - } + reg_histos_[regions_[i_reg]]->FillHistograms(thisHit, weight, j, I, TimeDiff, AmpDiff, stripID, hitc, hitl, otherTime); + } } } - //std::cout<getModule()); auto lay = std::to_string(thisHit->getLayer()); //swTag= mmapper_->getStringFromSw("ly"+lay+"_m"+mod); - std::string helper = mmapper_->getHwFromSw("ly"+lay+"_m"+mod); + std::string helper = mmapper_->getHwFromSw("ly" + lay + "_m" + mod); //std::cout<=2)and(word=="OneFit")){return;} - if((feb<2)and(word=="CTFit")){return;} - if(thisHit->getChiSq(0)<.85){return;} - - - - //std::cout<<"Feb "<getStrip()][0]<<" ,Baseline Feb>1: "<getStrip()][0]<<" More Baselines "<getStrip()][0]<<" ,Baseline Feb>1: "<getStrip()][0]<getStrip()); - }else{ - BigCount+=4096; - BigCount+=(feb-2)*2560+hyb*640+(int)(thisHit->getStrip()); + if ((feb>=2) && (word=="OneFit")) { return; } + if ((feb<2) && (word=="CTFit")) { return; } + if (thisHit->getChiSq(0) < .85) { return; } + + int count = 0; + if (feb <=1 ){ + count += feb*2048 + hyb*512 + (int)(thisHit->getStrip()); + } else { + count += 4096; + count += (feb-2)*2560 + hyb*640 + (int)(thisHit->getStrip()); } - //std::cout<<"READ HERE "<getStrip()<<" "<getStrip()][0]<getADCs(); - //std::cout<getADCs(); TF1* fitfunc = fourPoleFitFunction("Pulse 0",0); TF1* fitfunc2 = fourPoleFitFunction("Pulse 1",0); TF1* fitfunc3 = fourPoleFitFunction("Addition",1); - //TF1* baseline = new TF1("base","[0]",0.0,150.0); - //rETime((float(i))*24.0,T) - float TimeShift[2] = {-1*rETime(-(thisHit->getT0(0)),T),-1*rETime(-(thisHit->getT0(1)),T)}; - - fitfunc->FixParameter(0,TimeShift[0]); - fitfunc->FixParameter(3,thisHit->getAmp(0)); - if(feb<=1){ - fitfunc->FixParameter(1,times1_[feb][hyb][(int)thisHit->getStrip()][1]); - fitfunc->FixParameter(2,times1_[feb][hyb][(int)thisHit->getStrip()][2]); - fitfunc->FixParameter(4,baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); + float TimeShift[2] = {-1*rETime(-(thisHit->getT0(0)), T), -1*rETime(-(thisHit->getT0(1)), T)}; + + fitfunc->FixParameter(0, TimeShift[0]); + fitfunc->FixParameter(3, thisHit->getAmp(0)); + if (feb <=1 ){ + fitfunc->FixParameter(1, times1_[feb][hyb][(int)thisHit->getStrip()][1]); + fitfunc->FixParameter(2, times1_[feb][hyb][(int)thisHit->getStrip()][2]); + fitfunc->FixParameter(4, baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); //baseline->FixParameter(0,baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); - }else{ - fitfunc->FixParameter(1,times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - fitfunc->FixParameter(2,times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); - fitfunc->FixParameter(4,baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - //baseline->FixParameter(0,baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + } else { + fitfunc->FixParameter(1, times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + fitfunc->FixParameter(2, times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); + fitfunc->FixParameter(4, baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); } - if(thisHit->getFitN()==2){ - fitfunc2->FixParameter(0,TimeShift[1]); - fitfunc2->FixParameter(3,thisHit->getAmp(1)); - - fitfunc3->FixParameter(0,TimeShift[0]); - fitfunc3->FixParameter(3,thisHit->getAmp(0)); - fitfunc3->FixParameter(5,TimeShift[1]); - fitfunc3->FixParameter(8,thisHit->getAmp(1)); - - if(feb<=1){ - fitfunc2->FixParameter(1,times1_[feb][hyb][(int)thisHit->getStrip()][1]); - fitfunc2->FixParameter(2,times1_[feb][hyb][(int)thisHit->getStrip()][2]); - fitfunc2->FixParameter(4,baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); - - fitfunc3->FixParameter(1,times1_[feb][hyb][(int)thisHit->getStrip()][1]); - fitfunc3->FixParameter(2,times1_[feb][hyb][(int)thisHit->getStrip()][2]); - fitfunc3->FixParameter(4,baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); - - fitfunc3->FixParameter(6,times1_[feb][hyb][(int)thisHit->getStrip()][1]); - fitfunc3->FixParameter(7,times1_[feb][hyb][(int)thisHit->getStrip()][2]); - - }else{ - fitfunc2->FixParameter(1,times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - fitfunc2->FixParameter(2,times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); - fitfunc2->FixParameter(4,baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - - fitfunc3->FixParameter(1,times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - fitfunc3->FixParameter(2,times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); - fitfunc3->FixParameter(4,baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - - fitfunc3->FixParameter(6,times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); - fitfunc3->FixParameter(7,times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); + if (thisHit->getFitN()==2){ + fitfunc2->FixParameter(0, TimeShift[1]); + fitfunc2->FixParameter(3, thisHit->getAmp(1)); + + fitfunc3->FixParameter(0, TimeShift[0]); + fitfunc3->FixParameter(3, thisHit->getAmp(0)); + fitfunc3->FixParameter(5, TimeShift[1]); + fitfunc3->FixParameter(8, thisHit->getAmp(1)); + + if (feb <=1 ){ + fitfunc2->FixParameter(1, times1_[feb][hyb][(int)thisHit->getStrip()][1]); + fitfunc2->FixParameter(2, times1_[feb][hyb][(int)thisHit->getStrip()][2]); + fitfunc2->FixParameter(4, baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); + + fitfunc3->FixParameter(1, times1_[feb][hyb][(int)thisHit->getStrip()][1]); + fitfunc3->FixParameter(2, times1_[feb][hyb][(int)thisHit->getStrip()][2]); + fitfunc3->FixParameter(4, baseErr1_[feb][hyb][(int)thisHit->getStrip()][1]); + + fitfunc3->FixParameter(6, times1_[feb][hyb][(int)thisHit->getStrip()][1]); + fitfunc3->FixParameter(7, times1_[feb][hyb][(int)thisHit->getStrip()][2]); + + } else { + fitfunc2->FixParameter(1, times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + fitfunc2->FixParameter(2, times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); + fitfunc2->FixParameter(4, baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + + fitfunc3->FixParameter(1, times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + fitfunc3->FixParameter(2, times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); + fitfunc3->FixParameter(4, baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + + fitfunc3->FixParameter(6, times2_[feb-2][hyb][(int)thisHit->getStrip()][1]); + fitfunc3->FixParameter(7, times2_[feb-2][hyb][(int)thisHit->getStrip()][2]); } } - //for(int P=0;P<2;P++){for(int PP=0;PP<4;PP++){ std::cout<<"baseline: "<getStrip()][0]<getStrip()][0]<<" "<getStrip()][0]<vRange(0.0,3000.0,150.0,6000.0); - readout[K]++; - //auto gr = new TGraph(); - //auto gr2 = new TGraph(); - - std::string helper1="Feb: "+std::to_string(feb)+",Hyb: "+std::to_string(hyb)+",ch: "+std::to_string((int)thisHit->getStrip())+", chi_sqr value: "+std::to_string((float)thisHit->getChiSq(0)); + int Length = MatchList_.size(); + for (int k=0; kgetStrip()) + ", chi_sqr value: " + std::to_string((float)thisHit->getChiSq(0)); const char *thing1 = helper1.data(); - //gr2->SetPoint(0,0.,3000.);gr2->SetPoint(1,0.,6000.);gr2->SetPoint(2,150.0,6000.); - //gPad->DrawFrame(0.0,3000.0,150.0,6000.0); TCanvas *c1 = new TCanvas("c"); - c1->DrawFrame(0.0,3000.0,150.0,7000.0); + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); c1->SetTitle(thing1); - //auto gr = new TGraphErrors(); - //gr->SetName("ADCs"); float times[12]; float points[12]; float errors[12]; float zeroes[12]; - for(int i=0;i<6;i++){ - //std::cout<getStrip()][i]; - errors[i]=baseErr1_[feb][hyb][(int)thisHit->getStrip()][i+6]; - }else{ - //std::cout<<(float(i))*24.0-27.0<<" "<getStrip()][i]; - errors[i]=baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][i+6]; + for(int i=0; i<6; i++){ + zeroes[i] = 0.0; + if (feb <=1 ){ + times[i] = float(i)*24.0;//-27.0; + points[i] = adcs2[i]-baseErr1_[feb][hyb][(int)thisHit->getStrip()][i]; + errors[i] = baseErr1_[feb][hyb][(int)thisHit->getStrip()][i+6]; + } else { + times[i] = float(i)*24.0;//rETime((float(i))*24.0,T);//(float(i))*24.0-27.0; + points[i] = adcs2[i]-baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][i]; + errors[i] = baseErr2_[feb-2][hyb][(int)thisHit->getStrip()][i+6]; } } - auto gr = new TGraphErrors(6,times,points,zeroes,errors); + auto gr = new TGraphErrors(6, times, points, zeroes, errors); gr->SetName("ADCs"); gr->SetTitle(thing1); gr->GetYaxis()->SetTitle("ADC Counts"); gr->GetXaxis()->SetTitle("ns"); - gr->GetXaxis()->SetLimits(-10.0,130.); + gr->GetXaxis()->SetLimits(-10.0, 130.); gr->GetHistogram()->SetMaximum(2000.); gr->GetHistogram()->SetMinimum(-500.); - //gr2->Draw(""); gr->Draw("AL*"); - //TF1* helper1 = (TF1*)fitfunc->Clone("ff"); - - //baseline->SetLineColor(kBlue); - //baseline->Draw("same"); - + fitfunc->Draw("same"); - if(thisHit->getFitN()==2){ + if (thisHit->getFitN()==2){ fitfunc2->SetLineColor(kGreen); fitfunc2->Draw("same"); fitfunc3->SetLineColor(kOrange); fitfunc3->SetTitle(thing1); fitfunc3->Draw("same"); - - //TF1 *add = new TF1("ff+gg-base",); - //TF1* fitfunc3 = new TF1("ll","ff+gg+base"); - //TF1* h = new TF1("hh","ff+gg"); - //add->Draw("same"); - //auto ADD = new TF1(); } - //gr->Draw("same"); auto legend = new TLegend(0.1,0.7,.48,.9); - legend->AddEntry("gr","ADC counts"); - legend->AddEntry("base","Offline Baselines"); - legend->AddEntry("Pulse 0","First Pulse"); - if(thisHit->getFitN()==2){ - legend->AddEntry("Pulse 1","Second Pulse"); - legend->AddEntry("Addition","Summed Fit"); + legend->AddEntry("gr", "ADC counts"); + legend->AddEntry("base", "Offline Baselines"); + legend->AddEntry("Pulse 0", "First Pulse"); + if (thisHit->getFitN()==2){ + legend->AddEntry("Pulse 1", "Second Pulse"); + legend->AddEntry("Addition", "Summed Fit"); } legend->Draw("same"); - std::string helper2=word+std::to_string(readout[K]-1)+".png"; + std::string helper2 = word + std::to_string(readout[k]-1) + ".png"; const char *thing2 = helper2.data(); c1->SaveAs(thing2); - //gPad->Clear(); - - - //if(helper2=="OneFit8.png"){ - //std::cout<getStrip()<getStrip()][0]<<" "<getStrip()][1]<<" "<getStrip()][2]<<" "<getStrip()][3]<<" "<getStrip()][4]<<" "<getStrip()][5]<getAmp(0)<<" "<getT0(0)<<" "<getChiSq(0)<cd(); @@ -596,8 +492,6 @@ bool SvtRawDataAnaProcessor::process(IEvent* ievent) { std::string dirName = it->first; (it->second)->saveHistos(outF_,dirName); outF_->cd(dirName.c_str()); - //reg_selectors_[it->first]->getCutFlowHisto()->Scale(.5); - //reg_selectors_[it->first]->getCutFlowHisto()->Write(); } outF_->Close(); diff --git a/processors/src/TrackHitCompareAnaProcessor.cxx b/processors/src/TrackHitCompareAnaProcessor.cxx new file mode 100644 index 000000000..8dafe0b50 --- /dev/null +++ b/processors/src/TrackHitCompareAnaProcessor.cxx @@ -0,0 +1,775 @@ +/** + * @file TrackHitAnaProcessor.cxx + * @brief AnaProcessor used to compare two means of reconstruction directly. Fills histograms to study cluster reconstruction algorithms and dead channels. + * Does not feature a region selector or histomanager (i.e. no configurable json files), rather + * for the limited featured plots allows for more specific manipulation and control. + * @author Rory O'Dwyer and Cameron Bravo, SLAC National Accelerator Laboratory + */ +#include "TrackHitCompareAnaProcessor.h" +#include + +TrackHitCompareAnaProcessor::TrackHitCompareAnaProcessor(const std::string& name, Process& process) : Processor(name, process){ + mmapper_ = new ModuleMapper(2021); +} +TrackHitCompareAnaProcessor::~TrackHitCompareAnaProcessor(){} + +void TrackHitCompareAnaProcessor::configure(const ParameterSet& parameters) { + std::cout << "Configuring TrackHitAnaProcessor" << std::endl; + try + { + debug_ = parameters.getInteger("debug"); + layer_ = parameters.getInteger("layer"); + module_ = parameters.getInteger("module"); + isMC_ = parameters.getInteger("isMC"); + doingTracks_ = (parameters.getInteger("doTrack")==1); + pcut_ = (float)parameters.getDouble("cut"); + badchann_ = parameters.getString("badchannels"); + } + catch (std::runtime_error& error) + { + std::cout << error.what() << std::endl; + } + +} + + +void TrackHitCompareAnaProcessor::initialize(TTree* tree) { + fillDeads(); + tree_= tree; + if (isMC_ == 1){ + layers1_ = new TH1F("layers", "MC Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrk1_ = new TH1F("layersOnTrk", "MC Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrk1_ = new TH1F("layersOffTrk", "MC Strip Width for Clusters off Track", 12, 0.0, 12.0); + charges1_ = new TH1F("charges", "MC Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrk1_ = new TH1F("chargesOnTrk", "MC Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrk1_ = new TH1F("chargesOffTrk", "MC Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + layersNTD1_ = new TH1F("layersNTD", "MC Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrkNTD1_ = new TH1F("layersOnTrkNTD", "MC Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrkNTD1_ = new TH1F("layersOffTrkNTD", "MC Strip Width for Clusters off Track", 12, 0.0, 12.0); + chargesNTD1_ = new TH1F("chargesNTD", "MC Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrkNTD1_ = new TH1F("chargesOnTrkNTD", "MC Charge Distribution for On Track", 1000, 0.0, .000016); + chargesOffTrkNTD1_ = new TH1F("chargesOffTrkNTD", "MC Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + positions1_ = new TH1F("Positions", "MC Location of Cluster Hit;Layer;Hits", 14, 0.0, 14.0); + positionsOnTrk1_ = new TH1F("PositionsOnTrk", "MC Location of Cluster Hit for On Track", 14, 0.0, 14.0); + ClusDistances1_ = new TH1F("Minimum Cluster Difference", "MC Minimum Distance Between Clusters", 14, 0.0, 14.0); + ClusDistancesNTD1_ = new TH1F("Minimum Cluster Difference", "MC Minimum Distance Between Clusters", 14, 0.0, 14.0); + + times1_ = new TH1F("Times", "MC Time of Cluster Hit", 1000, -60.0, 60.0); + timesOnTrk1_ = new TH1F("TimesOnTrk", "MC Time of On Track Cluster Hit", 1000, -60.0, 60.0); + timesOffTrk1_ = new TH1F("TimesOffTrk", "MC Time of Off Cluster Hit", 1000, -60.0, 60.0); + timesNTD1_ = new TH1F("TimesNTD", "MC Time of Cluster Hit NTD", 1000, -60.0, 60.0); + timesOnTrkNTD1_ = new TH1F("TimesOnTrkNTD", "MC Time of On Track Cluster Hit NTD", 1000, -60.0, 60.0); + timesOffTrkNTD1_ = new TH1F("TimesOffTrkNTD", "MC Time of Off Cluster Hit NTD", 1000, -60.0, 60.0); + + tree_->SetBranchAddress("SiClusters", &Clusters_, &bClusters_); + tree_->SetBranchAddress("SiClustersOnTrack_KF", &ClustersKF_, &bClustersKF_); + tree_->SetBranchAddress("SVTRawTrackerHits", &svtraw_, &bsvtraw_); + return; + } + + // Instantiating the first layer + layers1_ = new TH1F("layers", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrk1_ = new TH1F("layersOnTrk", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrk1_ = new TH1F("layersOffTrk", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + charges1_ = new TH1F("charges", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrk1_ = new TH1F("chargesOnTrk", "Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrk1_ = new TH1F("chargesOffTrk", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + layersNTD1_ = new TH1F("layersNTD", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrkNTD1_ = new TH1F("layersOnTrkNTD", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrkNTD1_ = new TH1F("layersOffTrkNTD", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + chargesNTD1_ = new TH1F("chargesNTD", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrkNTD1_ = new TH1F("chargesOnTrkNTD", "Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrkNTD1_ = new TH1F("chargesOffTrkNTD", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + positions1_ = new TH1F("Positions", "Location of Cluster Hit;Layer;Hits", 14, 0.0, 14.0); + positionsOnTrk1_ = new TH1F("PositionsOnTrk", "Location of Cluster Hit for On Track", 14, 0.0, 14.0); + ClusDistances1_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + ClusDistancesNTD1_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + + times1_ = new TH1F("Times", "Time of Cluster Hit", 1000, -60.0, 60.0); + timesOnTrk1_ = new TH1F("TimesOnTrk", "Time of On Track Cluster Hit", 1000, -60.0, 60.0); + timesOffTrk1_ = new TH1F("TimesOffTrk", "Time of Off Cluster Hit", 1000, -60.0, 60.0); + timesNTD1_ = new TH1F("TimesNTD", "Time of Cluster Hit NTD", 1000, -60.0, 60.0); + timesOnTrkNTD1_ = new TH1F("TimesOnTrkNTD", "Time of On Track Cluster Hit NTD", 1000, -60.0, 60.0); + timesOffTrkNTD1_ = new TH1F("TimesOffTrkNTD", "Time of Off Cluster Hit NTD", 1000, -60.0, 60.0); + + if (doingTracks_){ + Z0VNShare2Hist1_ = new TH2F("Z0VNShare2Hist", "Z0 versus Number of Shared Hits No Cut", 100, 0, 3, 8, 0, 8); + Z0VNShare2HistCut1_ = new TH2F("Z0VNShare2HistCut", "Z0 versus Number of Shared Hits Momentum Cut", 100, 0, 3, 8, 0, 8); + SharedAmplitudes1_ = new TH1F("SharedAmplitudes", "The Amplitudes of Clusters Shared Between Tracks", 1000, 0.0, 0.000016); + UnSharedAmplitudes1_ = new TH1F("UnSharedAmplitudes", "The Amplitudes of Clusters Not Shared Between Tracks", 1000, 0.0, 0.000016); + SharedTimes1_ = new TH1F("SharedTimes", "The Times of Clusters Shared Between Tracks", 1000, -60.0, 60.0); + UnSharedTimes1_ = new TH1F("UnSharedTimes", "The Times of Clusters Not Shared Between Tracks", 1000, -60.0, 60.0); + TrackMomentumInTime1_ = new TH1F("TrackMomentumInTime", "The Momentum of In Time Tracks", 1000, 0.0, 7.0); + TrackMomentumOutTime1_ = new TH1F("TrackMomentumOutTime", "The Momentum of Out of Time Tracks", 1000, 0.0, 7.0); + TrackMomentumAllTime1_ = new TH1F("TrackMomentumAll", "The Momentum of All Tracks", 1000, 0.0, 7.0); + + TrackMomentumTInTime1_ = new TH1F("TrackMomentumTInTime", "The Transverse Momentum of In Time Tracks", 1000, 0.0, 7.0); + TrackMomentumTOutTime1_ = new TH1F("TrackMomentumTOutTime", "The Transverse Momentum of Out of Time Tracks", 1000, 0.0, 7.0); + TrackMomentumTAllTime1_ = new TH1F("TrackMomentumTAll", "The Transverse Momentum of All Tracks", 1000, 0.0, 7.0); + } + + + // Instantiating the second layer + layers2_ = new TH1F("layers", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrk2_ = new TH1F("layersOnTrk", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrk2_ = new TH1F("layersOffTrk", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + charges2_ = new TH1F("charges", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrk2_ = new TH1F("chargesOnTrk", "Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrk2_ = new TH1F("chargesOffTrk", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + layersNTD2_ = new TH1F("layersNTD", "Strip Width for All Clusters", 12, 0.0, 12.0); + layersOnTrkNTD2_ = new TH1F("layersOnTrkNTD", "Strip Width for Clusters on Track", 12, 0.0, 12.0); + layersOffTrkNTD2_ = new TH1F("layersOffTrkNTD", "Strip Width for Clusters off Track", 12, 0.0, 12.0); + chargesNTD2_ = new TH1F("chargesNTD", "Charge Distribution for All Clusters", 1000, 0.0, 0.000016); + chargesOnTrkNTD2_ = new TH1F("chargesOnTrkNTD", "Charge Distribution for On Track", 1000, 0.0, 0.000016); + chargesOffTrkNTD2_ = new TH1F("chargesOffTrkNTD", "Charge Distribution for Off Track", 1000, 0.0, 0.000016); + + times2_ = new TH1F("Times", "Time of Cluster Hit", 1000, -60.0, 60.0); + timesOnTrk2_ = new TH1F("TimesOnTrk", "Time of On Track Cluster Hit", 1000, -60.0, 60.0); + timesOffTrk2_ = new TH1F("TimesOffTrk", "Time of Off Cluster Hit", 1000, -60.0, 60.0); + timesNTD2_ = new TH1F("TimesNTD", "Time of Cluster Hit NTD", 1000, -60.0, 60.0); + timesOnTrkNTD2_ = new TH1F("TimesOnTrkNTD", "Time of On Track Cluster Hit NTD", 1000, -60.0, 60.0); + timesOffTrkNTD2_ = new TH1F("TimesOffTrkNTD", "Time of Off Cluster Hit NTD", 1000, -60.0, 60.0); + + positions2_ = new TH1F("Positions", "Location of Cluster Hit;Layer;Hits", 14, 0.0, 14.0); + positionsOnTrk2_ = new TH1F("PositionsOnTrk", "Location of Cluster Hit for On Track", 14, 0.0, 14.0); + ClusDistances2_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + ClusDistancesNTD2_ = new TH1F("Minimum Cluster Difference", "Minimum Distance Between Clusters", 14, 0.0, 14.0); + + if (doingTracks_){ + Z0VNShare2Hist2_ = new TH2F("Z0VNShare2Hist", "Z0 versus Number of Shared Hits No Cut", 100, 0, 3, 8, 0, 8); + Z0VNShare2HistCut2_ = new TH2F("Z0VNShare2HistCut", "Z0 versus Number of Shared Hits Momentum Cut", 100, 0, 3, 8, 0, 8); + SharedAmplitudes2_ = new TH1F("SharedAmplitudes", "The Amplitudes of Clusters Shared Between Tracks", 1000, 0.0, 0.000016); + UnSharedAmplitudes2_ = new TH1F("UnSharedAmplitudes", "The Amplitudes of Clusters Not Shared Between Tracks", 1000, 0.0, 0.000016); + SharedTimes2_ = new TH1F("SharedTimes", "The Times of Clusters Shared Between Tracks", 1000, -60.0, 60.0); + UnSharedTimes2_ = new TH1F("UnSharedTimes", "The Times of Clusters Not Shared Between Tracks", 1000, -60.0, 60.0); + + TrackMomentumInTime2_ = new TH1F("TrackMomentumInTime", "The Momentum of In Time Tracks", 1000, 0.0, 7.0); + TrackMomentumOutTime2_ = new TH1F("TrackMomentumOutTime", "The Momentum of Out of Time Tracks", 1000, 0.0, 7.0); + TrackMomentumAllTime2_ = new TH1F("TrackMomentumAll", "The Momentum of All Tracks", 1000, 0.0, 7.0); + + TrackMomentumTInTime2_ = new TH1F("TrackMomentumTInTime", "The Transverse Momentum of In Time Tracks", 1000, 0.0, 7.0); + TrackMomentumTOutTime2_ = new TH1F("TrackMomentumTOutTime", "The Transverse Momentum of Out of Time Tracks", 1000, 0.0, 7.0); + TrackMomentumTAllTime2_ = new TH1F("TrackMomentumTAll", "The Transverse Momentum of All Tracks", 1000, 0.0, 7.0); +} + tree_->SetBranchAddress("SiClusters", &Clusters_, &bClusters_); + tree_->SetBranchAddress("SiClustersOnTrack_KF", &ClustersKF_, &bClustersKF_); + tree_->SetBranchAddress("SVTRawTrackerHits", &svtraw_, &bsvtraw_); + tree_->SetBranchAddress("Identifier", &ident_, &bident_); + if (doingTracks_){ + tree_->SetBranchAddress("KalmanFullTracks", &tracks_, &btracks_); + } +} + +bool TrackHitCompareAnaProcessor::process(IEvent* ievent) { + if (doingTracks_){ + for (int i = 0; isize(); i++){ + Track* track = tracks_->at(i); + if (track->getTrackTime()*track->getTrackTime()<100.0){ + if (ident_ < 1.5){ + Z0VNShare2Hist1_->Fill(track->getZ0Err(), track->getNShared()); + } else { + Z0VNShare2Hist2_->Fill(track->getZ0Err(), track->getNShared()); + } + if (track->getP()Fill(track->getZ0Err(), track->getNShared()); + } else { + Z0VNShare2HistCut2_->Fill(track->getZ0Err(), track->getNShared()); + } + } + } + + // Track Momentum for In and Out of time hits (All not just transverse) + if (ident_ < 1.5) { + TrackMomentumAllTime1_->Fill(track->getP()); + } else { + TrackMomentumAllTime2_->Fill(track->getP()); + } + + if (track->getTrackTime()*track->getTrackTime() < 16.0){ + if (ident_ < 1.5){ + TrackMomentumInTime1_->Fill(track->getP()); + } else { + TrackMomentumInTime2_->Fill(track->getP()); + } + } else { + if (ident_ < 1.5) { + TrackMomentumOutTime1_->Fill(track->getP()); + } else { + TrackMomentumOutTime2_->Fill(track->getP()); + } + } + + // Transverse Track Momentum for same cuts + if (ident_ < 1.5){ + TrackMomentumTAllTime1_->Fill(track->getPt()); + } else { + TrackMomentumTAllTime2_->Fill(track->getPt()); + } + + if (track->getTrackTime()*track->getTrackTime() < 16.0){ + if (ident_ < 1.5){ + TrackMomentumTInTime1_->Fill(track->getPt()); + } else { + TrackMomentumTInTime2_->Fill(track->getPt()); + } + } else { + if (ident_ < 1.5){ + TrackMomentumTOutTime1_->Fill(track->getPt()); + } else { + TrackMomentumTOutTime2_->Fill(track->getPt()); + } + } + } + } + for(int i = 0; i < Clusters_->size(); i++){ + TrackerHit * clu = Clusters_->at(i); + Int_t layc = -1; + Int_t modc = -1; + + RawSvtHit * seed = (RawSvtHit*)(clu->getRawHits().At(0)); + + layc = clu->getLayer(); + modc = seed->getModule(); + + if (doingTracks_){ + bool isShared = false; + int increment = 0; + for (int i = 0; isize(); i++){ + Track* track = tracks_->at(i); + for (int j = 0; jgetSvtHits().GetEntries(); j++){ + TrackerHit * test = (TrackerHit *)(track->getSvtHits().At(j)); + if (clu->getTime()==test->getTime()){ + increment+=1; + } + } + } + if (increment > 1){ + isShared = true; + } + if (increment > 0){ + bool general = ((layer_==-1)||(module_==-1)); + if (((layc==layer_)&&(modc==module_))||(general)){ + if (isShared){ + if (ident_ < 1.5){ + SharedAmplitudes1_->Fill(clu->getCharge()); + } else { + SharedAmplitudes2_->Fill(clu->getCharge()); + } + if (ident_ < 1.5){ + SharedTimes1_->Fill(clu->getTime()); + } else { + SharedTimes2_->Fill(clu->getTime()); + } + } else { + if (ident_ < 1.5){ + UnSharedAmplitudes1_->Fill(clu->getCharge()); + } else { + UnSharedAmplitudes2_->Fill(clu->getCharge()); + } + if (ident_ < 1.5){ + UnSharedTimes1_->Fill(clu->getTime()); + } else { + UnSharedTimes2_->Fill(clu->getTime()); + } + } + } + } + } + + float seedStrip = (float)(seed->getStrip()); + float nLayers = (float)(clu->getRawHits().GetEntries()); + float ncharges = (float)(clu->getCharge()); + float ntimes = (float)(clu->getTime()); + bool onTrk = false; + bool NTD = false; + for (unsigned int j = 0; j < ClustersKF_->size(); j++){ + if (clu->getID()==(ClustersKF_->at(j)->getID())){ + onTrk = true; + } + } + std::string input = "ly" + std::to_string(layc+1) + "_m" + std::to_string(modc); + std::string helper = mmapper_->getHwFromSw(input); + int feb = std::stoi(helper.substr(1, 1)); + int hyb = std::stoi(helper.substr(3, 1)); + + int channelL = seedStrip-1; + int channelR = seedStrip+1; + if (channelL >= 0){ + NTD = (NTD)||(Deads_[GetStrip(feb, hyb, channelL)]==1); + } + if (((feb<=1)&&(channelR<=511))||((feb>1)&&(channelR<=639))){ + NTD = (NTD)||(Deads_[GetStrip(feb, hyb, channelR)]==1); + } + bool general = ((layer_==-1)||(module_==-1)); + if (((layc==layer_)&&(modc==module_))||(general)){ + // Now is the part where I fill the cluster distance histogram + float dist = 69420; + for (int p = 0; p < Clusters_->size(); p++){ + if (p==i) { continue; } + TrackerHit * clu2 = Clusters_->at(p); + RawSvtHit * seed2 = (RawSvtHit*)(clu2->getRawHits().At(0)); + float layc2 = clu->getLayer(); + float modc2 = seed->getModule(); + if ((!(layc2==layc))||(!(modc2==modc))) { continue; } + float new_dist = ((float)(seed2->getStrip())) - seedStrip; + if (new_dist < 0){ new_dist *= -1.0; } + if (new_dist < dist){ dist = new_dist; } + } + if (dist < 69420){ + if (ident_ < 1.5){ + ClusDistances1_->Fill(dist); + } else { + ClusDistances2_->Fill(dist); + } + if (NTD){ + if (ident_ < 1.5){ + ClusDistancesNTD1_->Fill(dist); + } else { + ClusDistancesNTD2_->Fill(dist); + } + } + } + if (ident_ < 1.5){ + layers1_->Fill(nLayers); + charges1_->Fill(ncharges); + positions1_->Fill(clu->getLayer()); + times1_->Fill(ntimes); + } else { + layers2_->Fill(nLayers); + charges2_->Fill(ncharges); + positions2_->Fill(clu->getLayer()); + times2_->Fill(ntimes); + } + if (onTrk){ + if (ident_ < 1.5){ + layersOnTrk1_->Fill(nLayers); + chargesOnTrk1_->Fill(ncharges); + timesOnTrk1_->Fill(ntimes); + positionsOnTrk1_->Fill(clu->getLayer()); + } else { + layersOnTrk2_->Fill(nLayers); + chargesOnTrk2_->Fill(ncharges); + timesOnTrk2_->Fill(ntimes); + positionsOnTrk2_->Fill(clu->getLayer()); + } + } else { + if (ident_ < 1.5){ + layersOffTrk1_->Fill(nLayers); + chargesOffTrk1_->Fill(ncharges); + timesOffTrk1_->Fill(ntimes); + } else { + layersOffTrk2_->Fill(nLayers); + chargesOffTrk2_->Fill(ncharges); + timesOffTrk2_->Fill(ntimes); + } + } + if (NTD){ + if (ident_ < 1.5){ + layersNTD1_->Fill(nLayers); + chargesNTD1_->Fill(ncharges); + timesNTD1_->Fill(ntimes); + } else { + layersNTD2_->Fill(nLayers); + chargesNTD2_->Fill(ncharges); + timesNTD2_->Fill(ntimes); + } + if (onTrk){ + if (ident_ < 1.5){ + layersOnTrkNTD1_->Fill(nLayers); + chargesOnTrkNTD1_->Fill(ncharges); + timesOnTrkNTD1_->Fill(ntimes); + } else { + layersOnTrkNTD2_->Fill(nLayers); + chargesOnTrkNTD2_->Fill(ncharges); + timesOnTrkNTD2_->Fill(ntimes); + } + } else { + if (ident_ < 1.5){ + layersOffTrkNTD1_->Fill(nLayers); + chargesOffTrkNTD1_->Fill(ncharges); + timesOffTrkNTD1_->Fill(ntimes); + } else { + layersOffTrkNTD2_->Fill(nLayers); + chargesOffTrkNTD2_->Fill(ncharges); + timesOffTrkNTD2_->Fill(ntimes); + } + } + } + } + } + return true; +} + +void TrackHitCompareAnaProcessor::fillDeads(){ + for (int i = 0; i<24576; i++){ + Deads_[i]=0.0; + } + std::string FILENAME = badchann_; + std::ifstream file(FILENAME.c_str()); + std::string line; + std::getline(file, line); + while (std::getline(file, line)) { + int value = std::atoi(line.c_str()); + Deads_[value] = 1.0; + } + file.close(); + return; +} + +int TrackHitCompareAnaProcessor::GetStrip(int feb, int hyb, int strip){ + int count = 0; + if (feb <= 1){ + count += feb*2048 + hyb*512 + strip; + } else { + count += 4096; + count += (feb-2)*2560 + hyb*640 + strip; + } + return count; +} + +void TrackHitCompareAnaProcessor::PlotClusterLayers(){ + TCanvas *c1 = new TCanvas("c"); + gPad->SetLogy(true); + c1->cd(); + layers1_->GetXaxis()->SetTitle("Width (strip)"); + layers1_->GetYaxis()->SetTitle("Hits"); + layersOnTrk1_->GetXaxis()->SetTitle("Width (strip)"); + layersOnTrk1_->GetYaxis()->SetTitle("Hits"); + layersOffTrk1_->GetXaxis()->SetTitle("Width (strip)"); + layersOffTrk1_->GetYaxis()->SetTitle("Hits"); + layers1_->SetLineColor(kRed); + layersOnTrk1_->SetLineColor(kRed); + layersOffTrk1_->SetLineColor(kRed); + + + layers2_->GetXaxis()->SetTitle("Width (strip)"); + layers2_->GetYaxis()->SetTitle("Hits"); + layersOnTrk2_->GetXaxis()->SetTitle("Width (strip)"); + layersOnTrk2_->GetYaxis()->SetTitle("Hits"); + layersOffTrk2_->GetXaxis()->SetTitle("Width (strip)"); + layersOffTrk2_->GetYaxis()->SetTitle("Hits"); + layers2_->SetLineColor(kBlue); + layersOnTrk2_->SetLineColor(kBlue); + layersOffTrk2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Layers for All Clusters"); + layers1_->Draw("e"); + layers2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(layers1_, "Layers Rec 1"); + legend1->AddEntry(layers2_, "Layers Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("allClusters.png"); + c1->Clear(); + + c1->SetTitle("Layers for On Track Clusters"); + layersOnTrk1_->Draw("e"); + layersOnTrk2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(layersOnTrk1_, "Layers On Track Rec 1"); + legend2->AddEntry(layersOnTrk2_, "Layers On Track Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("onClusters.png"); + c1->Clear(); + + c1->SetTitle("Layers for Off Track Clusters"); + layersOffTrk1_->Draw("e"); + layersOffTrk2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(layersOffTrk1_, "Layers Off Track Rec 1"); + legend3->AddEntry(layersOffTrk2_, "Layers Off Track Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("offClusters.png"); + c1->Clear(); + return; +} + +void TrackHitCompareAnaProcessor::PlotClusterCharges(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + + charges1_->GetXaxis()->SetTitle("Charge"); + charges1_->GetYaxis()->SetTitle("Hits"); + chargesOffTrk1_->GetXaxis()->SetTitle("Charge"); + chargesOffTrk1_->GetYaxis()->SetTitle("Hits"); + chargesOnTrk1_->GetXaxis()->SetTitle("Charge"); + chargesOnTrk1_->GetYaxis()->SetTitle("Hits"); + charges1_->SetLineColor(kRed); + chargesOnTrk1_->SetLineColor(kRed); + chargesOffTrk1_->SetLineColor(kRed); + + charges2_->GetXaxis()->SetTitle("Charge"); + charges2_->GetYaxis()->SetTitle("Hits"); + chargesOffTrk2_->GetXaxis()->SetTitle("Charge"); + chargesOffTrk2_->GetYaxis()->SetTitle("Hits"); + chargesOnTrk2_->GetXaxis()->SetTitle("Charge"); + chargesOnTrk2_->GetYaxis()->SetTitle("Hits"); + charges2_->SetLineColor(kBlue); + chargesOnTrk2_->SetLineColor(kBlue); + chargesOffTrk2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Charges for All Clusters"); + charges1_->Draw("e"); + charges2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(charges1_, "Charges Rec 1"); + legend1->AddEntry(charges2_, "Charges Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("allClustersCharge.png"); + c1->Clear(); + + c1->SetTitle("Charges for On Track Clusters"); + chargesOnTrk1_->Draw("e"); + chargesOnTrk2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(chargesOnTrk1_, "Charges On Track Rec 1"); + legend2->AddEntry(chargesOnTrk2_, "Charges On Track Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("onClustersCharge.png"); + c1->Clear(); + + c1->SetTitle("Charges for Off Track Clusters"); + chargesOffTrk1_->Draw("e"); + chargesOffTrk2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(chargesOffTrk1_, "Charges Off Track Rec 1"); + legend3->AddEntry(chargesOffTrk2_, "Charges Off Track Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("offClustersCharge.png"); + c1->Clear(); + return; +} + +void TrackHitCompareAnaProcessor::PlotClusterTimes(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + + times1_->GetXaxis()->SetTitle("Charge"); + times1_->GetYaxis()->SetTitle("Hits"); + timesOffTrk1_->GetXaxis()->SetTitle("Charge"); + timesOffTrk1_->GetYaxis()->SetTitle("Hits"); + timesOnTrk1_->GetXaxis()->SetTitle("Charge"); + timesOnTrk1_->GetYaxis()->SetTitle("Hits"); + times1_->SetLineColor(kRed); + timesOnTrk1_->SetLineColor(kRed); + timesOffTrk1_->SetLineColor(kRed); + + times2_->GetXaxis()->SetTitle("Charge"); + times2_->GetYaxis()->SetTitle("Hits"); + timesOffTrk2_->GetXaxis()->SetTitle("Charge"); + timesOffTrk2_->GetYaxis()->SetTitle("Hits"); + timesOnTrk2_->GetXaxis()->SetTitle("Charge"); + timesOnTrk2_->GetYaxis()->SetTitle("Hits"); + times2_->SetLineColor(kBlue); + timesOnTrk2_->SetLineColor(kBlue); + timesOffTrk2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Times for All Clusters"); + times1_->Draw("e"); + times2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(times1_, "Times Rec 1"); + legend1->AddEntry(times2_, "Times Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("alltimes.png"); + c1->Clear(); + + c1->SetTitle("Times for On Track Clusters"); + timesOnTrk1_->Draw("e"); + timesOnTrk2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(timesOnTrk1_, "Times On Track Rec 1"); + legend2->AddEntry(timesOnTrk2_, "Times On Track Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("ontimes.png"); + c1->Clear(); + + c1->SetTitle("Times for Off Track Clusters"); + timesOffTrk1_->Draw("e"); + timesOffTrk2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(timesOffTrk1_, "Times Off Track Rec 1"); + legend3->AddEntry(timesOffTrk2_, "Times Off Track Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("offtimes.png"); + c1->Clear(); + return; +} + +void TrackHitCompareAnaProcessor::PlotClusterPositions(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + + positions1_->GetXaxis()->SetTitle("Position"); + positions1_->GetYaxis()->SetTitle("Hits"); + positionsOnTrk1_->GetXaxis()->SetTitle("Position"); + positionsOnTrk1_->GetYaxis()->SetTitle("Hits"); + + positions2_->GetXaxis()->SetTitle("Position"); + positions2_->GetYaxis()->SetTitle("Hits"); + positionsOnTrk2_->GetXaxis()->SetTitle("Position"); + positionsOnTrk2_->GetYaxis()->SetTitle("Hits"); + + positions1_->SetLineColor(kRed); + positionsOnTrk1_->SetLineColor(kRed); + positions2_->SetLineColor(kBlue); + positionsOnTrk2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Position for All Clusters"); + positions1_->Draw("e"); + positions2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(positions1_, "Position Rec 1"); + legend1->AddEntry(positions2_, "Position Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("positions.png"); + c1->Clear(); + + c1->SetTitle("Position for On Track Clusters"); + positionsOnTrk1_->Draw("e"); + positionsOnTrk2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(positionsOnTrk1_, "Position On Track Rec 1"); + legend2->AddEntry(positionsOnTrk2_, "Position On Track Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("onpositions.png"); + c1->Clear(); + + c1->SetTitle("Cluster Distances"); + positionsOnTrk1_->Draw("e"); + positionsOnTrk2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(positionsOnTrk1_, "Cluster Distances Rec 1"); + legend3->AddEntry(positionsOnTrk2_, "Cluster Distances Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("clusDistances.png"); + c1->Clear(); + + return; +} + +void TrackHitCompareAnaProcessor::TrackMomenta(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + gPad->SetLogy(false); + TrackMomentumAllTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumAllTime1_->GetYaxis()->SetTitle("Momentum"); + TrackMomentumInTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumInTime1_->GetYaxis()->SetTitle("Momentum"); + TrackMomentumOutTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumOutTime1_->GetYaxis()->SetTitle("Momentum"); + + TrackMomentumAllTime1_->SetLineColor(kRed); + TrackMomentumInTime1_->SetLineColor(kRed); + TrackMomentumOutTime1_->SetLineColor(kRed); + TrackMomentumAllTime2_->SetLineColor(kBlue); + TrackMomentumInTime2_->SetLineColor(kBlue); + TrackMomentumOutTime2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Momentum for All Tracks"); + TrackMomentumAllTime1_->Draw("e"); + TrackMomentumAllTime2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(TrackMomentumAllTime1_, "Momenta Rec 1"); + legend1->AddEntry(TrackMomentumAllTime2_, "Momenta Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("allTracksMomenta.png"); + c1->Clear(); + + c1->SetTitle("Momentum for In Time Tracks"); + TrackMomentumInTime1_->Draw("e"); + TrackMomentumInTime2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(TrackMomentumInTime1_, "Momenta Rec 1"); + legend2->AddEntry(TrackMomentumInTime2_, "Momenta Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("inTimeTracksMomenta.png"); + c1->Clear(); + + gPad->SetLogy(true); + c1->SetTitle("Momentum for Out Time Tracks"); + TrackMomentumOutTime1_->Draw("e"); + TrackMomentumOutTime2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(TrackMomentumOutTime1_, "Momenta Rec 1"); + legend3->AddEntry(TrackMomentumOutTime2_, "Momenta Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("outTimeTracksMomenta.png"); + c1->Clear(); + return; +} + +void TrackHitCompareAnaProcessor::TrackTransverseMomenta(){ + TCanvas *c1 = new TCanvas("c"); + c1->cd(); + gPad->SetLogy(false); + TrackMomentumTAllTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumTAllTime1_->GetYaxis()->SetTitle("Momentum"); + TrackMomentumTInTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumTInTime1_->GetYaxis()->SetTitle("Momentum"); + TrackMomentumTOutTime1_->GetXaxis()->SetTitle("Track Number"); + TrackMomentumTOutTime1_->GetYaxis()->SetTitle("Momentum"); + + TrackMomentumTAllTime1_->SetLineColor(kRed); + TrackMomentumTInTime1_->SetLineColor(kRed); + TrackMomentumTOutTime1_->SetLineColor(kRed); + TrackMomentumTAllTime2_->SetLineColor(kBlue); + TrackMomentumTInTime2_->SetLineColor(kBlue); + TrackMomentumTOutTime2_->SetLineColor(kBlue); + + c1->DrawFrame(0.0, 3000.0, 150.0, 7000.0); + c1->SetTitle("Transverse Momentum for All Tracks"); + TrackMomentumTAllTime1_->Draw("e"); + TrackMomentumTAllTime2_->Draw("e same"); + auto legend1 = new TLegend(0.3, 0.8, .68, .9); + legend1->AddEntry(TrackMomentumTAllTime1_, "Momenta Rec 1"); + legend1->AddEntry(TrackMomentumTAllTime2_, "Momenta Rec 2"); + legend1->Draw("same e"); + c1->SaveAs("allTracksMomentaT.png"); + c1->Clear(); + + c1->SetTitle("Transverse Momentum for In Time Tracks"); + TrackMomentumTInTime1_->Draw("e"); + TrackMomentumTInTime2_->Draw("e same"); + auto legend2 = new TLegend(0.3, 0.8, .68, .9); + legend2->AddEntry(TrackMomentumTInTime1_, "Momenta Rec 1"); + legend2->AddEntry(TrackMomentumTInTime2_, "Momenta Rec 2"); + legend2->Draw("same e"); + c1->SaveAs("inTimeTracksMomentaT.png"); + c1->Clear(); + gPad->SetLogy(true); + c1->SetTitle("Transverse Momentum for Out Time Tracks"); + TrackMomentumTOutTime1_->Draw("e"); + TrackMomentumTOutTime2_->Draw("e same"); + auto legend3 = new TLegend(0.3, 0.8, .68, .9); + legend3->AddEntry(TrackMomentumTOutTime1_, "Momenta Rec 1"); + legend3->AddEntry(TrackMomentumTOutTime2_, "Momenta Rec 2"); + legend3->Draw("same e"); + c1->SaveAs("outTimeTracksMomentaT.png"); + c1->Clear(); + return; +} + +void TrackHitCompareAnaProcessor::finalize() { + PlotClusterLayers(); + PlotClusterCharges(); + PlotClusterTimes(); + PlotClusterPositions(); + if(doingTracks_){ + TrackMomenta(); + TrackTransverseMomenta(); + } + return; +} +DECLARE_PROCESSOR(TrackHitCompareAnaProcessor); diff --git a/scripts/AddIdentity.C b/scripts/AddIdentity.C new file mode 100644 index 000000000..abee9f9ca --- /dev/null +++ b/scripts/AddIdentity.C @@ -0,0 +1,18 @@ +/** + * Simple ROOT macro + * Associates to each event in a ROOT file an identifier which can be used post hadd + * to process a single file with multiple reconstruction algorithms used in a single processor. + **/ +void AddIdentity(){ + TFile *f = new TFile("rootcopy0.root", "update"); + TTree *tree = (TTree *)(f->Get("HPS_Event")); + Float_t helper = 1.0; + auto my_new_branch = tree->Branch("Identifier", &helper, "Identifier/F"); + for (Long64_t entry = 0; entry < tree->GetEntries(); entry++) { + //tree->GetEntry(); + /* something to compute my_local_variable */ + std::cout << entry << std::endl; + my_new_branch->Fill(); + } + tree->Write("", TObject::kOverwrite); +}