Kamil Sedlak
f6ccd6cc75
Added new subdirectory (analysis package) "musrSimAna". This is a general analysis program that allows one to plot histograms out of the simulated trees. The detector setup of the veto and coincidence relations with the positron/muon counters is defined in *.v1190 files. Also the histograms that should be plotted/saved are defined in this file. Thus (at least in principle) the user does not need to write an analysis program dedicated for his/her instrument. One just needs to modify the *.v1190 file. The "musrSimAna" is still in the development phase, no documentation is available at this stage.
1031 lines
47 KiB
C++
1031 lines
47 KiB
C++
#define musrAnalysis_cxx
|
|
#include "musrAnalysis.hh"
|
|
#include "musrTH.hh"
|
|
//#include <TH2.h>
|
|
#include <TStyle.h>
|
|
#include <TCanvas.h>
|
|
#include <algorithm>
|
|
#include <TF1.h>
|
|
//#include <math.h>
|
|
|
|
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; jentry<nentries;jentry++) {
|
|
Long64_t ientry = LoadTree(jentry);
|
|
if (ientry < 0) break;
|
|
nb = fChain->GetEntry(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()"<<std::endl;
|
|
|
|
// std::string V1190_FileName = "2.v1190";
|
|
// FILE *fSteeringFile=fopen(V1190_FileName.c_str(),"r");
|
|
FILE *fSteeringFile=fopen(charV1190FileName,"r");
|
|
if (fSteeringFile==NULL) {
|
|
if (strcmp(charV1190FileName,"Unset")==0) {
|
|
std::cout<<"musrDetectorConstruction: No input macro file specified"<<std::endl;
|
|
}
|
|
else {
|
|
std::cout << "E R R O R : Failed to open detector configuration file " << charV1190FileName << std::endl;
|
|
}
|
|
std::cout << "S T O P F O R C E D" << std::endl;
|
|
exit(1);
|
|
}
|
|
std::cout << "Configuration file \"" << charV1190FileName << "\" was opened."<<std::endl;
|
|
char line[1001];
|
|
|
|
// Write out the configuration file into the output file:
|
|
std::cout << "\n\n....oooOO0OOooo........oooOO0OOooo......Configuration file used for this run....oooOO0OOooo........oooOO0OOooo......"<<std::endl;
|
|
while (!feof(fSteeringFile)) {
|
|
fgets(line,1000,fSteeringFile);
|
|
// if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) // TS: Do not print comments
|
|
std::cout << line;
|
|
}
|
|
std::cout <<"....oooOO0OOooo........oooOO0OOooo......End of the configuration file.....oooOO0OOooo........oooOO0OOooo......\n\n"<<std::endl;
|
|
|
|
// Loop
|
|
rewind (fSteeringFile);
|
|
instrument = "Unset";
|
|
description = "Unset";
|
|
tdchwtype = "Unset";
|
|
tdcresolution = 195.3125*picosecond;
|
|
mdelay = 0;
|
|
pdelay = 800;
|
|
mcoincwin = 40;
|
|
pcoincwin = 40;
|
|
vcoincwin = 80;
|
|
muonRateFactor = 1.;
|
|
rewindTimeBins = 1000000000;
|
|
dataWindowMin = -0.2;
|
|
dataWindowMax = 10.0;
|
|
pileupWindowMin = -10.5;
|
|
pileupWindowMax = 10.5;
|
|
musrMode = 'D';
|
|
overallBinDelay = 0;
|
|
boolInfinitelyLowMuonRate = false;
|
|
// globTime = 0.;
|
|
|
|
int tmpivar;
|
|
float tmpfvar;
|
|
// int tmpCoincidenceMultimap[maxChannels][maxChannels];
|
|
std::multimap<int, int> tmpCoincidenceMultimap;
|
|
// for (int i=0; i<maxChannels; i++) {
|
|
// for (int j=0; j<maxChannels; j++) {
|
|
// // tmpCoincidenceMultimap[i][j]=0;
|
|
// tmpCoincidenceMultimap.insert(pair<int, int>("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 "<<tdcresolution/picosecond<<" ps)"<<std::endl;
|
|
std::cout<<" ===> too low time resolution might lead to numerical problems in musrCounter ==> S T O P"<<std::endl;
|
|
exit(1);
|
|
}
|
|
}
|
|
else if (strncmp(tmpString0,"MDELAY=",strlen("MDELAY="))==0) {
|
|
sscanf(&line[strlen("MDELAY=")],"%d",&tmpivar);
|
|
mdelay = tmpivar;
|
|
}
|
|
else if (strncmp(tmpString0,"PDELAY=",strlen("PDELAY="))==0) {
|
|
sscanf(&line[strlen("PDELAY=")],"%d",&tmpivar);
|
|
pdelay = tmpivar;
|
|
}
|
|
else if (strncmp(tmpString0,"MCOINCIDENCEW=",strlen("MCOINCIDENCEW="))==0) {
|
|
sscanf(&line[strlen("MCOINCIDENCEW=")],"%d",&tmpivar);
|
|
mcoincwin = tmpivar;
|
|
}
|
|
else if (strncmp(tmpString0,"PCOINCIDENCEW=",strlen("PCOINCIDENCEW="))==0) {
|
|
sscanf(&line[strlen("PCOINCIDENCEW=")],"%d",&tmpivar);
|
|
pcoincwin = tmpivar;
|
|
}
|
|
else if (strncmp(tmpString0,"VCOINCIDENCEW=",strlen("VCOINCIDENCEW="))==0) {
|
|
sscanf(&line[strlen("VCOINCIDENCEW=")],"%d",&tmpivar);
|
|
vcoincwin = tmpivar;
|
|
}
|
|
else if (strncmp(tmpString0,"MUONRATEFACTOR=",strlen("MUONRATEFACTOR="))==0) {
|
|
sscanf(&line[strlen("MUONRATEFACTOR=")],"%g",&tmpfvar);
|
|
muonRateFactor = tmpfvar;
|
|
}
|
|
else if (strncmp(tmpString0,"INFINITELYLOWMUONRATE=",strlen("INFINITELYLOWMUONRATE"))==0) {
|
|
boolInfinitelyLowMuonRate=true;
|
|
}
|
|
else if (strncmp(tmpString0,"REWINDTIMEBINS=",strlen("REWINDTIMEBINS="))==0) {
|
|
sscanf(&line[strlen("REWINDTIMEBINS=")],"%d",&tmpivar);
|
|
rewindTimeBins = tmpivar; std::cout<<" JUHELAK: rewindTimeBins="<<rewindTimeBins<<std::endl;
|
|
}
|
|
else if (strncmp(tmpString0,"DATAWINDOWMIN=",strlen("DATAWINDOWMIN="))==0) {
|
|
sscanf(&line[strlen("DATAWINDOWMIN=")],"%g",&tmpfvar);
|
|
dataWindowMin = tmpfvar;
|
|
}
|
|
else if (strncmp(tmpString0,"DATAWINDOWMAX=",strlen("DATAWINDOWMAX="))==0) {
|
|
sscanf(&line[strlen("DATAWINDOWMAX=")],"%g",&tmpfvar);
|
|
dataWindowMax = tmpfvar;
|
|
}
|
|
else if (strncmp(tmpString0,"PILEUPWINDOWMIN=",strlen("PILEUPWINDOWMIN="))==0) {
|
|
sscanf(&line[strlen("PILEUPWINDOWMIN=")],"%g",&tmpfvar);
|
|
pileupWindowMin = tmpfvar;
|
|
}
|
|
else if (strncmp(tmpString0,"PILEUPWINDOWMAX=",strlen("PILEUPWINDOWMAX="))==0) {
|
|
sscanf(&line[strlen("PILEUPWINDOWMAX=")],"%g",&tmpfvar);
|
|
pileupWindowMax = tmpfvar;
|
|
}
|
|
else if (strncmp(tmpString0,"MUSRMODE=",strlen("MUSRMODE="))==0) {
|
|
sscanf(&line[strlen("MUSRMODE=")],"%s",&tmpString1);
|
|
musrMode = tmpString1[0];
|
|
}
|
|
else if (strncmp(tmpString0,"musrTH",strlen("musrTH"))==0) {
|
|
// Definition of the histograms - either musrTH1D or musrTH2D
|
|
int beginningOfHistoTitle=0, endOfHistoTitle =0;
|
|
for (Int_t i=0; i<(int)strlen(line); i++) {
|
|
if (line[i]=='"') {
|
|
if (beginningOfHistoTitle==0) beginningOfHistoTitle=i+1;
|
|
else endOfHistoTitle=i-1;
|
|
}
|
|
}
|
|
char histoName[200];
|
|
sscanf(&line[0],"%*s %s",histoName);
|
|
char histoTitle[501];
|
|
CopySubstring(line,beginningOfHistoTitle,endOfHistoTitle,histoTitle);
|
|
int nrBinsX, nrBinsY;
|
|
float minX, maxX, minY, maxY;
|
|
char varXName[500]; char varYName[500];
|
|
char furtherOption[200];
|
|
float rot_ref_frequency, rot_ref_phase;
|
|
if (strcmp(tmpString0,"musrTH1D")==0) {
|
|
sscanf(&line[endOfHistoTitle+2],"%d %g %g %s %s %g %g",&nrBinsX,&minX,&maxX,varXName,furtherOption,&rot_ref_frequency,&rot_ref_phase);
|
|
}
|
|
else if (strcmp(tmpString0,"musrTH2D")==0) {
|
|
sscanf(&line[endOfHistoTitle+2],"%d %g %g %d %g %g %s %s",&nrBinsX,&minX,&maxX,&nrBinsY,&minY,&maxY,varXName,varYName);
|
|
}
|
|
else {
|
|
std::cout<<" !!! ERROR: keyword \""<<tmpString0<<"\" in the steering file not understood ==> S T O P !!!"<<std::endl;
|
|
exit(1);
|
|
}
|
|
std::string stringVariableName1 = varXName;
|
|
std::string stringVariableName2 = varYName;
|
|
variableMapType::iterator iter1 = variableMap.find(stringVariableName1);
|
|
variableMapType::iterator iter2;
|
|
if (iter1==variableMap.end()) {
|
|
std::cout<<" !!! ERROR: Variable with the name \""<<stringVariableName1<<"\" not defined in variableMap ==> ignored, histogram not defined! S T O P !!!"<<std::endl;
|
|
exit(1);
|
|
}
|
|
if (strcmp(tmpString0,"musrTH2D")==0) {
|
|
iter2 = variableMap.find(stringVariableName2);
|
|
if (iter2==variableMap.end()) {
|
|
std::cout<<" !!! ERROR: Variable with the name \""<<stringVariableName2<<"\" not defined in variableMap ==> ignored, histogram not defined! S T O P !!!"<<std::endl;
|
|
exit(1);
|
|
}
|
|
}
|
|
|
|
Double_t* pointerToVariable1 = iter1->second;
|
|
if (strcmp(tmpString0,"musrTH1D")==0) {
|
|
// std::cout<<"pointerToVariable1="<<pointerToVariable1<<" pointer to muDecayPosZ="<<&muDecayPosZ<<std::endl;
|
|
musrTH* myTH = new musrTH("1D",histoName,histoTitle,nrBinsX,minX,maxX,nrBinsY,minY,maxY,pointerToVariable1,NULL);
|
|
listOfAllHistograms1D.push_back(myTH);
|
|
if (strcmp(furtherOption,"rotreference")==0) {
|
|
myTH->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="<<motherHistoName<<")"<<std::endl;
|
|
for(std::list<musrTH*>::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!"<<std::endl;
|
|
std::cout<<" "<<line<<std::endl;
|
|
std::cout<<" ==> S T O P "<<std::endl;
|
|
exit(1);
|
|
}
|
|
if (motherOfHumanDecayPileupHistograms==NULL) {
|
|
std::cout<<"\n\nmusrAnalysis::ReadInInputParameters Mother histogram (motherOfHumanDecayPileupHistograms)"
|
|
<<"of the \"Pileup human decay histogram\" not found!"<<std::endl;
|
|
std::cout<<" "<<line<<std::endl;
|
|
std::cout<<" ==> S T O P "<<std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
char* pch0 = strstr(line,motherHistoName)+strlen(motherHistoName);
|
|
char* pch = strstr(pch0,motherHistoPileupName)+strlen(motherHistoPileupName);
|
|
|
|
int nscan; int N1; char NAME1[100];
|
|
nscan = sscanf(pch,"%d %s",&N1,NAME1);
|
|
// indexHuman=0;
|
|
Int_t iBinHuman=0;
|
|
Int_t nBinHuman=0;
|
|
do {
|
|
Bool_t nameDefined=false;
|
|
for (humanDecayMapType::const_iterator it = humanDecayMap.begin(); it!=humanDecayMap.end(); ++it) {
|
|
if ((it->second)==(std::string) NAME1) {
|
|
nameDefined = true;
|
|
iBinHuman = it->first;
|
|
}
|
|
}
|
|
if (!nameDefined) {
|
|
nBinHuman++;
|
|
iBinHuman=nBinHuman;
|
|
humanDecayMap.insert(std::pair<Int_t,std::string>(iBinHuman,NAME1));
|
|
}
|
|
humanDecayMultimap.insert(std::pair<Int_t,Int_t>(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="<<nrBinsHumanHist<<std::endl;
|
|
// humanDecayHistograms = new musrTH("1D","humanDecayHistograms","Where the muons stop and decay",indexHuman,0,float(indexHuman),0,0,0,NULL,NULL);
|
|
humanDecayHistograms = new musrTH("1D","humanDecayHistograms","Where the muons stop and decay",nrBinsHumanHist,0,float(nrBinsHumanHist),0,0,0,NULL,NULL);
|
|
humanDecayPileupHistograms = new musrTH("1D","humanDecayPileupHistograms","Where the muons contr. to pileup stop and decay",nrBinsHumanHist,0,float(nrBinsHumanHist),0,0,0,NULL,NULL);
|
|
// for (Int_t k=0; k<indexHuman; k++) {
|
|
// humanDecayHistograms->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 "<<iConditionTMP<<")"<< std::endl;
|
|
exit(1);
|
|
}
|
|
if (iConditionTMP>=nrConditions) {
|
|
std::cout<<" !!! ERROR: Condition number must be < "<<nrConditions<<" (not "<<iConditionTMP<<")"<< std::endl;
|
|
exit(1);
|
|
}
|
|
if (strcmp(conditionNameTMP,"alwaysTrue")==0) conditionMap[iConditionTMP]=&alwaysTrue;
|
|
else if (strcmp(conditionNameTMP,"oncePerEvent")==0) conditionMap[iConditionTMP]=&oncePerEvent;
|
|
else if (strcmp(conditionNameTMP,"muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&muonDecayedInSample_gen;
|
|
else if (strcmp(conditionNameTMP,"muonTriggered_gen")==0) conditionMap[iConditionTMP]=&muonTriggered_gen;
|
|
else if (strcmp(conditionNameTMP,"muonTriggered_det")==0) conditionMap[iConditionTMP]=&muonTriggered_det;
|
|
else if (strcmp(conditionNameTMP,"positronHit_det")==0) conditionMap[iConditionTMP]=&positronHit_det;
|
|
else if (strcmp(conditionNameTMP,"goodEvent_det")==0) conditionMap[iConditionTMP]=&goodEvent_det;
|
|
else if (strcmp(conditionNameTMP,"goodEvent_gen")==0) conditionMap[iConditionTMP]=&goodEvent_gen;
|
|
else if (strcmp(conditionNameTMP,"goodEvent_det_AND_goodEvent_gen")==0) conditionMap[iConditionTMP]=&goodEvent_det_AND_goodEvent_gen;
|
|
else if (strcmp(conditionNameTMP,"pileupEventCandidate")==0) conditionMap[iConditionTMP]=&pileupEventCandidate;
|
|
else if (strcmp(conditionNameTMP,"pileupEvent")==0) conditionMap[iConditionTMP]=&pileupEvent;
|
|
else if (strcmp(conditionNameTMP,"goodEvent_det_AND_muonDecayedInSample_gen")==0) conditionMap[iConditionTMP]=&goodEvent_det_AND_muonDecayedInSample_gen;
|
|
else {
|
|
std::cout<<" !!! ERROR: Condition of the name \""<<conditionNameTMP<<"\" not predefined ==> Add it in the musrAnalysis.cxx S T O P !!!"<<std::endl;
|
|
exit(1);
|
|
}
|
|
}
|
|
else if (strcmp(tmpString0,"draw")==0) {
|
|
// Which histograms should be drawn
|
|
int iConditionTMP;
|
|
char histoNameTMP[200];
|
|
sscanf(&line[0],"%*s %s %d",histoNameTMP,&iConditionTMP);
|
|
if (iConditionTMP<-1) {
|
|
std::cout<<" !!! ERROR: draw: condition number must be >= -1 (not "<<iConditionTMP<<")"<< std::endl;
|
|
exit(1);
|
|
}
|
|
if (iConditionTMP>=nrConditions) {
|
|
std::cout<<" !!! ERROR: draw: condition number must be < "<<nrConditions<<" (not "<<iConditionTMP<<")"<< std::endl;
|
|
exit(1);
|
|
}
|
|
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
|
|
if ((*it)->IsThisHistoNamed(histoNameTMP)) {
|
|
(*it)->SetDrawListOfHistograms(iConditionTMP);
|
|
}
|
|
}
|
|
for(std::list<musrTH*>::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<int,Double_t>(N1,PHASE_SHIFT/180*pi));
|
|
std::cout<<"N1="<<N1<<", PHASE_SHIFT="<<PHASE_SHIFT/180<<" pi"<<std::endl;
|
|
nscan = sscanf(pch,"%s %s",NAME1,NAME2);
|
|
char* pch1 = strstr(pch ,NAME1)+strlen(NAME1);
|
|
char* pch2 = strstr(pch1,NAME2)+strlen(NAME2);
|
|
pch=pch2;
|
|
nscan = sscanf(pch,"%d %g",&N1,&PHASE_SHIFT);
|
|
} while (nscan==2);
|
|
}
|
|
else if (strcmp(tmpString0,"fit")==0) {
|
|
char histoName[100]; char functionName[100]; float xMin; float xMax; float p0, p1, p2, p3, p4, p5, p6;
|
|
sscanf(&line[0],"%*s %s %s %g %g %g %g %g %g %g",histoName,functionName,&xMin,&xMax,&p0,&p1,&p2,&p3,&p4,&p5,&p6);
|
|
|
|
TF1 *funct;
|
|
|
|
if (strcmp(functionName,"funct1")==0) {
|
|
funct = new TF1("funct1","[3]*exp(([4]-x)/2.19703)*(1+[2]*cos(x*[0]+[1]))");
|
|
funct -> 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 \""<<functionName<<"\" not defined! ==> S T O P"<<std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
musrTH* tmpHistograms;
|
|
for(std::list<musrTH*>::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!"<<std::endl;
|
|
std::cout<<" "<<line<<std::endl;
|
|
std::cout<<" ==> S T O P "<<std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
tmpHistograms->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 ("<<tmpChannel<<") or channel coincidence required ("<<iCoinc
|
|
//ckdel 20.9.2010 <<") larger than maxChannels ("<<maxChannels<<") ==> S T O P"<<std::endl;
|
|
//ckdel 20.9.2010 exit(1);
|
|
//ckdel 20.9.2010 }
|
|
tmpCoincidenceMultimap.insert(std::pair<int, int>(tmpChannel,iCoinc));
|
|
// tmpCoincidenceMultimap[tmpChannel][abs(iCoinc)]=iCoinc;
|
|
} while (jj<endOfCoincidenceArray);
|
|
// Read in the last part of the configuration line (after the coincidence definition)
|
|
if (tmpType[0]=='P') {
|
|
char tmpHistoName[500]; strcpy(tmpHistoName, "Unset");
|
|
char tmpHistoNameAdd[500]; strcpy(tmpHistoNameAdd,"Unset");
|
|
int tmpT0=0, tmpT1=0, tmpT2=0, tmpHistoaddNr=0;
|
|
sscanf(&line[endOfCoincidenceArray],"%s %d %d %d %d %s",tmpHistoName,&tmpT0,&tmpT1,&tmpT2,&tmpHistoaddNr,tmpHistoNameAdd);
|
|
counter->SetTDChistogram(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="<<tdcresolution<<" ps"<<std::endl;
|
|
std::cout << "MDELAY="<<mdelay<<std::endl;
|
|
std::cout << "PDELAY="<<pdelay<<std::endl;
|
|
std::cout << "MCOINCIDENCEW="<<mcoincwin<<std::endl;
|
|
std::cout << "PCOINCIDENCEW="<<pcoincwin<<std::endl;
|
|
std::cout << "VCOINCIDENCEW="<<vcoincwin<<std::endl;
|
|
std::cout << "DATAWINDOWMIN="<<dataWindowMin<<std::endl;
|
|
std::cout << "DATAWINDOWMAX="<<dataWindowMax<<std::endl;
|
|
std::cout << "PILEUPWINDOWMIN="<<pileupWindowMin<<std::endl;
|
|
std::cout << "PILEUPWINDOWMAX="<<pileupWindowMax<<std::endl;
|
|
std::cout << "REWINDTIMEBINS="<<rewindTimeBins<<std::endl;
|
|
|
|
|
|
// Loop over all counters and set maximum coincidence window (maxCoincidenceTimeWindow):
|
|
// The "K" counters are tricky - they have to be set according to the detector, in which they
|
|
// are in coincidence.
|
|
|
|
pileupWindowBinMin = Long64_t(pileupWindowMin/tdcresolution);
|
|
pileupWindowBinMax = Long64_t(pileupWindowMax/tdcresolution);
|
|
dataWindowBinMin = Long64_t( dataWindowMin/tdcresolution);
|
|
dataWindowBinMax = Long64_t( dataWindowMax/tdcresolution);
|
|
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
|
|
char DetectorType = (it->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<int,int>::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 ("<<iChn2<<") not found! This should never happen!!! ==> S T O P"<<std::endl;
|
|
exit(1);
|
|
}
|
|
else {
|
|
(it->second)->SetCoincidenceCounter(itTwo->second,iChn2);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// for (int j=0; j<maxChannels; j++) {
|
|
// for (multimap<int,int>::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 ("<<j<<") not found! This should never happen!!! ==> S T O P"<<std::endl;
|
|
// exit(1);
|
|
// }
|
|
// else {
|
|
// // int iChn = tmpCoincidenceMultimap[iChanNr][j];
|
|
// int iChn = itCoinc->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()"<<std::endl;
|
|
// safeTimeWindow=2*(std::max(datawindow,pileupwindow));
|
|
safeTimeWindow = (dataWindowMax>pileupWindowMax)? 3*dataWindowMax: 3*pileupWindowMax;
|
|
std::cout<<"musrAnalysis::CreateHistograms(): safeTimeWindow = "<<safeTimeWindow<<" (bins="<<Long64_t(safeTimeWindow/tdcresolution)<<")"<<std::endl;
|
|
currentTime=0.;
|
|
// nextEventTime=0.;
|
|
nextUnfilledEventTime=0.;
|
|
numberOfRewinds=0;
|
|
numberOfGoodMuons=0;
|
|
lastPreprocessedEntry=-1;
|
|
// currentEventID =0;
|
|
fChain->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="<<fieldValue<<"T, omega="<<omega<<std::endl;
|
|
}
|
|
//================================================================
|
|
|
|
void musrAnalysis::SaveHistograms(char* runChar, char* v1190FileName) {
|
|
|
|
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
|
|
// (*it)->ListHistograms();
|
|
(*it)->FitHistogramsIfRequired(omega);
|
|
(*it)->DrawTH1DdrawList("");
|
|
// (*it)->DrawTH1D("hist",0);
|
|
// (*it)->DrawTH1D("hist",1);
|
|
}
|
|
for(std::list<musrTH*>::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()"<<std::endl;
|
|
TDirectory *dir = gDirectory;
|
|
char HistogramOutputFileName[500];
|
|
sprintf(HistogramOutputFileName,"data/his_%s_%s.root",runChar,v1190FileName);
|
|
TFile *fHistograms = new TFile(HistogramOutputFileName,"recreate");
|
|
TIter next(dir->GetList());
|
|
while ( TObject *obj = next() ) {
|
|
if (obj->InheritsFrom(TH1::Class()) ) {obj->Write();}
|
|
}
|
|
fHistograms->Close();
|
|
delete fHistograms;
|
|
|
|
|
|
//==============================
|
|
|
|
std::cout<<"==================================================================="<<std::endl;
|
|
std::cout<<"Number of raw trigger counts (ignoring pileup gate, vetos and coincidences): "<<numberOfMuonCandidates<<std::endl;
|
|
std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl;
|
|
Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
|
|
std::cout<<"Duration of the \"experiment\": "<<durationOfExperiment<<" second"<<std::endl;
|
|
// std::cout<<" (numberOfRewinds="<<numberOfRewinds<<")"<<std::endl;
|
|
std::cout<<"In this run, MUONRATEFACTOR was set to "<<muonRateFactor<<", and the raw trigger rate was "
|
|
<<numberOfMuonCandidates/durationOfExperiment<<" muons/second"<<std::endl;
|
|
std::cout<<" To get event rate of 30000 events/second, set MUONRATEFACTOR to "
|
|
<< muonRateFactor*( (numberOfMuonCandidates/durationOfExperiment)/30000 )<<std::endl;
|
|
std::cout<<"========================== E N D =================================="<<std::endl;
|
|
}
|
|
|
|
//================================================================
|
|
|
|
void musrAnalysis::AnalyseEvent(Long64_t iiiEntry) {
|
|
// std::cout<<"\n______________________________________________________________________________________________\n";
|
|
// std::cout<<"musrAnalysis::AnalyseEvent("<<iiiEntry<<") currentTime="<<currentTime<<" (bin="<<Long64_t(currentTime/tdcresolution)<<")"<<std::endl;
|
|
|
|
|
|
// Reset some initial variables
|
|
|
|
// Loop over several next event and preprocess them (i.e. fill
|
|
// them into the lists/maps of the class musrCounter).
|
|
while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
|
|
Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
|
|
nextUnfilledEventTime+=deltaT;
|
|
};
|
|
|
|
fChain->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<Long64_t>::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="<<timeBinLimit<<std::endl;
|
|
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
|
|
(*it).second->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()"<<std::endl;
|
|
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
|
|
(*it).second->RewindHitsInCounter(timeBinsToRewind);
|
|
}
|
|
}
|
|
|
|
//================================================================
|
|
|
|
void musrAnalysis::FillHistograms(Int_t iiiEntry) {
|
|
// std::cout<<"musrAnalysis::FillHistograms() event="<<eventID<<" , bool="<<generatedInfo<<","<<detectorInfo<<std::endl;
|
|
|
|
oncePerEvent = true;
|
|
Bool_t mCounterHitExistsForThisEventID = false;
|
|
Bool_t pCounterHitExistsForThisEventID = false;
|
|
Long64_t timeBinOfThePreviouslyProcessedHit = -100000000;
|
|
Long64_t timeBin0 = -100000000;
|
|
Long64_t timeBin1 = -100000000;
|
|
Long64_t timeBin2 = -100000000;
|
|
Int_t kEntry = 0;
|
|
Int_t posEntry = 0;
|
|
Int_t idetM = 0;
|
|
Int_t idetM_ID = 0;
|
|
Double_t idetM_edep = 0.;
|
|
Int_t idetP = 0;
|
|
Int_t idetP_ID = 0;
|
|
Double_t idetP_edep = 0.;
|
|
Bool_t doubleHitM = false;
|
|
Bool_t doubleHitP = false;
|
|
|
|
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
|
|
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
|
|
mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
|
|
timeBinOfThePreviouslyProcessedHit = timeBin0;
|
|
|
|
//___________________________________________________________
|
|
|
|
do {
|
|
// std::cout<<" timeBin0 ="<<timeBin0<<std::endl;
|
|
// Check whether there was good hit in the Positron counter
|
|
// Long64_t dataBinMin = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMin : timeBinOfThePreviouslyProcessedHit-100000000;
|
|
// Long64_t dataBinMax = (mCounterHitExistsForThisEventID) ? timeBin0+dataWindowBinMax : timeBinOfThePreviouslyProcessedHit+100000000;
|
|
pileup_muDecayDetID_double = -1001;
|
|
if (mCounterHitExistsForThisEventID) {
|
|
numberOfGoodMuons++;
|
|
Long64_t dataBinMin = timeBin0+dataWindowBinMin;
|
|
Long64_t dataBinMax = timeBin0+dataWindowBinMax;
|
|
pCounterHitExistsForThisEventID = PositronCounterHit(eventID,dataBinMin,dataBinMax,timeBin1,timeBin2,posEntry,idetP,idetP_ID,idetP_edep,doubleHitP);
|
|
//cDEL if (pCounterHitExistsForThisEventID) std::cout<<" timeBin1-timeBin2 = "<<timeBin1<<"-"<<timeBin2<<"="<<timeBin1-timeBin2 <<std::endl;
|
|
if (pCounterHitExistsForThisEventID && (kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) {
|
|
// This must be a pileup event (positron counter hit comes from the different event than the muon counter hit)
|
|
fChain->GetEntry(posEntry);
|
|
pileup_muDecayDetID_double = muDecayDetID;
|
|
// if (pileup_muDecayDetID_double==-1000) {
|
|
// std::cout<<"DEBUG: pileup_muDecayDetID_double==-1000, posEntry="<<posEntry<<", eventID="<<eventID<<", idetP_edep="<<idetP_edep<<", idetP="<<idetP<<", idetP_ID="<<idetP_ID<<std::endl;
|
|
// }
|
|
fChain->GetEntry(iiiEntry);
|
|
}
|
|
}
|
|
// if (doubleHitP) {
|
|
// std::cout<<"yyyyyyyyyyyyyyyyy Total double hit nr="<<testIVar1++<<std::endl;
|
|
// }
|
|
|
|
// CALCULATE VARIABLES
|
|
// Initial positron
|
|
pos_Trans_Momentum = sqrt(posIniMomX*posIniMomX+posIniMomY*posIniMomY);
|
|
pos_Momentum = sqrt(pos_Trans_Momentum*pos_Trans_Momentum+posIniMomZ*posIniMomZ);
|
|
pos_Radius = pos_Trans_Momentum/(-BFieldAtDecay_Bz)/0.3;
|
|
pos_Theta = myAtan2(pos_Trans_Momentum,posIniMomZ);
|
|
pos_Phi = myAtan2(posIniMomY,posIniMomX);
|
|
// Detector info
|
|
det_m0edep = (mCounterHitExistsForThisEventID) ? idetM_edep : -1000.;
|
|
// det_time0 = timeBin0*tdcresolution;
|
|
// det_time1 = timeBin1*tdcresolution;
|
|
det_time10 = ((timeBin0==-100000000)|| (timeBin1==-100000000)) ? -100000000 : (timeBin1-timeBin0)*tdcresolution;
|
|
det_time20 = ((timeBin0==-100000000)|| (timeBin2==-100000000)) ? -100000000 : (timeBin2-timeBin0)*tdcresolution;
|
|
//cks if ((det_time10>-1)&&(det_time10<0.01)) {std::cout<<" eventID="<<eventID<<" det_time10="<< det_time10<<std::endl;}
|
|
gen_time10 = muDecayTime-muM0Time;
|
|
det_time10_MINUS_gen_time10 = (det_time10 - gen_time10)/picosecond;
|
|
det_posEdep = (pCounterHitExistsForThisEventID) ? idetP_edep : -1000.;
|
|
// std::cout<<eventID<<" det_time10="<<det_time10<<" t1="<<timeBin1*tdcresolution<<" t0="<<timeBin0*tdcresolution<<" gen_time10="<<gen_time10<<std::endl;
|
|
|
|
// CALCULATE CONDITIONS
|
|
alwaysTrue = true;
|
|
muonDecayedInSample_gen = muDecayDetID==201;
|
|
muonTriggered_gen = muM0Time > -1000;
|
|
muonTriggered_det = mCounterHitExistsForThisEventID;
|
|
positronHit_det = pCounterHitExistsForThisEventID;
|
|
goodEvent_det = muonTriggered_det && positronHit_det;
|
|
goodEvent_gen = (muDecayTime>-999)&&(muM0Time>-999);
|
|
goodEvent_det_AND_goodEvent_gen = goodEvent_det && goodEvent_gen;
|
|
pileupEventCandidate = ((kEntry>0)&&(posEntry>0)&&(kEntry!=posEntry)) ? true:false;
|
|
pileupEvent = pileupEventCandidate&&goodEvent_det;
|
|
goodEvent_det_AND_muonDecayedInSample_gen = goodEvent_det && muonDecayedInSample_gen;
|
|
// if (pileupEvent) std::cout<<"pileupEvent (eventID="<<eventID<<std::endl;
|
|
|
|
// Fill pileup-variables, but only if positron comes from different muon than the trigger signal
|
|
// if (mCounterHitExistsForThisEventID&&pCounterHitExistsForThisEventID)
|
|
// std::cout<<eventID<<" "<<mCounterHitExistsForThisEventID<<": iiiEntry (trigger) = "<<iiiEntry<<" kEntry (muon hit)="<<kEntry<<" "<<pCounterHitExistsForThisEventID<<" posEntry (positron hit)="<<posEntry<<std::endl;
|
|
// if (goodEvent_det) {
|
|
// std::cout<<"wwwwwwwwwwwwww \nEvent Nr: "<<eventID<<" muDecayTime="<< muDecayTime <<" muM0Time="<<muM0Time<<" det_time_start[idetP]="<<det_time_start[idetP]<<" t1-t0="<<(timeBin1-timeBin0)*tdcresolution<<std::endl;
|
|
// }
|
|
for (Int_t i=0; i<nrConditions; i++) { condition[i]=false; }
|
|
// MyPrintConditions();
|
|
for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
|
|
int i = it->first;
|
|
condition[i]=*(it->second);
|
|
}
|
|
|
|
// A special condition can me added here:
|
|
// condition[29] = muIniMomZ>27;
|
|
|
|
// Fill all "Detector Histograms"
|
|
for(std::list<musrTH*>::const_iterator it = listOfAllHistograms1D.begin(); it != listOfAllHistograms1D.end(); ++it) {
|
|
(*it)->FillTH1D(wght,&condition[0]);
|
|
}
|
|
for(std::list<musrTH*>::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!!!"<<std::endl;
|
|
std::cout<<" idetP="<<idetP<<" idetP_ID="<<idetP_ID<< std::endl;
|
|
}
|
|
else{
|
|
(itPcounterHitInThisEvent->second)->FillTDChistogram(float(timeBin1-timeBin0+overallBinDelay)+0.5,wght);
|
|
}
|
|
}
|
|
|
|
// Check whether there is a good trigger in the next event
|
|
// std::cout<<" FillHistograms: timeBinOfThePreviouslyProcessedHit = "<<timeBinOfThePreviouslyProcessedHit<<std::endl;
|
|
// mCounterHitExistsForThisEventID = mCounter->GetNextGoodHitInThisEvent(eventID,timeBinOfThePreviouslyProcessedHit,0,'M',timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
|
|
mCounterHitExistsForThisEventID = mCounter->GetNextGoodMuon(eventID,timeBinOfThePreviouslyProcessedHit,timeBin0,kEntry,idetM,idetM_ID,idetM_edep,doubleHitM);
|
|
timeBinOfThePreviouslyProcessedHit = timeBin0;
|
|
// if (mCounterHitExistsForThisEventID) std::cout<<" YYYYYYYYYYYYYYYYYYY check this : LOOOPING AGAIN"<<std::endl;
|
|
oncePerEvent=false;
|
|
} while(mCounterHitExistsForThisEventID);
|
|
|
|
}
|
|
|
|
//================================================================
|
|
void musrAnalysis::InitialiseEvent() {
|
|
runID_double = runID;
|
|
eventID_double = eventID;
|
|
muDecayDetID_double = muDecayDetID;
|
|
det_n_double = det_n;
|
|
muDecayPosR = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
|
|
|
|
wght = ( fChain->GetBranch("weight") ) ? weight : 1.;
|
|
}
|
|
|
|
//================================================================
|
|
Double_t musrAnalysis::PreprocessEvent(Long64_t iEn) {
|
|
// std::cout<<"musrAnalysis::PreprocessEvent()"<<std::endl;
|
|
fChain->GetEntry(iEn); InitialiseEvent();
|
|
// std::cout<<"\n musrAnalysis::PreprocessEvent() Filling event "<<eventID<<std::endl;
|
|
// MyPrintTree();
|
|
|
|
Double_t globTime = nextUnfilledEventTime;
|
|
for (Int_t i=0; i<det_n; i++) {
|
|
// // Int_t detID=det_ID[i];
|
|
std::map<int,musrCounter*>::iterator it;
|
|
it = allCounterMap.find(det_ID[i]);
|
|
if (it==allCounterMap.end()) {
|
|
// uncomment std::cout<<"Active detector with det_ID="<<det_ID[i]<<" not found !!!"<<std::endl;
|
|
}
|
|
else {
|
|
// Double_t omega = 851.372*fieldNomVal[0];
|
|
Double_t dPhaseShift = (omega==0) ? 0:phaseShiftMap[det_ID[i]]/omega;
|
|
//cDEL if (det_ID[i]<20) std::cout<<"phaseShift = "<<phaseShiftMap[det_ID[i]]<<" dPhaseShift="<<dPhaseShift<<" tdcresolution="<<tdcresolution<<std::endl;
|
|
Long64_t timeBin = Long64_t( (globTime+det_time_start[i] )/tdcresolution );
|
|
Long64_t timeBin2 = Long64_t( (globTime+det_time_start[i]+dPhaseShift)/tdcresolution );
|
|
//cDEL if (det_ID[i]<20) std::cout<<" timeBin = "<<timeBin<<" ("<<globTime+det_time_start[i]<<")"<<std::endl;
|
|
//cDEL if (det_ID[i]<20) std::cout<<" timeBin2 = "<<timeBin2<<" ("<<globTime+det_time_start[i]+dPhaseShift<<")"<<std::endl;
|
|
(*it).second->FillHitInCounter(det_edep[i],timeBin,timeBin2,iEn,eventID,i,det_ID[i]);
|
|
}
|
|
}
|
|
|
|
// std::cout<<"lastPreprocessedEntry+1="<<lastPreprocessedEntry+1<<" iEn="<<iEn<<" eventID="<<eventID<<std::endl;
|
|
if ((lastPreprocessedEntry+1)!=iEn) {std::cout<<"ERROR PreprocessEvent - should never happen!!!"<<std::endl; }
|
|
lastPreprocessedEntry = iEn;
|
|
return timeToNextEvent*muonRateFactor;
|
|
}
|
|
|
|
//================================================================
|
|
Bool_t musrAnalysis::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) {
|
|
if (pCounterMap.empty()) return false;
|
|
|
|
Bool_t positronHitFound = false;
|
|
if (musrMode=='D') {
|
|
// Loop over all positron counters
|
|
for (counterMapType::const_iterator it = pCounterMap.begin(); it!=pCounterMap.end(); ++it) {
|
|
// std::cout<<" ===POSITRON==="<< pCounterMap.size()<<std::endl;
|
|
// Bool_t thereWasHit = (it->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="<<tBin1<<" tBin2="<<tBin2<<std::endl;
|
|
if (doubleHit) {return false;} // There were two hits in the same positron counter
|
|
if (thereWasHit) {
|
|
if (positronHitFound) {doubleHit = true; return false;} // There were two hits in two different positron counters
|
|
positronHitFound = true;
|
|
}
|
|
}
|
|
}
|
|
return positronHitFound;
|
|
}
|
|
|
|
//================================================================
|
|
void musrAnalysis::CopySubstring(char* inputChar, int iStart, int iEnd, char* outputChar) {
|
|
int iiiEnd = std::min(iEnd,(int)strlen(inputChar));
|
|
int lastChar=0;
|
|
for (Int_t i=0; i<=(iiiEnd-iStart); i++) {
|
|
outputChar[i]=inputChar[i+iStart];
|
|
lastChar=i;
|
|
}
|
|
outputChar[lastChar+1]='\0';
|
|
// std::cout<<"iStart="<<iStart<<" iiiEnd="<<iiiEnd<<std::endl;
|
|
}
|
|
|
|
//================================================================
|
|
void musrAnalysis::MyPrintTree() {
|
|
std::cout<<"==========\n det_n="<<det_n<<std::endl;
|
|
for (Int_t jj=0; jj<det_n; jj++) {
|
|
std::cout<<"\t det_ID="<<det_ID[jj]<<"\t det_edep="<<det_edep[jj]<<"\t det_time_start="<<det_time_start[jj]<<std::endl;
|
|
}
|
|
}
|
|
|
|
//================================================================
|
|
void musrAnalysis::MyPrintConditions() {
|
|
std::cout<<"CONDITION NR = ";
|
|
for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
|
|
std::cout<<" "<<it->first;
|
|
}
|
|
std::cout<<std::endl<<" ";
|
|
for (conditionMapType::const_iterator it = conditionMap.begin(); it!=conditionMap.end(); ++it) {
|
|
std::cout<<" "<<*(it->second);
|
|
}
|
|
std::cout<<std::endl;
|
|
}
|
|
//================================================================
|
|
Double_t musrAnalysis::myAtan2(Double_t y, Double_t x) {
|
|
if ( (y<-999.99999) && (y>-1000.00001) ) return -1000.;
|
|
if ( (x<-999.99999) && (x>-1000.00001) ) return -1000.;
|
|
if ( (x==0) && (y==0) ) return -1000.;
|
|
return atan2(y,x);
|
|
}
|