diff --git a/musrSimAna/40041.v1190 b/musrSimAna/40041.v1190 new file mode 100644 index 0000000..e34dccd --- /dev/null +++ b/musrSimAna/40041.v1190 @@ -0,0 +1,128 @@ +# TDC V1190 Set up file generated by Deltat : 16-Apr-2008 10:30:27 +# bins are always in TDC channel units (195.3125 ps/channel) + +INSTRUMENT=GPS + +DESCRIPTION=No Veto -- 1port +TYPE=TDCV1190 +RESOLUTION=100 + +MDELAY=0 +PDELAY=2000 +#REWINDTIMEBINS=10000000 + +MCOINCIDENCEW=50 +PCOINCIDENCEW=50 +VCOINCIDENCEW=100 + +MUONRATEFACTOR=0.0991035 +#MUONRATEFACTOR=1000 +#INFINITELYLOWMUONRATE + +DATAWINDOWMIN=-5 +DATAWINDOWMAX=10.0 +PILEUPWINDOWMIN=-10.2 +PILEUPWINDOWMAX=10.2 + +102; "M up"; M; 0.4; 2005; ; + 1; "B1"; P; 0.4; 2005; ; B1; 1485; 1515; 50995; + 2; "B2"; P; 0.4; 2005; ; B2; 1485; 1515; 50995; + 11; "F1"; P; 0.4; 2005; ; F11; 1485; 1515; 50995; + 12; "F1"; P; 0.4; 2005; ; F12; 1485; 1515; 50995; + 13; "F1"; P; 0.4; 2005; ; F13; 1485; 1515; 50995; +! 3; "F right"; K; 0.4; 1263; +! 4; "B left"; P; 0.2; 1311; 5 -10 -11 -12 -13; Back; 1485; 1515; 50995; +! 5; "B right"; K; 0.2; 1262; +! 6; Up; P; 0.2; 1306; -10 -11 -12 -13; Up; 1490; 1520; 51000; +! 7; Down; P; 0.2; 1264; -10 -11 -12 -13; Down; 1525; 1555; 51035; +! 8; "R int"; P; 0.2; 1309; -10 -11 -12 -13; R_int; 1480; 1510; 50990; +! 9; "R ext"; P; 0.3; 1267; -10 -11 -12 -13; R_ext; 1520; 1550; 51020; 8; Rite +!10; "Bveto up"; V; 0.3; 2085; +!11; "Bveto down"; V; 0.3; 2089; +!12; "Bveto left"; V; 0.3; 2080; +!13; "Bveto right"; V; 0.3; 2079; +!14; "F center"; P; 0.2; 2073; -10 -11 -12 -13; F_cntr; 720; 750; 50230; 2; Forw +!15; Clock; C; +!16; Rejected; R; +!17; Ch17; N; +!18; Ch18; N; +!19; Ch19; N; +!20; Ch20; N; +!21; Ch21; N; +!22; Ch22; N; +!23; Ch23; N; +!24; Ch24; N; +!25; Ch25; N; +!26; Ch26; N; +!27; Ch27; N; +!28; Ch28; N; +!29; Ch29; N; +!30; Ch30; N; +!31; Ch31; N; + +counterPhaseShifts 1 0 2 0 11 180 12 180 13 180 + +musrTH2D hMuDecayMap "Muon decay map;z [mm]; r[mm]" 100 -80. 20. 40 0. 20. muDecayPosZ muDecayPosR +musrTH1D hMuDecayPosZ "Penetration of muons into the sample;z[mm];N" 100 -5.0 5. muDecayPosZ +musrTH1D hMuDecayPosX "X of decayed muons;x[mm];N" 100 -5.0 5. muDecayPosX +musrTH1D hdet_m0edep "Energy deposited in the trigger;E[MeV];N" 100 0. 2 det_m0edep +musrTH1D hdet_posEdep "Energy deposited in the positron Counters;E[MeV];N" 100 0. 2 det_posEdep +#musrTH1D hdet_time0 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time0 +#musrTH1D hdet_time1 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time1 +musrTH1D hdet_time10 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time10 +musrTH1D hgen_time10 "generated muSR spectra;time[#mus];N" 220 -11. 11. gen_time10 +musrTH1D hdet_time10_MINUS_gen_time10 "detected - generated ;time_{det}-time_{gen} [ps];N" 100 -500. 500. det_time10_MINUS_gen_time10 +#musrTH1D hMuDecayDetID "Detector ID where the muons decay;Detector ID; N" 4000 -2000. 2000. muDecayDetID +musrTH1D hdet_time20 "detected muSR spectra (phase shifted);time[#mus];N" 700 -2.0 12.0 det_time20 +musrTH1D hMuDecayDetID "Detector ID where the muons decay;Detector ID; N" 2002 -1001. 1001. muDecayDetID +musrTH1D hMuDecayDetIDpileup "Detector ID where the pileup muons decay;Detector ID; N" 2002 -1001. 1001. pileup_muDecayDetID +humanDecayHistograms hMuDecayDetID hMuDecayDetIDpileup 201 Sample 211 SampleCell 300 Holder 301 Holder 302 Holder 303 Holder 304 Holder 255 PbCollim 251 CuCollim 261 Collim1 102 M0 1 BC 2 BC 11 FC 12 FC 13 FC -1 World 231 World -1000 Escaped + +condition 0 alwaysTrue +condition 1 oncePerEvent +condition 2 muonDecayedInSample_gen +condition 3 muonTriggered_gen +condition 4 muonTriggered_det +condition 5 positronHit_det +condition 6 goodEvent_det +condition 7 goodEvent_gen +condition 8 goodEvent_det_AND_goodEvent_gen +condition 9 pileupEvent + +#fit hdet_time20 pol0 -1.9 -0.1 0 +#fit hdet_time20 funct1 0.2 9.8 0 0 0.3 100 +#fit hdet_time20 funct2 0.2 9.8 0 0 0.3 100 0 +fit hdet_time20 funct4 0.2 9.8 0 0 0.3 100 0 + +#draw hMuDecayPosZ 0 +draw hMuDecayPosZ 1 +#draw hMuDecayPosZ 2 +#draw hMuDecayPosZ 3 +#draw hMuDecayPosZ 4 +###draw hMuDecayMap 0 +###draw hMuDecayMap 1 +draw hdet_m0edep 0 +draw hdet_m0edep 1 +draw hdet_m0edep 2 +draw hdet_m0edep 3 +draw hdet_m0edep 4 +draw hdet_m0edep 5 +draw hdet_m0edep 6 +draw hdet_time10 6 +draw hgen_time10 7 +draw hdet_time10_MINUS_gen_time10 8 +draw hdet_posEdep 5 +draw hMuDecayDetID 1 +draw hdet_m0edep 9 +draw hdet_time10 9 +draw hMuDecayDetIDpileup 1 +draw hMuDecayDetIDpileup 9 +draw humanDecayHistograms 1 +draw humanDecayHistograms 6 +draw humanDecayHistograms 9 +draw humanDecayPileupHistograms 1 +draw humanDecayPileupHistograms 6 +draw humanDecayPileupHistograms 9 +draw hdet_time20 4 +draw hdet_time20 6 +$ diff --git a/musrSimAna/40042.v1190 b/musrSimAna/40042.v1190 new file mode 100644 index 0000000..db64a8b --- /dev/null +++ b/musrSimAna/40042.v1190 @@ -0,0 +1,128 @@ +# TDC V1190 Set up file generated by Deltat : 16-Apr-2008 10:30:27 +# bins are always in TDC channel units (195.3125 ps/channel) + +INSTRUMENT=GPS + +DESCRIPTION=No Veto -- 1port +TYPE=TDCV1190 +RESOLUTION=100 + +MDELAY=0 +PDELAY=2000 +#REWINDTIMEBINS=10000000 + +MCOINCIDENCEW=50 +PCOINCIDENCEW=50 +VCOINCIDENCEW=100 + +MUONRATEFACTOR=0.102712 +#MUONRATEFACTOR=1000 +#INFINITELYLOWMUONRATE + +DATAWINDOWMIN=-5 +DATAWINDOWMAX=10.0 +PILEUPWINDOWMIN=-10.2 +PILEUPWINDOWMAX=10.2 + +102; "M up"; M; 0.4; 2005; ; + 1; "B1"; P; 0.4; 2005; ; B1; 1485; 1515; 50995; + 2; "B2"; P; 0.4; 2005; ; B2; 1485; 1515; 50995; + 11; "F1"; P; 0.4; 2005; ; F11; 1485; 1515; 50995; + 12; "F1"; P; 0.4; 2005; ; F12; 1485; 1515; 50995; + 13; "F1"; P; 0.4; 2005; ; F13; 1485; 1515; 50995; +! 3; "F right"; K; 0.4; 1263; +! 4; "B left"; P; 0.2; 1311; 5 -10 -11 -12 -13; Back; 1485; 1515; 50995; +! 5; "B right"; K; 0.2; 1262; +! 6; Up; P; 0.2; 1306; -10 -11 -12 -13; Up; 1490; 1520; 51000; +! 7; Down; P; 0.2; 1264; -10 -11 -12 -13; Down; 1525; 1555; 51035; +! 8; "R int"; P; 0.2; 1309; -10 -11 -12 -13; R_int; 1480; 1510; 50990; +! 9; "R ext"; P; 0.3; 1267; -10 -11 -12 -13; R_ext; 1520; 1550; 51020; 8; Rite +!10; "Bveto up"; V; 0.3; 2085; +!11; "Bveto down"; V; 0.3; 2089; +!12; "Bveto left"; V; 0.3; 2080; +!13; "Bveto right"; V; 0.3; 2079; +!14; "F center"; P; 0.2; 2073; -10 -11 -12 -13; F_cntr; 720; 750; 50230; 2; Forw +!15; Clock; C; +!16; Rejected; R; +!17; Ch17; N; +!18; Ch18; N; +!19; Ch19; N; +!20; Ch20; N; +!21; Ch21; N; +!22; Ch22; N; +!23; Ch23; N; +!24; Ch24; N; +!25; Ch25; N; +!26; Ch26; N; +!27; Ch27; N; +!28; Ch28; N; +!29; Ch29; N; +!30; Ch30; N; +!31; Ch31; N; + +counterPhaseShifts 1 0 2 0 11 180 12 180 13 180 + +musrTH2D hMuDecayMap "Muon decay map;z [mm]; r[mm]" 100 -80. 20. 40 0. 20. muDecayPosZ muDecayPosR +musrTH1D hMuDecayPosZ "Penetration of muons into the sample;z[mm];N" 100 -5.0 5. muDecayPosZ +musrTH1D hMuDecayPosX "X of decayed muons;x[mm];N" 100 -5.0 5. muDecayPosX +musrTH1D hdet_m0edep "Energy deposited in the trigger;E[MeV];N" 100 0. 2 det_m0edep +musrTH1D hdet_posEdep "Energy deposited in the positron Counters;E[MeV];N" 100 0. 2 det_posEdep +#musrTH1D hdet_time0 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time0 +#musrTH1D hdet_time1 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time1 +musrTH1D hdet_time10 "detected muSR spectra;time[#mus];N" 220 -11. 11. det_time10 +musrTH1D hgen_time10 "generated muSR spectra;time[#mus];N" 220 -11. 11. gen_time10 +musrTH1D hdet_time10_MINUS_gen_time10 "detected - generated ;time_{det}-time_{gen} [ps];N" 100 -500. 500. det_time10_MINUS_gen_time10 +#musrTH1D hMuDecayDetID "Detector ID where the muons decay;Detector ID; N" 4000 -2000. 2000. muDecayDetID +musrTH1D hdet_time20 "detected muSR spectra (phase shifted);time[#mus];N" 700 -2.0 12.0 det_time20 +musrTH1D hMuDecayDetID "Detector ID where the muons decay;Detector ID; N" 2002 -1001. 1001. muDecayDetID +musrTH1D hMuDecayDetIDpileup "Detector ID where the pileup muons decay;Detector ID; N" 2002 -1001. 1001. pileup_muDecayDetID +humanDecayHistograms hMuDecayDetID hMuDecayDetIDpileup 201 Sample 211 SampleCell 300 Holder 301 Holder 302 Holder 303 Holder 304 Holder 255 PbCollim 251 CuCollim 261 Collim1 102 M0 1 BC 2 BC 11 FC 12 FC 13 FC -1 World 231 World -1000 Escaped + +condition 0 alwaysTrue +condition 1 oncePerEvent +condition 2 muonDecayedInSample_gen +condition 3 muonTriggered_gen +condition 4 muonTriggered_det +condition 5 positronHit_det +condition 6 goodEvent_det +condition 7 goodEvent_gen +condition 8 goodEvent_det_AND_goodEvent_gen +condition 9 pileupEvent + +#fit hdet_time20 pol0 -1.9 -0.1 0 +#fit hdet_time20 funct1 0.2 9.8 0 0 0.3 100 +#fit hdet_time20 funct2 0.2 9.8 0 0 0.3 100 0 +fit hdet_time20 funct4 1.5 9.8 0 0 0.3 100 0 + +#draw hMuDecayPosZ 0 +draw hMuDecayPosZ 1 +#draw hMuDecayPosZ 2 +#draw hMuDecayPosZ 3 +#draw hMuDecayPosZ 4 +###draw hMuDecayMap 0 +###draw hMuDecayMap 1 +draw hdet_m0edep 0 +draw hdet_m0edep 1 +draw hdet_m0edep 2 +draw hdet_m0edep 3 +draw hdet_m0edep 4 +draw hdet_m0edep 5 +draw hdet_m0edep 6 +draw hdet_time10 6 +draw hgen_time10 7 +draw hdet_time10_MINUS_gen_time10 8 +draw hdet_posEdep 5 +draw hMuDecayDetID 1 +draw hdet_m0edep 9 +draw hdet_time10 9 +draw hMuDecayDetIDpileup 1 +draw hMuDecayDetIDpileup 9 +draw humanDecayHistograms 1 +draw humanDecayHistograms 6 +draw humanDecayHistograms 9 +draw humanDecayPileupHistograms 1 +draw humanDecayPileupHistograms 6 +draw humanDecayPileupHistograms 9 +draw hdet_time20 4 +draw hdet_time20 6 +$ diff --git a/musrSimAna/Makefile b/musrSimAna/Makefile new file mode 100644 index 0000000..040294a --- /dev/null +++ b/musrSimAna/Makefile @@ -0,0 +1,20 @@ +#include /home/daquser/temp/SctRodDaq/build/Makefile.include +#include /home/daquser/temp/SctRodDaq/build/Makefile.commands +# +# +CPP=g++ +#CPPFLAGS= -I$(ROOTSYS)/include -I/userdisk/sedlak/myLCG/CondDBMySQL_new/CondDBMySQL/include -I/usr/include/mysql +#LDFLAGS= $(shell root-config --glibs) -L/userdisk/sedlak/myLCG/CondDBMySQL_new/CondDBMySQL/src/.libs -lconddb -L/usr/lib/mysql -lmysqlclient +CPPFLAGS= -I$(ROOTSYS)/include +LDFLAGS= $(shell root-config --glibs) -lMinuit + +all: musrRoot + +%.o:%.cxx + $(CPP) -c $(CPPFLAGS) $^ + +musrRoot: musrSimAna.o musrAnalysis.o musrCounter.o musrTH.o + $(CPP) -o musrSimAna $(LDFLAGS) $^ + +clean: + -rm -f *.o musrSimAna diff --git a/musrSimAna/musrAnalysis.cxx b/musrSimAna/musrAnalysis.cxx new file mode 100644 index 0000000..76dae6a --- /dev/null +++ b/musrSimAna/musrAnalysis.cxx @@ -0,0 +1,1030 @@ +#define musrAnalysis_cxx +#include "musrAnalysis.hh" +#include "musrTH.hh" +//#include +#include +#include +#include +#include +//#include + +void musrAnalysis::Loop(char* runChar, char* v1190FileName, Int_t nrEvents) +{ +// In a ROOT session, you can do: +// Root > .L musrAnalysis.C +// Root > musrAnalysis t +// Root > t.GetEntry(12); // Fill t data members with entry number 12 +// Root > t.Show(); // Show values of entry 12 +// Root > t.Show(16); // Read and show values of entry 16 +// Root > t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + hInfo = new TH1D("hInfo","Different values that need to be passed to ploting program",10000,1.,10001.); + ReadInInputParameters(v1190FileName); + CreateHistograms(); + + nentries = fChain->GetEntriesFast(); + if (nrEvents!=0) nentries = nrEvents; + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); InitialiseEvent(); nbytes += nb; + // if (Cut(ientry) < 0) continue; + if (jentry%100000==0) std::cout << "Analysing event nr. jentry=" << jentry << std::endl; + AnalyseEvent(jentry); + } + SaveHistograms(runChar,v1190FileName); +} + + + +//================================================================ + +void musrAnalysis::ReadInInputParameters(char* charV1190FileName) { + std::cout<<"musrAnalysis::ReadInInputParameters()"< tmpCoincidenceMultimap; + // for (int i=0; i("a", 1)); + // + // } + // } + while (!feof(fSteeringFile)) { + fgets(line,1000,fSteeringFile); + if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')&&(line[0]!='$')&&(line[0]!='!')) { + char tmpString0[200]="Unset"; + char tmpString1[200]="Unset"; + sscanf(&line[0],"%s",tmpString0); + if (strncmp(tmpString0,"INSTRUMENT=",strlen("INSTRUMENT="))==0) { + instrument=&line[strlen("INSTRUMENT=")]; + } + else if (strncmp(tmpString0,"DESCRIPTION=",strlen("DESCRIPTION="))==0) { + description=&line[strlen("DESCRIPTION=")]; + } + else if (strncmp(tmpString0,"TYPE=",strlen("TYPE="))==0) { + tdchwtype=&line[strlen("TYPE=")]; + } + else if (strncmp(tmpString0,"RESOLUTION=",strlen("RESOLUTION="))==0) { + sscanf(&line[strlen("RESOLUTION=")],"%g",&tmpfvar); + tdcresolution = tmpfvar * picosecond; + if (tdcresolution<1*picosecond) { + std::cout<<"\n !!! ERROR !!! Time resolution of TDC must be larger than 1 ps (you requested "< too low time resolution might lead to numerical problems in musrCounter ==> S T O P"< S T O P !!!"< ignored, histogram not defined! S T O P !!!"< ignored, histogram not defined! S T O P !!!"<second; + if (strcmp(tmpString0,"musrTH1D")==0) { + // std::cout<<"pointerToVariable1="<SetRotatingReferenceFrame(rot_ref_frequency, rot_ref_phase); + hInfo->Fill(21, rot_ref_frequency); // value may be overwritten - just the last rot. ref. frame will be saved + hInfo->Fill(22, rot_ref_phase); // value may be overwritten - just the last rot. ref. frame will be saved + } + } + else { + Double_t* pointerToVariable2 = iter2->second; + musrTH* myTH = new musrTH("2D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,pointerToVariable2); + listOfAllHistograms2D.push_back(myTH); + } + } + else if (strncmp(tmpString0,"humanDecayHistograms",strlen("humanDecayHistograms"))==0) { + char motherHistoName[200]; + char motherHistoPileupName[200]; + sscanf(&line[0],"%*s %s %s",motherHistoName,motherHistoPileupName); + std::cout<<" Defining humanDecayHistograms (motherHistoName="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(motherHistoName)) { + motherOfHumanDecayHistograms = *it; + } + if ((*it)->IsThisHistoNamed(motherHistoPileupName)) { + motherOfHumanDecayPileupHistograms = *it; + } + } + if (motherOfHumanDecayHistograms==NULL) { + std::cout<<"\n\nmusrAnalysis::ReadInInputParameters Mother histogram (motherOfHumanDecayHistograms)" + <<"of the \"Human decay histogram\" not found!"< S T O P "< S T O P "<second)==(std::string) NAME1) { + nameDefined = true; + iBinHuman = it->first; + } + } + if (!nameDefined) { + nBinHuman++; + iBinHuman=nBinHuman; + humanDecayMap.insert(std::pair(iBinHuman,NAME1)); + } + humanDecayMultimap.insert(std::pair(iBinHuman,N1)); + char* pch2 = strstr(pch,NAME1)+strlen(NAME1); + pch=pch2; + nscan = sscanf(pch,"%d %s",&N1,NAME1); + // indexHuman++; + } while (nscan==2); + + Int_t nrBinsHumanHist=humanDecayMap.size(); + // Int_t nrBinsHumanHist=0; + // std::string oldString="blaBlaoldStringUndefiNED"; + // for (humanDecayMultimapType::const_iterator it = humanDecayMultimap.begin(); it!=humanDecayMultimap.end(); ++it) { + // if (oldString!=(it->first)) nrBinsHumanHist++; + // } + std::cout<<"\nnrBinsHumanHist="<SetBinLabel1D(k+1,sHumanDecayArray[k]); + // } + // oldString="blaBlaoldStringUndefiNED"; + for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) { + humanDecayHistograms -> SetBinLabel1D(it->first,it->second); + humanDecayPileupHistograms -> SetBinLabel1D(it->first,it->second); + } + } + else if (strncmp(tmpString0,"condition",strlen("condition"))==0) { + // Selection of the predefined conditions to be applied + int iConditionTMP; + char conditionNameTMP[200]; + sscanf(&line[0],"%*s %d %s",&iConditionTMP,conditionNameTMP); + if (iConditionTMP<0) { + std::cout<<" !!! ERROR: Condition number must be >= 0 (not "<=nrConditions) { + std::cout<<" !!! ERROR: Condition number must be < "< Add it in the musrAnalysis.cxx S T O P !!!"<= -1 (not "<=nrConditions) { + std::cout<<" !!! ERROR: draw: condition number must be < "<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoNameTMP)) { + (*it)->SetDrawListOfHistograms(iConditionTMP); + } + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoNameTMP)) { + (*it)->SetDrawListOfHistograms(iConditionTMP); + } + } + if (humanDecayHistograms->IsThisHistoNamed(histoNameTMP)) { + humanDecayHistograms->SetDrawListOfHistograms(iConditionTMP); + } + if (humanDecayPileupHistograms->IsThisHistoNamed(histoNameTMP)) { + humanDecayPileupHistograms->SetDrawListOfHistograms(iConditionTMP); + } + } + else if (strcmp(tmpString0,"counterPhaseShifts")==0) { + int nscan; int N1; float PHASE_SHIFT; char NAME1[100]; char NAME2[100]; + char *pch = line + strlen("counterPhaseShifts"); + do { + nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT); + phaseShiftMap.insert(std::pair(N1,PHASE_SHIFT/180*pi)); + std::cout<<"N1="< SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + } + else if (strcmp(functionName,"funct2")==0) { + funct = new TF1("funct2","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + funct -> SetParameter(5,p4); + } + else if (strcmp(functionName,"funct3")==0) { + funct = new TF1("funct3","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))+[5]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> FixParameter(4,xMin); + funct -> FixParameter(5,p4); + } + else if (strcmp(functionName,"funct4")==0) { + funct = new TF1("funct4","[3]*exp(-x/2.19703)*(1+[2]*cos(x*[0]+[1]))+[4]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + funct -> SetParameter(3,p3); + funct -> SetParameter(4,p4); + } + else if (strcmp(functionName,"pol0")==0) { + funct = new TF1("pol0","pol0"); + } + else if (strcmp(functionName,"simpleExpoPLUSconst")==0) { + funct = new TF1("simpleExpoPLUSconst","[0]*exp(-x/2.19703)+[1]"); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + } + else if (strcmp(functionName,"rotFrameTime20")==0) { + funct = new TF1("rotFrameTime20","[2]*exp(-x/2.19703)*cos(x*[0]+[1]) "); + funct -> SetParameter(0,p0); + funct -> SetParameter(1,p1); + funct -> SetParameter(2,p2); + } + else { + std::cout<<"musrAnalysis::ReadInInputParameters: function \""< S T O P"<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + if ((*it)->IsThisHistoNamed(histoName)) { + tmpHistograms = *it; + } + } + if (tmpHistograms==NULL) { + std::cout<<"\n\nmusrAnalysis::ReadInInputParameters: Fit required, but the corresponding histogram " + <<"not found!"< S T O P "<AssignFunction(funct, xMin, xMax); + } + + else { + //------------------------------------------------------------------ + // Second part of the configuration file + // + // Replace spaces between quotation marks by '_' + // and remove quatation marks and semicolon by ' '. + // Remember position of the beginning and end of the defintion of the coincidence counters + int nrOfSemicolons=0, beginningOfCoincidenceArray=0, endOfCoincidenceArray =0; + for (Int_t i=0; i<(int)strlen(line); i++) { + if (line[i]=='"') { + line[i] = ' '; + Int_t ii=i; + do { + ii++; + if (isspace(line[ii])) line[ii] = '_'; + } while (line[ii]!='"'); + line[ii]= ' '; + } + else if (line[i] == ';') { + line[i] = ' '; + nrOfSemicolons++; + if (nrOfSemicolons==5) beginningOfCoincidenceArray=i; + if (nrOfSemicolons==6) endOfCoincidenceArray=i; + } + } + + // Read in the values from the configuration file + int tmpChannel=-1, tmpTimeShift=0; + char tmpName[200], tmpType[200]; + float tmpThreshold =0.; + sscanf(&line[0],"%d %s %s %g %d",&tmpChannel,tmpName,tmpType,&tmpThreshold,&tmpTimeShift); + musrCounter* counter = new musrCounter(tmpChannel,tmpName,tmpType[0],tmpThreshold,tmpTimeShift); + if (tmpType[0]=='M') {mCounter =counter; allCounterMap[tmpChannel]=counter; overallBinDelay=tmpTimeShift;} + else if (tmpType[0]=='P') {pCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} + else if (tmpType[0]=='K') {kCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} + else if (tmpType[0]=='V') {vCounterMap[tmpChannel]=counter; allCounterMap[tmpChannel]=counter;} + + // Read in the values for the coincidence detectors + // (Tricky while the number of detectors in coincidence is not known) + if ((endOfCoincidenceArray>0)&&(endOfCoincidenceArray!=beginningOfCoincidenceArray)) { + char tmpCoinc[100]; strcpy(tmpCoinc,"987654321"); + int jj=beginningOfCoincidenceArray; + do { + do {jj++;} while (isspace(line[jj])); // skip space characters; + if (jj>endOfCoincidenceArray) {break;} + sscanf(&line[jj],"%s",tmpCoinc); + jj+=strlen(tmpCoinc); + int iCoinc=atoi(tmpCoinc); + if (iCoinc==987654321) continue; +//ckdel 20.9.2010 if ((maxChannels<=tmpChannel)||(maxChannels<=iCoinc)) { +//ckdel 20.9.2010 std::cout<<" Channel number ("< S T O P"<(tmpChannel,iCoinc)); + // tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc; + } while (jjSetTDChistogram(tmpHistoName,tmpT0,tmpT1,tmpT2,tmpHistoaddNr,tmpHistoNameAdd); + } + } + } + // std::cout << line; + } + } + std::cout << "INSTRUMENT=" << instrument << std::endl; + std::cout << "DESCRIPTION=" << description << std::endl; + std::cout << "TYPE=" << tdchwtype << std::endl; + std::cout << "RESOLUTION="<second)->GetCounterType(); + if (DetectorType=='M') { + (it->second)->SetMaxCoincidenceTimeWindow(pileupWindowBinMin); + (it->second)->SetCoincidenceTimeWindow(pileupWindowBinMin,pileupWindowBinMax); + (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-mcoincwin,-mcoincwin,mcoincwin); + } + else if (DetectorType=='P') { + (it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin); + (it->second)->SetCoincidenceTimeWindow(dataWindowBinMin,dataWindowBinMax); + (it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors(-pcoincwin,-pcoincwin,pcoincwin); + } + else if (DetectorType=='V') { + (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin); + (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin); + } + int iChanNr = (it->second)->GetCounterNr(); + for (std::multimap::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { + if ((*itCoinc).first == iChanNr) { + int iChn2 = itCoinc->second; + counterMapType::const_iterator itTwo = allCounterMap.find(abs(iChn2)); + if (itTwo==allCounterMap.end()) { + std::cout<<" Pointer to coincidence counter ("< S T O P"<second)->SetCoincidenceCounter(itTwo->second,iChn2); + } + } + } + } + + // for (int j=0; j::iterator itCoinc = tmpCoincidenceMultimap.begin(); itCoinc != tmpCoincidenceMultimap.end(); ++itCoinc) { + // std::cout << " houby [" << (*itCoinc).first << ", " << (*itCoinc).second << "]" << std::endl; + + // + // if (tmpCoincidenceMultimap[iChanNr][j]!=0) { +// int j = (*itCoinc).first; +// counterMapType::const_iterator itTwo = allCounterMap.find(j); +// if (itTwo==allCounterMap.end()) { +// std::cout<<" Pointer to coincidence counter ("< S T O P"<second; +// (it->second)->SetCoincidenceCounter(itTwo->second,iChn); +// } +// // } +// } +// } + + + + // + // + // char tmpString1[100]="Unset",tmpString2[100]="Unset",tmpString3[100]="Unset"; + // sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2); + +} + +//================================================================ + +void musrAnalysis::CreateHistograms() { + std::cout<<"musrAnalysis::CreateHistograms()"<pileupWindowMax)? 3*dataWindowMax: 3*pileupWindowMax; + std::cout<<"musrAnalysis::CreateHistograms(): safeTimeWindow = "<GetEntry(1); Double_t fieldValue = fieldNomVal[0]; if (nFieldNomVal>1) fieldValue += fieldNomVal[1]; + // omega = 851.372*fabs(fieldValue); + // omega = 851.372*fieldValue; // value set originally in this program + // omega = 850.62*fieldValue; // value used in Geant ? + omega = 851.610*fieldValue; // value from the fits of the data + + + hInfo->Fill(1, fieldValue); + hInfo->Fill(6, runID); + hInfo->Fill(7, hGeantParameters->GetBinContent(7)); + std::cout<<"musrAnalysis::CreateHistograms(): fieldValue="<DrawTH1D("hist",0); + // (*it)->DrawTH1D("hist",1); + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + (*it)->DrawTH2DdrawList("hist"); + } + for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) { + (it->second)->DrawTDChistogram(); + } + if (humanDecayHistograms) { + // humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,indexHuman,nHumanDecayArray); + humanDecayHistograms->FillHumanDecayArray(motherOfHumanDecayHistograms,humanDecayMap,humanDecayMultimap); + humanDecayHistograms->DrawTH1DdrawList("text"); + humanDecayPileupHistograms->FillHumanDecayArray(motherOfHumanDecayPileupHistograms,humanDecayMap,humanDecayMultimap); + humanDecayPileupHistograms->DrawTH1DdrawList("text"); + } + + Long64_t numberOfMuonCandidates = mCounter->GetNumberOfMuonCandidates(); + hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers" + + //============================== + // Write out the histograms + std::cout<<"musrAnalysis::SaveHistograms()"<GetList()); + while ( TObject *obj = next() ) { + if (obj->InheritsFrom(TH1::Class()) ) {obj->Write();} + } + fHistograms->Close(); + delete fHistograms; + + + //============================== + + std::cout<<"==================================================================="<lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)GetEntry(iiiEntry); InitialiseEvent(); + + // Loop over all interesting "moments", which are: + // 1) any hit in the muon counter + // 2) any new event (even if nothing was detected in the muon counter) + + Long64_t currentEventID = eventID; + Double_t timeOfM0Hit; + Long64_t evID; + Bool_t thisEventAlreadyFilled = false; + + // std::cout<<"___________________\n"; + // for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) { + // (*it).second->myPrintThisCounter(eventID); + // } + FillHistograms(iiiEntry); + + + currentTime+=timeToNextEvent*muonRateFactor; + Long64_t currentTimeBin = Long64_t(currentTime/tdcresolution); + if (!boolInfinitelyLowMuonRate) RemoveOldHitsFromCounters(currentTimeBin-1); // Remove all hits that can not influence the next event + else RemoveOldHitsFromCounters(std::numeric_limits::max()/2); + if (currentTimeBin>rewindTimeBins) { // When the time variables start to be too large, rewind the time info everywhere; + RewindAllTimeInfo(rewindTimeBins); + } +} + +//================================================================ + +//void musrAnalysis::RemoveOldHitsFromCounters(Double_t timeLimit) { +void musrAnalysis::RemoveOldHitsFromCounters(Long64_t timeBinLimit) { + // Loop over all Counters and remove hits that happen at or before the "timeBinLimit". + // std::cout<<"musrAnalysis::RemoveOldHitsFromCounters: timeBinLimit="<RemoveHitsInCounter(timeBinLimit); + } +} + +//================================================================ + +//void musrAnalysis::RewindAllTimeInfo(Double_t timeToRewind) { +void musrAnalysis::RewindAllTimeInfo(Long64_t timeBinsToRewind) { + currentTime -= timeBinsToRewind*tdcresolution; + nextUnfilledEventTime -= timeBinsToRewind*tdcresolution; + numberOfRewinds++; + // Loop over all timing information and rewind it by "timeToRewind"; + // std::cout<<"musrAnalysis::RewindAllTimeInfo()"<RewindHitsInCounter(timeBinsToRewind); + } +} + +//================================================================ + +void musrAnalysis::FillHistograms(Int_t iiiEntry) { + // std::cout<<"musrAnalysis::FillHistograms() event="<::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) { + (*it)->FillTH1D(wght,&condition[0]); + } + for(std::list::const_iterator it = listOfAllHistograms2D.begin(); it != listOfAllHistograms2D.end(); ++it) { + (*it)->FillTH2D(wght,&condition[0]); + } + if ((goodEvent_det)&&(pCounterHitExistsForThisEventID)) { + // musrCounter* pCounterHitInThisEvent + counterMapType::const_iterator itPcounterHitInThisEvent = pCounterMap.find(idetP_ID); + if (itPcounterHitInThisEvent==pCounterMap.end()) { + std::cout<<" ERROR! musrAnalysis: itPcounterHitInThisEvent not found ==> This should never happen!!!"<GetBranch("weight") ) ? weight : 1.; +} + +//================================================================ +Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) { + // std::cout<<"musrAnalysis::PreprocessEvent()"<GetEntry(iEn); InitialiseEvent(); + // std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<::iterator it; + it = allCounterMap.find(det_ID[i]); + if (it==allCounterMap.end()) { + // uncomment std::cout<<"Active detector with det_ID="<FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]); + } + } + + // std::cout<<"lastPreprocessedEntry+1="<second)->GetNextGoodHitInThisEvent(evID,dataBinMin,dataBinMax,'P',tBin1,kEntry,idetP,idetP_ID,idetP_edep,doubleHit); + Bool_t thereWasHit = (it->second)->GetNextGoodPositron(evID,dataBinMin,dataBinMax,tBin1,tBin2,kEntry,idetP,idetP_ID,idetP_edep,doubleHit); + // std::cout<<"000000000000 tBin1="<first; + } + std::cout<second); + } + std::cout<-1000.00001) ) return -1000.; + if ( (x<-999.99999) && (x>-1000.00001) ) return -1000.; + if ( (x==0) && (y==0) ) return -1000.; + return atan2(y,x); +} diff --git a/musrSimAna/musrAnalysis.hh b/musrSimAna/musrAnalysis.hh new file mode 100644 index 0000000..9fcf924 --- /dev/null +++ b/musrSimAna/musrAnalysis.hh @@ -0,0 +1,567 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Thu Mar 25 15:00:51 2010 by ROOT version 5.24/00 +// from TTree t1/a simple Tree with simple variables +// found on file: data/musr_40003.root +////////////////////////////////////////////////////////// + +#ifndef musrAnalysis_h +#define musrAnalysis_h + +#include +#include +#include +#include +#include +#include +#include +#include +#include "musrCounter.hh" +//#include "musrTH.hh" +#include + +class musrTH; + +class musrAnalysis { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + Int_t fCurrent; //!current Tree number in a TChain + + // Declaration of leaf types + Int_t runID; + Int_t eventID; + Double_t weight; + Double_t timeToNextEvent; + Double_t BFieldAtDecay_Bx; + Double_t BFieldAtDecay_By; + Double_t BFieldAtDecay_Bz; + Double_t BFieldAtDecay_B3; + Double_t BFieldAtDecay_B4; + Double_t BFieldAtDecay_B5; + Double_t muIniTime; + Double_t muIniPosX; + Double_t muIniPosY; + Double_t muIniPosZ; + Double_t muIniMomX; + Double_t muIniMomY; + Double_t muIniMomZ; + Double_t muIniPolX; + Double_t muIniPolY; + Double_t muIniPolZ; + Int_t muDecayDetID; + Double_t muDecayPosX; + Double_t muDecayPosY; + Double_t muDecayPosZ; + Double_t muDecayTime; + Double_t muDecayPolX; + Double_t muDecayPolY; + Double_t muDecayPolZ; + Double_t muTargetTime; + Double_t muTargetPolX; + Double_t muTargetPolY; + Double_t muTargetPolZ; + Double_t muTargetMomX; + Double_t muTargetMomY; + Double_t muTargetMomZ; + Double_t muM0Time; + Double_t muM0PolX; + Double_t muM0PolY; + Double_t muM0PolZ; + Double_t posIniMomX; + Double_t posIniMomY; + Double_t posIniMomZ; + Int_t nFieldNomVal; + Double_t fieldNomVal[50]; //[nFieldNomVal] + Int_t det_n; + Int_t det_ID[500]; //[det_n] + Double_t det_edep[500]; //[det_n] + Double_t det_edep_el[500]; //[det_n] + Double_t det_edep_pos[500]; //[det_n] + Double_t det_edep_gam[500]; //[det_n] + Double_t det_edep_mup[500]; //[det_n] + Int_t det_nsteps[500]; //[det_n] + Double_t det_length[500]; //[det_n] + Double_t det_time_start[500]; //[det_n] + Double_t det_time_end[500]; //[det_n] + Double_t det_x[500]; //[det_n] + Double_t det_y[500]; //[det_n] + Double_t det_z[500]; //[det_n] + Double_t det_kine[500]; //[det_n] + Double_t det_VrtxKine[500]; //[det_n] + Double_t det_VrtxX[500]; //[det_n] + Double_t det_VrtxY[500]; //[det_n] + Double_t det_VrtxZ[500]; //[det_n] + Int_t det_VrtxVolID[500]; //[det_n] + Int_t det_VrtxProcID[500]; //[det_n] + Int_t det_VrtxTrackID[500]; //[det_n] + Int_t det_VrtxParticleID[500]; //[det_n] + Double_t det_VvvKine[500]; //[det_n] + Double_t det_VvvX[500]; //[det_n] + Double_t det_VvvY[500]; //[det_n] + Double_t det_VvvZ[500]; //[det_n] + Int_t det_VvvVolID[500]; //[det_n] + Int_t det_VvvProcID[500]; //[det_n] + Int_t det_VvvTrackID[500]; //[det_n] + Int_t det_VvvParticleID[500]; //[det_n] + + // List of branches + TBranch *b_runID; //! + TBranch *b_eventID; //! + TBranch *b_weight; //! + TBranch *b_timeToNextEvent; //! + TBranch *b_BFieldAtDecay; //! + TBranch *b_muIniTime; //! + TBranch *b_muIniPosX; //! + TBranch *b_muIniPosY; //! + TBranch *b_muIniPosZ; //! + TBranch *b_muIniMomX; //! + TBranch *b_muIniMomY; //! + TBranch *b_muIniMomZ; //! + TBranch *b_muIniPolX; //! + TBranch *b_muIniPolY; //! + TBranch *b_muIniPolZ; //! + TBranch *b_muDecayDetID; //! + TBranch *b_muDecayPosX; //! + TBranch *b_muDecayPosY; //! + TBranch *b_muDecayPosZ; //! + TBranch *b_muDecayTime; //! + TBranch *b_muDecayPolX; //! + TBranch *b_muDecayPolY; //! + TBranch *b_muDecayPolZ; //! + TBranch *b_muTargetTime; //! + TBranch *b_muTargetPolX; //! + TBranch *b_muTargetPolY; //! + TBranch *b_muTargetPolZ; //! + TBranch *b_muTargetMomX; //! + TBranch *b_muTargetMomY; //! + TBranch *b_muTargetMomZ; //! + TBranch *b_muM0Time; //! + TBranch *b_muM0PolX; //! + TBranch *b_muM0PolY; //! + TBranch *b_muM0PolZ; //! + TBranch *b_posIniMomX; //! + TBranch *b_posIniMomY; //! + TBranch *b_posIniMomZ; //! + TBranch *b_nFieldNomVal; //! + TBranch *b_fieldNomVal; //! + TBranch *b_det_n; //! + TBranch *b_det_ID; //! + TBranch *b_det_edep; //! + TBranch *b_det_edep_el; //! + TBranch *b_det_edep_pos; //! + TBranch *b_det_edep_gam; //! + TBranch *b_det_edep_mup; //! + TBranch *b_det_nsteps; //! + TBranch *b_det_length; //! + TBranch *b_det_time_start; //! + TBranch *b_det_time_end; //! + TBranch *b_det_x; //! + TBranch *b_det_y; //! + TBranch *b_det_z; //! + TBranch *b_det_kine; //! + TBranch *b_det_VrtxKine; //! + TBranch *b_det_VrtxX; //! + TBranch *b_det_VrtxY; //! + TBranch *b_det_VrtxZ; //! + TBranch *b_det_VrtxVolID; //! + TBranch *b_det_VrtxProcID; //! + TBranch *b_det_VrtxTrackID; //! + TBranch *b_det_VrtxParticleID; //! + TBranch *b_det_VvvKine; //! + TBranch *b_det_VvvX; //! + TBranch *b_det_VvvY; //! + TBranch *b_det_VvvZ; //! + TBranch *b_det_VvvVolID; //! + TBranch *b_det_VvvProcID; //! + TBranch *b_det_VvvTrackID; //! + TBranch *b_det_VvvParticleID; //! + + typedef std::map variableMapType; + variableMapType variableMap; + Double_t runID_double; + Double_t eventID_double; + Double_t muDecayDetID_double; + Double_t det_n_double; + Double_t muDecayPosR; + Double_t wght; + Double_t det_m0edep; + Double_t det_posEdep; + Double_t pos_Trans_Momentum; + Double_t pos_Momentum; + Double_t pos_Radius; + Double_t pos_Theta; + Double_t pos_Phi; + // Double_t det_time0; + // Double_t get_time0; + // Double_t det_time1; + // Double_t gen_time1; + Double_t det_time10; + Double_t gen_time10; + Double_t det_time10_MINUS_gen_time10; + Double_t det_time20; + Double_t pileup_muDecayDetID_double; + + + typedef std::map conditionMapType; + conditionMapType conditionMap; + Bool_t alwaysTrue; + Bool_t oncePerEvent; + Bool_t muonDecayedInSample_gen; + Bool_t muonTriggered_gen; + Bool_t muonTriggered_det; + Bool_t positronHit_det; + Bool_t goodEvent_det; + Bool_t goodEvent_gen; + Bool_t goodEvent_det_AND_goodEvent_gen; + Bool_t pileupEventCandidate; + Bool_t pileupEvent; + Bool_t goodEvent_det_AND_muonDecayedInSample_gen; + + musrAnalysis(TTree *tree=0); + virtual ~musrAnalysis(); + virtual Int_t Cut(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry); + virtual Long64_t LoadTree(Long64_t entry); + virtual void Init(TTree *tree); + virtual void Loop(char* runChar, char* v1190FileName, Int_t nrEvents); + virtual Bool_t Notify(); + virtual void Show(Long64_t entry = -1); + virtual void ReadInInputParameters(char* charV1190FileName); + virtual void CreateHistograms(); + virtual void AnalyseEvent(Long64_t iiiEntry); + virtual void FillHistograms(Int_t iiiEntry); + virtual void SaveHistograms(char* runChar,char* v1190FileName); + virtual void RemoveOldHitsFromCounters(Long64_t timeBinLimit); + // virtual void RewindAllTimeInfo(Double_t timeToRewind); + virtual void RewindAllTimeInfo(Long64_t timeBinsToRewind); + virtual void InitialiseEvent(); + virtual Double_t PreprocessEvent(Long64_t iEn); + virtual Bool_t PositronCounterHit(Int_t evID, Long64_t dataBinMin, Long64_t dataBinMax, Long64_t& tBin1, Long64_t& tBin2, Int_t& kEntry, Int_t& idetP, Int_t& idetP_ID, Double_t& idetP_edep, Bool_t& doubleHit); + void CopySubstring(char* inputChar,int iStart,int iEnd,char* outputChar); + void MyPrintTree(); + void MyPrintConditions(); + Double_t myAtan2(Double_t y, Double_t x); + + typedef std::map phaseShiftMapType; + phaseShiftMapType phaseShiftMap; + TH1D* hGeantParameters; + TH1D* hInfo; + + static const Double_t pi=3.14159265358979324; + + private: + std::string instrument; + std::string description; + std::string tdchwtype; + Double_t tdcresolution; + Long64_t nentries; + Int_t mdelay; + Int_t pdelay; + Int_t mcoincwin; + Int_t pcoincwin; + Int_t vcoincwin; + Double_t muonRateFactor; + Long64_t rewindTimeBins; + Double_t dataWindowMin; + Double_t dataWindowMax; + Double_t pileupWindowMin; + Double_t pileupWindowMax; + Long64_t pileupWindowBinMin; + Long64_t pileupWindowBinMax; + Long64_t dataWindowBinMin; + Long64_t dataWindowBinMax; + Int_t overallBinDelay; + Bool_t boolInfinitelyLowMuonRate; + + char musrMode; // D = time diferential; I = time integral; + Double_t safeTimeWindow; + Double_t currentTime; + // Double_t nextEventTime; + Double_t nextUnfilledEventTime; + Long64_t numberOfRewinds; + Long64_t numberOfGoodMuons; + // Int_t currentEventID; + Long64_t lastPreprocessedEntry; + musrCounter* mCounter; + typedef std::map counterMapType; + counterMapType pCounterMap; + counterMapType kCounterMap; + counterMapType vCounterMap; + counterMapType allCounterMap; + Int_t testIVar1; + Double_t omega; + + static const Double_t microsecond=1.; + static const Double_t nanosecond=0.001; + static const Double_t picosecond=0.000001; + // static const Double_t rewindTime=1000000.; + // static const Double_t rewindTime=1000.; + // static const Long64_t rewindTimeBins=1000000000000000; // Max Long64_t can be +9,223,372,036,854,775,807 + // static const Long64_t rewindTimeBins = 1000000000; + // static const Double_t dataWindowMin = -0.2; + // static const Double_t dataWindowMax = 10.0; + // static const Double_t pileupWindowMin = -10.5; + // static const Double_t pileupWindowMax = 10.5; +//ckdel 20.9.2010 static const int maxChannels=32; + +public: + static const Int_t nrConditions = 31; + Bool_t condition[nrConditions]; + +private: + // HISTOGRAMS: + std::list listOfAllHistograms1D; + std::list listOfAllHistograms2D; + TH1D** hTargetZ; + TH2D** hMuDecayMap; + // + musrTH* humanDecayHistograms; + musrTH* motherOfHumanDecayHistograms; + musrTH* humanDecayPileupHistograms; + musrTH* motherOfHumanDecayPileupHistograms; + // int indexHuman; + // int nHumanDecayArray[50]; + // std::string sHumanDecayArray[50]; + + typedef std::map humanDecayMapType; + typedef std::multimap humanDecayMultimapType; + humanDecayMapType humanDecayMap; + humanDecayMultimapType humanDecayMultimap; +}; + +#endif + +#ifdef musrAnalysis_cxx +musrAnalysis::musrAnalysis(TTree *tree) +{ + variableMap["muDecayPosX"]=&muDecayPosX; + variableMap["muDecayPosY"]=&muDecayPosY; + variableMap["muDecayPosZ"]=&muDecayPosZ; + variableMap["runID"]=&runID_double; + variableMap["eventID"]=&eventID_double; + variableMap["weight"]=&weight; + variableMap["timeToNextEvent"]=&timeToNextEvent; + variableMap["BFieldAtDecay_Bx"]=&BFieldAtDecay_Bx; + variableMap["BFieldAtDecay_By"]=&BFieldAtDecay_By; + variableMap["BFieldAtDecay_Bz"]=&BFieldAtDecay_Bz; + variableMap["BFieldAtDecay_B3"]=&BFieldAtDecay_B3; + variableMap["BFieldAtDecay_B4"]=&BFieldAtDecay_B4; + variableMap["BFieldAtDecay_B5"]=&BFieldAtDecay_B5; + variableMap["muIniTime"]=&muIniTime; + variableMap["muIniPosX"]=&muIniPosX; + variableMap["muIniPosY"]=&muIniPosY; + variableMap["muIniPosZ"]=&muIniPosZ; + variableMap["muIniMomX"]=&muIniMomX; + variableMap["muIniMomY"]=&muIniMomY; + variableMap["muIniMomZ"]=&muIniMomZ; + variableMap["muIniPolX"]=&muIniPolX; + variableMap["muIniPolY"]=&muIniPolY; + variableMap["muIniPolZ"]=&muIniPolZ; + variableMap["muDecayDetID"]=&muDecayDetID_double; + variableMap["muDecayPosX"]=&muDecayPosX; + variableMap["muDecayPosY"]=&muDecayPosY; + variableMap["muDecayPosZ"]=&muDecayPosZ; + variableMap["muDecayTime"]=&muDecayTime; + variableMap["muDecayPolX"]=&muDecayPolX; + variableMap["muDecayPolY"]=&muDecayPolY; + variableMap["muDecayPolZ"]=&muDecayPolZ; + variableMap["muTargetTime"]=&muTargetTime; + variableMap["muTargetPolX"]=&muTargetPolX; + variableMap["muTargetPolY"]=&muTargetPolY; + variableMap["muTargetPolZ"]=&muTargetPolZ; + variableMap["muTargetMomX"]=&muTargetMomX; + variableMap["muTargetMomY"]=&muTargetMomY; + variableMap["muTargetMomZ"]=&muTargetMomZ; + variableMap["muM0Time"]=&muM0Time; + variableMap["muM0PolX"]=&muM0PolX; + variableMap["muM0PolY"]=&muM0PolY; + variableMap["muM0PolZ"]=&muM0PolZ; + variableMap["posIniMomX"]=&posIniMomX; + variableMap["posIniMomY"]=&posIniMomY; + variableMap["posIniMomZ"]=&posIniMomZ; + // variableMap["nFieldNomVal"]=&nFieldNomVal_double; + // variableMap["fieldNomVal0"]=...; //[nFieldNomVal] + variableMap["det_n"]=&det_n_double; + // + variableMap["muDecayPosR"]=&muDecayPosR; + variableMap["wght"]=&wght; + variableMap["det_m0edep"]=&det_m0edep; + variableMap["det_posEdep"]=&det_posEdep; + variableMap["pos_Trans_Momentum"]=&pos_Trans_Momentum; + variableMap["pos_Momentum"]=&pos_Momentum; + variableMap["pos_Radius"]=&pos_Radius; + variableMap["pos_Theta"]=&pos_Theta; + variableMap["pos_Phi"]=&pos_Phi; + // variableMap["det_time0"]=&det_time0; + // variableMap["gen_time0"]=&gen_time0; + // variableMap["det_time1"]=&det_time1; + // variableMap["gen_time1"]=&gen_time1; + variableMap["det_time10"]=&det_time10; + variableMap["gen_time10"]=&gen_time10; + variableMap["det_time10_MINUS_gen_time10"]=&det_time10_MINUS_gen_time10; + variableMap["pileup_muDecayDetID"]=&pileup_muDecayDetID_double; + variableMap["det_time20"]=&det_time20; + + testIVar1=0; + humanDecayHistograms=NULL; + motherOfHumanDecayHistograms=NULL; + humanDecayPileupHistograms=NULL; + motherOfHumanDecayPileupHistograms=NULL; + +// if parameter tree is not specified (or zero), connect the file +// used to generate this class and read the Tree. + if (tree == 0) { + TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("data/musr_40003.root"); + if (!f) { + f = new TFile("data/musr_40003.root"); + } + tree = (TTree*)gDirectory->Get("t1"); + + } + hGeantParameters = (TH1D*) gDirectory->Get("hGeantParameters"); + Init(tree); +} + +musrAnalysis::~musrAnalysis() +{ + if (!fChain) return; + delete fChain->GetCurrentFile(); +} + +Int_t musrAnalysis::GetEntry(Long64_t entry) +{ +// Read contents of entry. + if (!fChain) return 0; + return fChain->GetEntry(entry); +} + +Long64_t musrAnalysis::LoadTree(Long64_t entry) +{ +// Set the environment to read one entry + if (!fChain) return -5; + Long64_t centry = fChain->LoadTree(entry); + if (centry < 0) return centry; + if (!fChain->InheritsFrom(TChain::Class())) return centry; + TChain *chain = (TChain*)fChain; + if (chain->GetTreeNumber() != fCurrent) { + fCurrent = chain->GetTreeNumber(); + Notify(); + } + return centry; +} + +void musrAnalysis::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fCurrent = -1; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("runID", &runID, &b_runID); + fChain->SetBranchAddress("eventID", &eventID, &b_eventID); + fChain->SetBranchAddress("weight", &weight, &b_weight); + fChain->SetBranchAddress("timeToNextEvent", &timeToNextEvent, &b_timeToNextEvent); + fChain->SetBranchAddress("BFieldAtDecay", &BFieldAtDecay_Bx, &b_BFieldAtDecay); + fChain->SetBranchAddress("muIniTime", &muIniTime, &b_muIniTime); + fChain->SetBranchAddress("muIniPosX", &muIniPosX, &b_muIniPosX); + fChain->SetBranchAddress("muIniPosY", &muIniPosY, &b_muIniPosY); + fChain->SetBranchAddress("muIniPosZ", &muIniPosZ, &b_muIniPosZ); + fChain->SetBranchAddress("muIniMomX", &muIniMomX, &b_muIniMomX); + fChain->SetBranchAddress("muIniMomY", &muIniMomY, &b_muIniMomY); + fChain->SetBranchAddress("muIniMomZ", &muIniMomZ, &b_muIniMomZ); + fChain->SetBranchAddress("muIniPolX", &muIniPolX, &b_muIniPolX); + fChain->SetBranchAddress("muIniPolY", &muIniPolY, &b_muIniPolY); + fChain->SetBranchAddress("muIniPolZ", &muIniPolZ, &b_muIniPolZ); + fChain->SetBranchAddress("muDecayDetID", &muDecayDetID, &b_muDecayDetID); + fChain->SetBranchAddress("muDecayPosX", &muDecayPosX, &b_muDecayPosX); + fChain->SetBranchAddress("muDecayPosY", &muDecayPosY, &b_muDecayPosY); + fChain->SetBranchAddress("muDecayPosZ", &muDecayPosZ, &b_muDecayPosZ); + fChain->SetBranchAddress("muDecayTime", &muDecayTime, &b_muDecayTime); + fChain->SetBranchAddress("muDecayPolX", &muDecayPolX, &b_muDecayPolX); + fChain->SetBranchAddress("muDecayPolY", &muDecayPolY, &b_muDecayPolY); + fChain->SetBranchAddress("muDecayPolZ", &muDecayPolZ, &b_muDecayPolZ); + fChain->SetBranchAddress("muTargetTime", &muTargetTime, &b_muTargetTime); + fChain->SetBranchAddress("muTargetPolX", &muTargetPolX, &b_muTargetPolX); + fChain->SetBranchAddress("muTargetPolY", &muTargetPolY, &b_muTargetPolY); + fChain->SetBranchAddress("muTargetPolZ", &muTargetPolZ, &b_muTargetPolZ); + fChain->SetBranchAddress("muTargetMomX", &muTargetMomX, &b_muTargetMomX); + fChain->SetBranchAddress("muTargetMomY", &muTargetMomY, &b_muTargetMomY); + fChain->SetBranchAddress("muTargetMomZ", &muTargetMomZ, &b_muTargetMomZ); + fChain->SetBranchAddress("muM0Time", &muM0Time, &b_muM0Time); + fChain->SetBranchAddress("muM0PolX", &muM0PolX, &b_muM0PolX); + fChain->SetBranchAddress("muM0PolY", &muM0PolY, &b_muM0PolY); + fChain->SetBranchAddress("muM0PolZ", &muM0PolZ, &b_muM0PolZ); + fChain->SetBranchAddress("posIniMomX", &posIniMomX, &b_posIniMomX); + fChain->SetBranchAddress("posIniMomY", &posIniMomY, &b_posIniMomY); + fChain->SetBranchAddress("posIniMomZ", &posIniMomZ, &b_posIniMomZ); + fChain->SetBranchAddress("nFieldNomVal", &nFieldNomVal, &b_nFieldNomVal); + fChain->SetBranchAddress("fieldNomVal", &fieldNomVal, &b_fieldNomVal); + fChain->SetBranchAddress("det_n", &det_n, &b_det_n); + fChain->SetBranchAddress("det_ID", det_ID, &b_det_ID); + fChain->SetBranchAddress("det_edep", det_edep, &b_det_edep); + fChain->SetBranchAddress("det_edep_el", det_edep_el, &b_det_edep_el); + fChain->SetBranchAddress("det_edep_pos", det_edep_pos, &b_det_edep_pos); + fChain->SetBranchAddress("det_edep_gam", det_edep_gam, &b_det_edep_gam); + fChain->SetBranchAddress("det_edep_mup", det_edep_mup, &b_det_edep_mup); + fChain->SetBranchAddress("det_nsteps", det_nsteps, &b_det_nsteps); + fChain->SetBranchAddress("det_length", det_length, &b_det_length); + fChain->SetBranchAddress("det_time_start", det_time_start, &b_det_time_start); + fChain->SetBranchAddress("det_time_end", det_time_end, &b_det_time_end); + fChain->SetBranchAddress("det_x", det_x, &b_det_x); + fChain->SetBranchAddress("det_y", det_y, &b_det_y); + fChain->SetBranchAddress("det_z", det_z, &b_det_z); + fChain->SetBranchAddress("det_kine", det_kine, &b_det_kine); + fChain->SetBranchAddress("det_VrtxKine", det_VrtxKine, &b_det_VrtxKine); + fChain->SetBranchAddress("det_VrtxX", det_VrtxX, &b_det_VrtxX); + fChain->SetBranchAddress("det_VrtxY", det_VrtxY, &b_det_VrtxY); + fChain->SetBranchAddress("det_VrtxZ", det_VrtxZ, &b_det_VrtxZ); + fChain->SetBranchAddress("det_VrtxVolID", det_VrtxVolID, &b_det_VrtxVolID); + fChain->SetBranchAddress("det_VrtxProcID", det_VrtxProcID, &b_det_VrtxProcID); + fChain->SetBranchAddress("det_VrtxTrackID", det_VrtxTrackID, &b_det_VrtxTrackID); + fChain->SetBranchAddress("det_VrtxParticleID", det_VrtxParticleID, &b_det_VrtxParticleID); + fChain->SetBranchAddress("det_VvvKine", det_VvvKine, &b_det_VvvKine); + fChain->SetBranchAddress("det_VvvX", det_VvvX, &b_det_VvvX); + fChain->SetBranchAddress("det_VvvY", det_VvvY, &b_det_VvvY); + fChain->SetBranchAddress("det_VvvZ", det_VvvZ, &b_det_VvvZ); + fChain->SetBranchAddress("det_VvvVolID", det_VvvVolID, &b_det_VvvVolID); + fChain->SetBranchAddress("det_VvvProcID", det_VvvProcID, &b_det_VvvProcID); + fChain->SetBranchAddress("det_VvvTrackID", det_VvvTrackID, &b_det_VvvTrackID); + fChain->SetBranchAddress("det_VvvParticleID", det_VvvParticleID, &b_det_VvvParticleID); + Notify(); +} + +Bool_t musrAnalysis::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void musrAnalysis::Show(Long64_t entry) +{ +// Print contents of entry. +// If entry is not specified, print current entry + if (!fChain) return; + fChain->Show(entry); +} +Int_t musrAnalysis::Cut(Long64_t entry) +{ +// This function may be called from Loop. +// returns 1 if entry is accepted. +// returns -1 otherwise. + return 1; +} +#endif // #ifdef musrAnalysis_cxx diff --git a/musrSimAna/musrCounter.cxx b/musrSimAna/musrCounter.cxx new file mode 100644 index 0000000..61eea88 --- /dev/null +++ b/musrSimAna/musrCounter.cxx @@ -0,0 +1,277 @@ +//#include +#include "musrCounter.hh" +#include "TCanvas.h" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +musrCounter::musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT) { + // pointerToAnalysis=this; + counterNr = CHANNEL_NR; + strcpy(counterName,CHANNEL_NAME); + counterType = CHANNEL_TYPE; + couterEThreshold = E_THRESH; + counterTimeShift = (Long64_t) TIME_SHIFT; + std::cout<<"musrCounter::musrCounter: Creating counter "<Fill(variable,vaha); +} + +//================================================================ +void musrCounter::DrawTDChistogram() { + char canvasName[501]; + sprintf(canvasName,"c%s",TDC_histoName); + TCanvas* cTmp = new TCanvas(canvasName,canvasName); + histTDC->Draw(); +} + +//================================================================ +void musrCounter::FillHitInCounter(Double_t edep, Long64_t timeBin, Long64_t timeBin2, Int_t kEntry, Int_t eveID, Int_t iDet, Int_t detectorID){ + //cDEL std::cout<<"FillHitInCounter: timeBin="<=couterEThreshold) { + + hitInfo* hInfo = new hitInfo(kEntry,eveID,iDet,detectorID,edep,timeBin2-counterTimeShift); + //cDEL std::cout<first)>(timeBinLimit+maxCoincidenceTimeWindow-counterTimeShift)) return; //note that maxCoincidenceTimeWindow is usually negative number + else { + // std::cout<<" Deleting hit from counter "<first; + // int tempEvnr= it->second; + // hitMap_TMP.insert( std::pair(tempBinT-timeBinsToRewind,tempEvnr) ); + hitInfo* tempEvnr= it->second; + tempEvnr->RewindTimeBin2(timeBinsToRewind); + hitMap_TMP.insert( std::pair(tempBinT-timeBinsToRewind,tempEvnr) ); + } + hitMap.swap(hitMap_TMP); +} + +//================================================================ +Bool_t musrCounter::IsInCoincidence(Long64_t timeBin, Bool_t ignoreHitsAtBinZero, Long64_t timeBinMinimum, Long64_t timeBinMaximum){ + // timeBin ... time bin, at which the coincidence is searched + // counterTimeShiftOfRequestingCounter ... time shift (in bin units) of the counter, for which the coincidence is searched + // ignoreHitsAtBinZero ... if "true", hits at timeBin will be ignored (needed for searching of coincidence of M counter + // with other M counters or P counters with other P counters) + // "false" should be used for coincidence detectors and vetos. + if (hitMap.empty()) return false; + // Long64_t timeBinToTest = timeBin + counterTimeShift - counterTimeShiftOfRequestingCounter; + Long64_t timeBinToTest = timeBin; + // If timeBinMinimum and timeBinMaximum are not specified, use internal time window of the detector (koincidence or veto detectors). + // Otherwise use timeBinMinimum and timeBinMaximum (e.g.coincidence of a positron counter with other positron counters). + Long64_t timeBinMin = (timeBinMinimum==-123456789) ? timeBinToTest + coincidenceTimeWindowMin : timeBinToTest + timeBinMinimum; + Long64_t timeBinMax = (timeBinMaximum==-123456789) ? timeBinToTest + coincidenceTimeWindowMax : timeBinToTest + timeBinMaximum; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + Long64_t timeBinOfCount_tmp = it->first; + if ((timeBinOfCount_tmp >= timeBinMin) && (timeBinOfCount_tmp <= timeBinMax)) { + if ((timeBin!=timeBinOfCount_tmp)||(!ignoreHitsAtBinZero)) { + return true; + } + } + else if (timeBinOfCount_tmp > timeBinMax) return false; + } + return false; +} + +//================================================================ +Bool_t musrCounter::GetNextGoodMuon(Int_t evtID, Long64_t timeBinMin, Long64_t& timeBinOfNextHit, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) { + // This function searches for a good muon, i.e. the muon + // 1) belongs to the currently analysed event + // 2) is in coincidence with all required coincidence detectors + // 3) is not in coincidence with veto detectors + // 4) is not in coincidence with other muons (more precisely - with hits in m-counter). + // INPUT PARAMETERS: evtID, timeBinMin + // OUTPUT PARAMETERS: timeBinOfNextHit + // + // Loop over the hits in the counter + // std::cout<<" musrCounter::GetNextGoodMuon timeBinMin="< S T O P !!!\n"; exit(1);} + doubleHitFound = false; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + + Long64_t timeBinOfCount_tmp = it->first; + timeBinOfNextHit = timeBinOfCount_tmp; + if (timeBinOfCount_tmp <= timeBinMin) continue; // This hit was already processed previously ==> skip it + + Int_t eventNumber = (it->second)->eventIDnumber; + if (eventNumber!=evtID) continue; // This trigger hit does not correspond to the currently processed event + // ==> skip it, becuase it was already proceesed or will be processed in future + numberOfMuonCandidates++; + // std::cout<<"*** "<second)->eventEntry; + idet = (it->second)->det_i; + idetID = (it->second)->det_id; + idetEdep = (it->second)->det_edep; + // timeBinOfNextHit = timeBinOfCount_tmp; + return true; + + endOfThisHit: + continue; + } + return false; +} + +//================================================================ +Bool_t musrCounter::GetNextGoodPositron(Int_t evtID, Long64_t timeBinMin, Long64_t timeBinMax, Long64_t& timeBinOfNextGoodHit, Long64_t& timeBinOfNextGoodHit_phaseShifted, Int_t& kEntry, Int_t& idet, Int_t& idetID, Double_t& idetEdep, Bool_t& doubleHitFound) { + // INPUT PARAMETERS: evtID, timeBinMin + // OUTPUT PARAMETERS: timeBinOfNextGoodHit + // + // Loop over the hits in the counter + if (hitMap.empty()) return false; + if (counterType!='P') {std::cout<<"\n!!! FATAL ERROR !!! musrCounter::GetNextGoodPositron: not the positron counter! ==> S T O P !!!\n"; exit(1);} + doubleHitFound = false; + for (hitMap_TYPE::iterator it = hitMap.begin(); it != hitMap.end(); ++it) { + Int_t eventNumber = (it->second)->eventIDnumber; + Long64_t timeBinOfCount_tmp = it->first; + if ((timeBinOfCount_tmp <= timeBinMin) || (timeBinOfCount_tmp > timeBinMax)) continue; // This hit is out of the data interval ==> skip it + + // Hit candidate was found. Now check its coincidences and vetos + for (counterMapType::const_iterator itCounter = koincidenceCounterMap.begin(); itCounter!=koincidenceCounterMap.end(); ++itCounter) { + // if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) )) goto endOfThisHit; // no coincidence found ==> skip hit + if (!( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) )) goto endOfThisHit; // no coincidence found ==> skip hit + } + for (counterMapType::const_iterator itCounter = vetoCounterMap.begin(); itCounter!=vetoCounterMap.end(); ++itCounter) { + // if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift) ) goto endOfThisHit; // coincidence with veto found ==> skip hit + if ( (itCounter->second)->IsInCoincidence(timeBinOfCount_tmp) ) goto endOfThisHit; // coincidence with veto found ==> skip hit + } + // Check coincidences with other P counters + // Coincidences with other P counters must be checked in musrAnalysis.cxx + // if (this->IsInCoincidence(timeBinOfCount_tmp,counterTimeShift,true) ) { // coincidence with another P-counter hit ==> skip hit + // hovno + // ADD HERE THE CHECK THAT THE POSITRON IS FOUND WITHIN THE DATA INTERVAL, i.e. within . + // AND ALSO CHECK THAT THIS IS TRUE FOR THE COINCIDENCE WITH OTHER POSITRON COUNTERS!!! + if (this->IsInCoincidence(timeBinOfCount_tmp,true,timeBinMin,timeBinMax) ) { // coincidence with another P-counter hit ==> skip hit + doubleHitFound = true; + // std::cout<<"tttttttttttttttttttt doubleHitN="<second)->eventEntry; + idet = (it->second)->det_i; + idetID = (it->second)->det_id; + idetEdep = (it->second)->det_edep; + timeBinOfNextGoodHit = timeBinOfCount_tmp; + timeBinOfNextGoodHit_phaseShifted = (it->second) -> timeBin2; + //cDEL std::cout<<"timeBinOfNextGoodHit ="< S T O P !!!"<second)->SetMaxCoincidenceTimeWindow(maxCoinc); + (it->second)->SetCoincidenceTimeWindow(min,max); + } +} + +//================================================================ +void musrCounter::SetCoincidenceTimeWindowOfAllVetoDetectors(Long64_t maxCoinc, Long64_t min, Long64_t max) { + for (counterMapType::const_iterator it = vetoCounterMap.begin(); it!=vetoCounterMap.end(); ++it) { + if ( ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=0) && + ( ((it->second)->GetMaxCoincidenceTimeWindow()) !=maxCoinc) ) { + std::cout<<" !!!! ERROR SetCoincidenceTimeWindowOfAllVetoDetectors : coincidenceTimeWindow set multiple times! ==> S T O P !!!"<second)->SetMaxCoincidenceTimeWindow(maxCoinc); + (it->second)->SetCoincidenceTimeWindow(min,max); + } +} + +//================================================================ +void musrCounter::myPrintThisCounter(Int_t evtID) { + Bool_t eventMixing=false; + std::cout<<"musrCounter::myPrintThisCounter: counterNr="<first<<" "<<(it->second)->eventIDnumber<<","; + if (evtID != (it->second)->eventIDnumber) {eventMixing=true;} + } + if (eventMixing) {std::cout<<" Potential event mmmmmmmmixing";} + std::cout< +#include +#include +#include +#include +#include + +class hitInfo { +public: + hitInfo(Int_t kEntry, Int_t evID, Int_t deI, Int_t detectorID, Double_t deEDEP, Long64_t timeBIN2) {eventEntry=kEntry; eventIDnumber=evID; det_i=deI; det_id = detectorID; det_edep=deEDEP; timeBin2=timeBIN2;} + ~hitInfo() {} + void RewindTimeBin2(Long64_t timeBinsToRewind) {timeBin2-=timeBinsToRewind;} + Int_t eventEntry; + Int_t eventIDnumber; + Int_t det_i; + Int_t det_id; + Double_t det_edep; + Long64_t timeBin2; +}; + +class musrCounter { + public: + musrCounter(int CHANNEL_NR, char CHANNEL_NAME[200], char CHANNEL_TYPE, float E_THRESH, int TIME_SHIFT); + ~musrCounter(); + int GetCounterNr() {return counterNr;} + char GetCounterType() {return counterType;} + void SetCoincidenceCounter(musrCounter* c, int icNr) { + int cNr = abs(icNr); + if (icNr>0) { + std::cout<<"SetCoincidenceCounter: Adding counter ="< counterMapType; + counterMapType koincidenceCounterMap; + counterMapType vetoCounterMap; + char TDC_histoName[200]; + int TDC_t0, TDC_t1, TDC_t2, TDC_histoNrAdd; + char TDC_histoNameAdd[200]; + TH1D* histTDC; + Long64_t coincidenceTimeWindowMin; + Long64_t coincidenceTimeWindowMax; + Long64_t maxCoincidenceTimeWindow; + // typedef std::map hitMap_TYPE; // Long64_t = timeBin, Int_t=eventID + typedef std::map hitMap_TYPE; // Long64_t = timeBin, hitInfo = eventID and det_i + hitMap_TYPE hitMap; + // std::list timeOfHitsList; + + Int_t doubleHitN; + Long64_t numberOfMuonCandidates; +}; + +#endif diff --git a/musrSimAna/musrSimAna.cxx b/musrSimAna/musrSimAna.cxx new file mode 100644 index 0000000..90ebd5b --- /dev/null +++ b/musrSimAna/musrSimAna.cxx @@ -0,0 +1,106 @@ +#include "musrAnalysis.hh" +// ROOT libraries: +#include +#include +#include +#include + +int main(int argc,char** argv) { + // if (strcmp(argv[2],"nographics")!=0) { + // TApplication class is used for the ROOT graphical interface + TApplication* myapp=new TApplication("myapp",0,0); + // std::cout<<"1 argc = "< Stop!"< Stop!"<Get("t1"); + if (tree==NULL) { + std::cout<<"\n Error!!! Tree \"t1\" not found! ==> Stop!"<4) { + for (int i=3; iRun()"<Run(kTRUE); + // std::cout<<" after myapp->Run()"<=4)&&(strcmp(argv[3],"nographic")==0)) {;} + else { + std::cout<<" argc = "<Run(kTRUE); + } + std::cout<<" after myapp->Run()"<Sumw2(); + histArray2D[i] = hTemp; + // std::cout<<"histogram hTemp defined, name="<Sumw2(); + histArray1D[i] = hTemp; + // std::cout<<"histogram hTemp defined, name="<cd(); tempC->Update(); + char chopt[1000]; + strcpy(chopt,option); + Int_t i_first=0, i_last=musrAnalysis::nrConditions; if (idHist>=0) {i_first=idHist; i_last=idHist+1;} + for (Int_t i=i_first; i=0) {i_first=idHist; i_last=idHist+1;} + for (Int_t i=i_first; i::const_iterator it = drawList.begin(); it != drawList.end(); ++it) { + int iii = *it; + std::cout<<" Drawing histog"<::const_iterator it = drawList.begin(); it != drawList.end(); ++it) { + int iii = *it; + DrawTH2D(option, iii); + } +} +//============================================================================================== +void musrTH::SetBinLabel1D(Int_t iBin,std::string slabel) { + for (Int_t i=0; iGetXaxis()->SetBinLabel(iBin,slabel.c_str()); + } +} +//============================================================================================== +Double_t musrTH::GetBinContent1D(Int_t i, Int_t jBin) { + // std::cout<<"musrTH::GetBinContent1D: i="<first; + Int_t iDecayHistoBin = it->second; + // if (iDecayHistoBin==-123456789) {iHumanBinForAllTheRest=iHumanBin; continue;} + Double_t value = decayMapHistos->GetBinContent1D(i,iDecayHistoBin); + Double_t value2 = histArray1D[i]-> GetBinContent(iHumanBin); + //decayMapHistos->GetBinContent1D(i,iDecayHistoBin); + histArray1D[i]-> SetBinContent(iHumanBin,value+value2); + } + + // // if (iHumanBinForAllTheRest != -1) { + for (Int_t j=1; j<= (decayMapHistos->GetNbinsX1D()); j++) { + Double_t value = decayMapHistos->GetBinContent1D(i,j); + Bool_t thisBinWasAssigned=false; + if (value!=0) { + for (humanDecayMultimapType::const_iterator it = myMultimap.begin(); it!=myMultimap.end(); ++it) { + Int_t iDecayHistoBin = it->second; + if (iDecayHistoBin==j) thisBinWasAssigned=true; + } + } + if ((!thisBinWasAssigned)&&(value!=0)) { + std::cout<<"musrHT.cxx: Some muons stoped and decayed in detector no. "< GetBinContent(XXXX); + } + } + + + } +} +//============================================================================================== +void musrTH::AssignFunction(TF1* function, Double_t xMin, Double_t xMax) { + funct = function; + funct_xMin = xMin; + funct_xMax = xMax; + +} +//============================================================================================== +void musrTH::ListHistograms() { + std::cout<<"musrTH::ListHistograms():"<GetName()<<"\", funct_xMin="<FixParameter(0,omega);} + if (strcmp(funct->GetName(),"funct2")==0) {funct->FixParameter(0,omega);} + if (strcmp(funct->GetName(),"funct3")==0) {funct->FixParameter(0,omega);} + if (strcmp(funct->GetName(),"funct4")==0) {funct->FixParameter(0,omega);} + if (strcmp(funct->GetName(),"rotFrameTime20")==0) { + if (funct->GetParameter(0)==0) { + funct->SetParameter(0,omega); std::cout<<" FUNKCE rotFrameTime20"<<"omega initialy at "<GetNumberFreeParameters()); i++) std::cout<GetParameter(i)<<", "; + std::cout<Fit(funct,"LL","",funct_xMin,funct_xMax); + } + // if (strcmp(funct->GetName(),"simpleExpoPLUSconst")==0) { + // N0_FromLastFit=funct->GetParameter(0); std::cout<<"N0_FromLastFit="<GetBinWidth(2); std::cout<<"BinWidth_FromLastFit="<SetLineWidth(2); + // funct->SetLineColor(2); +} +//============================================================================================== +//void musrTH::FillHumanDecayArray(musrTH* decayMapHistos, const int nBins, const int* iBins) { +// for (Int_t i=0; iGetBinContent1D(i,jBin); +// histArray1D[i]-> SetBinContent(j+1,value); +// } +// } +//} +//============================================================================================== diff --git a/musrSimAna/musrTH.hh b/musrSimAna/musrTH.hh new file mode 100644 index 0000000..7f6b9c0 --- /dev/null +++ b/musrSimAna/musrTH.hh @@ -0,0 +1,62 @@ +#ifndef musrTH_h +#define musrTH_h 1 +#include +#include +#include +#include +#include +#include +//#include +#include +#include + +class musrAnalysis; + +class musrTH { +public: + typedef std::map humanDecayMapType; + typedef std::multimap humanDecayMultimapType; + + musrTH(char* dimension, char* histoName, char* histoTitle, Int_t nrOfBinsX, Double_t minBoundaryX, Double_t maxBoundaryX, Int_t nrOfBinsX, Double_t minBoundaryY, Double_t maxBoundaryY, Double_t* variableName1, Double_t* variableName2); + ~musrTH(); + void FillTH1D(Double_t vaha, Bool_t* cond); + void FillTH2D(Double_t vaha, Bool_t* cond); + void DrawTH1D(Option_t* option, Int_t idHist); + void DrawTH1DdrawList(Option_t* option); + void DrawTH2D(Option_t* option, Int_t idHist); + void DrawTH2DdrawList(Option_t* option); + void SetDrawListOfHistograms(int i); + Bool_t IsThisHistoNamed(char* someName) {return (strcmp(someName,histogramName)==0);} + Double_t GetBinContent1D(Int_t i, Int_t jBin); + void SetBinLabel1D(Int_t iBin,std::string slabel); + // void FillHumanDecayArray(musrTH* decayMapHistos, const int nBins, const int* iBins); + void FillHumanDecayArray(musrTH* decayMapHistos, humanDecayMapType myMap, humanDecayMultimapType myMultimap); + // TH1D** GetHistArray1D(std::string& varToBeFilled) {varToBeFilled = variableToBeFilled_X; return histArray1D;} + TH1D** GetHistArray1D() {return histArray1D;} + // TH2D** GetHistArray2D() {return histArray2D;} + Int_t GetNbinsX1D(); + Int_t GetNbinsX2D(); + Int_t GetNbinsY2D(); + void AssignFunction(TF1* function, Double_t xMin, Double_t xMax); + void FitHistogramsIfRequired(Double_t omega); + void SetRotatingReferenceFrame(Double_t frequency, Double_t phase) {bool_rotating_reference_frame=true; + rot_ref_frequency=frequency; rot_ref_phase=phase;} + void ListHistograms(); + +private: + char histogramName[501]; + TH1D** histArray1D; + TH2D** histArray2D; + Double_t* variableToBeFilled_X; + Double_t* variableToBeFilled_Y; + std::list drawList; + TF1* funct; + Double_t funct_xMin, funct_xMax; + Bool_t bool_rotating_reference_frame; + Double_t rot_ref_frequency, rot_ref_phase; + // Double_t N0_FromLastFit; + // Double_t BinWidth_FromLastFit; + // Double_t funct_p0, funct_p1, funct_p2, funct_p3, funct_p4, funct_p5, funct_p6, funct_p7, funct_p8, funct_p9; +}; + +#endif