LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
musrsim/musrSimAna/musrAnalysis.cxx
Kamil Sedlak cb18486947 6.5.2011 Kamil Sedlak 
1) Adding/improving the simulation of light signals and APD
2) Many small changes and improvements
3) Adding manual to musrSimAna to the svn repository
4) Adding some example files for musrSim
2011-05-06 15:16:49 +00:00

1305 lines
62 KiB
C++
Raw Blame History

#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[10001];
// 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,10000,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,10000,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 (strcmp(tmpString0,"CLONECHANNEL")==0) {
int ichannel_orig_tmp, ichannel_new_tmp;
sscanf(&line[0],"%*s %d %d",&ichannel_orig_tmp,&ichannel_new_tmp);
bool_clonedChannelsMultimap_NotEmpty = true;
clonedChannelsMultimap.insert(std::pair<int,int>(ichannel_orig_tmp,ichannel_new_tmp));
}
else if (strcmp(tmpString0,"DEBUGEVENT")==0) {
int ieventToDebug_tmp, iLevelToDebug_tmp;
sscanf(&line[0],"%*s %d %d",&ieventToDebug_tmp,&iLevelToDebug_tmp);
bool_debugingRequired=true;
debugEventMap.insert(std::pair<int,int>(ieventToDebug_tmp,iLevelToDebug_tmp));
}
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 if (strcmp(conditionNameTMP,"goodEvent_F_det")==0) conditionMap[iConditionTMP]=&goodEvent_F_det;
else if (strcmp(conditionNameTMP,"goodEvent_B_det")==0) conditionMap[iConditionTMP]=&goodEvent_B_det;
else if (strcmp(conditionNameTMP,"goodEvent_U_det")==0) conditionMap[iConditionTMP]=&goodEvent_U_det;
else if (strcmp(conditionNameTMP,"goodEvent_D_det")==0) conditionMap[iConditionTMP]=&goodEvent_D_det;
else if (strcmp(conditionNameTMP,"goodEvent_L_det")==0) conditionMap[iConditionTMP]=&goodEvent_L_det;
else if (strcmp(conditionNameTMP,"goodEvent_R_det")==0) conditionMap[iConditionTMP]=&goodEvent_R_det;
else if (strcmp(conditionNameTMP,"goodEvent_F_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_F_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_B_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_B_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_U_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_U_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_D_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_D_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_L_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_L_det_AND_pileupEvent;
else if (strcmp(conditionNameTMP,"goodEvent_R_det_AND_pileupEvent")==0) conditionMap[iConditionTMP]=&goodEvent_R_det_AND_pileupEvent;
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,"counterGrouping")==0) {
int nscan; int N1; char NAME[100];
char *pch = line + strlen("counterGrouping");
char counterGroupName[100];
nscan = sscanf(pch,"%s",counterGroupName);
char* pch1 = strstr(pch,counterGroupName)+strlen(counterGroupName);
pch = pch1;
do {
nscan = sscanf(pch,"%d",&N1);
if (strcmp(counterGroupName,"F")==0) F_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"B")==0) B_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"U")==0) U_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"D")==0) D_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"L")==0) L_posCounterList.push_back(N1);
else if (strcmp(counterGroupName,"R")==0) R_posCounterList.push_back(N1);
else {
std::cout<<"\n\n UNKNOWN COUNTER GROUP REQUIRED !!! =====> S T O P F O R C E D"<<std::endl;
std::cout<<line<<std::endl;
exit(1);
}
std::cout<<"counterGroupName="<<counterGroupName<<" N1="<<N1<<std::endl;
nscan = sscanf(pch,"%s",NAME);
char* pch2 = strstr(pch ,NAME)+strlen(NAME);
pch=pch2;
nscan = sscanf(pch,"%d",&N1);
} while (nscan==1);
}
else if (strcmp(tmpString0,"setSpecialAnticoincidenceTimeWindow")==0) {
int channelNumber;
double min, max;
char unit[100]="Undefined";
Long64_t lmin, lmax;
sscanf(&line[0],"%*s %d %lf %lf %s",&channelNumber,&min,&max,unit);
if ((strcmp(unit,"bin")==0)||(strcmp(unit,"bins")==0)) {lmin = Long64_t (min); lmax = Long64_t (max);}
else if ((strcmp(unit,"ns") ==0)||(strcmp(unit,"nanosecond") ==0)) {
lmin = Long64_t (min*nanosecond/tdcresolution);
lmax = Long64_t (max*nanosecond/tdcresolution);
}
else if ((strcmp(unit,"mus") ==0)||(strcmp(unit,"microsecond") ==0)) {
lmin = Long64_t (min*microsecond/tdcresolution);
lmax = Long64_t (max*microsecond/tdcresolution);
}
else {
std::cout<<" ERROR! musrAnalysis: unknown unit in setSpecialAnticoincidenceTimeWindow !!! (unit="<<unit<<")"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
std::cout<<" lmin="<<lmin<<" bins"<<std::endl;
std::cout<<" lmax="<<lmax<<" bins"<<std::endl;
if (allCounterMap.find(channelNumber)==allCounterMap.end()) {
std::cout<<" ERROR! musrAnalysis: unknown counter ("<<channelNumber<<") in setSpecialAnticoincidenceTimeWindow"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
musrCounter* counter = (allCounterMap.find(channelNumber))->second;
counter->SetAntiCoincidenceTimeWindowMin(lmin);
counter->SetAntiCoincidenceTimeWindowMax(lmax);
}
else if (strcmp(tmpString0,"artificiallyChangeMuDecayTime")==0) {
float min, max, mmmin, mmmax;
sscanf(&line[0],"%*s %g %g %g %g %g %g %g",&min,&max,&mmmin,&mmmax);
bool_muDecayTimeTransformation = true;
muDecayTime_t_min = min;
muDecayTime_t_max = max;
muDecayTime_Transformation_min = mmmin;
muDecayTime_Transformation_max = mmmax;
if ((muDecayTime_t_max <= muDecayTime_t_min) || (muDecayTime_Transformation_max <= muDecayTime_Transformation_min)) {
std::cout<<" ERROR! musrAnalysis: error when setting the \"artificiallyChangeMuDecayTime\" parameters! Min > Max !!!"<<std::endl;
std::cout<<" ==> S T O P "<<std::endl;
exit(1);
}
}
else if (strcmp(tmpString0,"fit")==0) {
char histoName[100]; char functionName[100]; char functOption[100]; float xMin; float xMax; float p0, p1, p2, p3, p4, p5, p6;
sscanf(&line[0],"%*s %s %s %s %g %g %g %g %g %g %g",histoName,functionName,functOption,&xMin,&xMax,&p0,&p1,&p2,&p3,&p4,&p5,&p6);
if (strcmp(functOption,"\"\"")==0) strcpy(functOption,"");
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 = new TF1("rotFrameTime20","[2]*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, functOption, 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) {
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 (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_M(pileupWindowBinMin,pileupWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('M',-mcoincwin+pileupWindowBinMin,-mcoincwin,mcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-mcoincwin+pileupWindowBinMin,-vcoincwin,vcoincwin);
}
else if (DetectorType=='P') {
(it->second)->SetMaxCoincidenceTimeWindow(dataWindowBinMin);
(it->second)->SetCoincidenceTimeWindow_P(dataWindowBinMin,dataWindowBinMax);
(it->second)->SetCoincidenceTimeWindowOfAllCoincidenceDetectors('P',-pcoincwin+dataWindowBinMin,-pcoincwin,pcoincwin);
(it->second)->SetCoincidenceTimeWindowOfAllVetoDetectors(-pcoincwin+dataWindowBinMin,-vcoincwin,vcoincwin);
}
// else if (DetectorType=='V') {
// (it->second)->SetMaxCoincidenceTimeWindow(-vcoincwin);
// (it->second)->SetCoincidenceTimeWindow(-vcoincwin,vcoincwin);
// }
}
// 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();
Long64_t numberOfMuonCandidatesAfterVK = mCounter->GetNumberOfMuonCandidatesAfterVK();
Double_t durationOfExperiment = (numberOfRewinds*rewindTimeBins*tdcresolution+currentTime)/1000000;
hInfo->Fill(5,(Double_t) numberOfMuonCandidates); // number of "triggers"
hInfo->Fill(4,(Double_t) numberOfMuonCandidatesAfterVK); // number of "triggers"
hInfo->Fill(3,durationOfExperiment);
//==============================
// 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 trigger counts (after vetos and coincidences but not pile-up rejected): "<<numberOfMuonCandidatesAfterVK<<std::endl;
std::cout<<"Number of triggered events (i.e. only good \"muons\"): "<<numberOfGoodMuons<<std::endl;
std::cout<<"Duration of the \"experiment\": "<<durationOfExperiment<<" second"<<std::endl;
// std::cout<<" (numberOfRewinds="<<numberOfRewinds<<")"<<std::endl;
std::cout<<"MUONRATEFACTOR was set to "<<muonRateFactor<<std::endl;
std::cout<<" Raw trigger rate = "<<numberOfMuonCandidates/durationOfExperiment<<" muons/second";
std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = "
<< muonRateFactor*( (numberOfMuonCandidates/durationOfExperiment)/30000 )<<")"<<std::endl;
std::cout<<" Trigger rate after V & K = "<<numberOfMuonCandidatesAfterVK/durationOfExperiment<<" muons/second";
std::cout<<" (to get 30000 events/second, set MUONRATEFACTOR = "
<< muonRateFactor*( (numberOfMuonCandidatesAfterVK/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).
if (bool_debugingRequired) {
if (debugEventMap[eventID]>0) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"PreprocessEvent\"_________"<<std::endl;}
}
while (((iiiEntry>lastPreprocessedEntry)||(((nextUnfilledEventTime-currentTime)<safeTimeWindow))&&(!boolInfinitelyLowMuonRate)) && (lastPreprocessedEntry+1<nentries)) {
Double_t deltaT = PreprocessEvent(lastPreprocessedEntry+1);
nextUnfilledEventTime+=deltaT;
};
fChain->GetEntry(iiiEntry); InitialiseEvent();
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) PrintHitsInAllCounters();
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"PreprocessEvent\"_________"<<std::endl;}
}
// 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);
// }
if (bool_debugingRequired) {
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(before \"FillHistograms\"_________"<<std::endl;}
}
FillHistograms(iiiEntry);
if (bool_debugingRequired) {
if (debugEventMap[eventID]>1) {std::cout<<"DEBUGEVENT "<<eventID<<"_________________(after \"FillHistograms\"_________"<<std::endl;}
}
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::PrintHitsInAllCounters() {
std::cout<<"___________________\n";
for (counterMapType::const_iterator it = allCounterMap.begin(); it!=allCounterMap.end(); ++it) {
(*it).second->myPrintThisCounter(eventID,0);
}
}
//================================================================
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_eventID = -1001;
pileup_muDecayDetID_double = -1001;
pileup_muDecayPosX = -1000000000;
pileup_muDecayPosY = -1000000000;
pileup_muDecayPosZ = -1000000000;
pileup_muDecayPosR = -1000000000;
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_eventID = eventID;
pileup_muDecayDetID_double = muDecayDetID;
pileup_muDecayPosX = muDecayPosX;
pileup_muDecayPosY = muDecayPosY;
pileup_muDecayPosZ = muDecayPosZ;
pileup_muDecayPosR = sqrt(muDecayPosX*muDecayPosX+muDecayPosY*muDecayPosY);
// 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;
// posCounterList_Iterator = find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID);
// goodEvent_F_det = posCounterList_Iterator != F_posCounterList.end()
goodEvent_F_det = goodEvent_det && ( (find(F_posCounterList.begin(), F_posCounterList.end(), idetP_ID)) != F_posCounterList.end() );
goodEvent_B_det = goodEvent_det && ( (find(B_posCounterList.begin(), B_posCounterList.end(), idetP_ID)) != B_posCounterList.end() );
goodEvent_U_det = goodEvent_det && ( (find(U_posCounterList.begin(), U_posCounterList.end(), idetP_ID)) != U_posCounterList.end() );
goodEvent_D_det = goodEvent_det && ( (find(D_posCounterList.begin(), D_posCounterList.end(), idetP_ID)) != D_posCounterList.end() );
goodEvent_L_det = goodEvent_det && ( (find(L_posCounterList.begin(), L_posCounterList.end(), idetP_ID)) != L_posCounterList.end() );
goodEvent_R_det = goodEvent_det && ( (find(R_posCounterList.begin(), R_posCounterList.end(), idetP_ID)) != R_posCounterList.end() );
// std::cout<<"goodEvent_F_det="<<goodEvent_F_det<<std::endl;
if (pileupEvent&&goodEvent_F_det) {
// std::cout<<" DEBUG: Pileup Event: eventID = "<<eventID<<" pileup_eventID = "<<pileup_eventID<<" det_time10 = "<<det_time10<<std::endl;
// debugEventMap.insert(std::pair<int,int>(eventID,10));
}
goodEvent_F_det_AND_pileupEvent = goodEvent_F_det && pileupEvent;
goodEvent_B_det_AND_pileupEvent = goodEvent_B_det && pileupEvent;
goodEvent_U_det_AND_pileupEvent = goodEvent_U_det && pileupEvent;
goodEvent_D_det_AND_pileupEvent = goodEvent_D_det && pileupEvent;
goodEvent_L_det_AND_pileupEvent = goodEvent_L_det && pileupEvent;
goodEvent_R_det_AND_pileupEvent = goodEvent_R_det && pileupEvent;
// if (bool_debugingRequired && muonTriggered_det) {
// std::cout<<"DEBUG: goodEvent_det: eventID="<<eventID<<std::endl;
// if (goodEvent_det) std::cout<<" ___DETECTED___"<<std::endl;
// MyPrintTree();
// MyPrintConditions();
// }
// if (pileupEvent) {
// std::cout<<"\n NEW: pileupEvent: eventID="<<eventID<<", kEntry="<<kEntry<<", posEntry="<<posEntry<< std::endl;
// std::cout<<"det_time10 = "<<det_time10<<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();
// Clone some channels into different one, if requested by user
// (This is usefull when e.g. user splits a signal from a veto
// and uses it in two different ways - e.g. once for vetoing
// muons, and second (with a different threshold) for validating
// a positron candidate. This is initiated by the
// keyword "CLONECHANNEL" in the *.v1190 file
if (bool_clonedChannelsMultimap_NotEmpty) {
// std::cout<<"det_n="<<det_n<<std::endl;
Int_t det_n_OLD=det_n;
for (Int_t i=0; i<det_n_OLD; i++) {
// std::cout<<" det_ID["<<i<<"]="<<det_ID[i]<<" edep="<<det_edep[i]<<std::endl;
clonedChannelsMultimapType::const_iterator it = clonedChannelsMultimap.find(det_ID[i]);
// std::cout<<" clonedChannelsMultimap[i]="<<clonedChannelsMultimap[i]<<std::endl;
if (it!=clonedChannelsMultimap.end()) {
int chNumTMP = it->first;
// std::cout<<"DEBUG: chNumTMP="<<chNumTMP<<" ";
std::pair<clonedChannelsMultimapType::const_iterator,clonedChannelsMultimapType::const_iterator> ret = clonedChannelsMultimap.equal_range(chNumTMP);
for (clonedChannelsMultimapType::const_iterator ittt=ret.first; ittt!=ret.second; ++ittt) {
// std::cout << " ittt->second=" << ittt->second;
int chNumNewTMP = ittt->second;
det_ID[det_n] = chNumNewTMP;
det_edep[det_n] = det_edep[i];
det_edep_el[det_n] = det_edep_el[i];
det_edep_pos[det_n] = det_edep_pos[i];
det_edep_gam[det_n] = det_edep_gam[i];
det_edep_mup[det_n] = det_edep_mup[i];
det_nsteps[det_n] = det_nsteps[i];
det_length[det_n] = det_length[i];
det_time_start[det_n] = det_time_start[i];
det_time_end[det_n] = det_time_end[i];
det_x[det_n] = det_x[i];
det_y[det_n] = det_y[i];
det_z[det_n] = det_z[i];
det_kine[det_n] = det_kine[i];
det_VrtxKine[det_n] = det_VrtxKine[i];
det_VrtxX[det_n] = det_VrtxX[i];
det_VrtxY[det_n] = det_VrtxY[i];
det_VrtxZ[det_n] = det_VrtxZ[i];
det_VrtxVolID[det_n] = det_VrtxVolID[i];
det_VrtxProcID[det_n] = det_VrtxProcID[i];
det_VrtxTrackID[det_n] = det_VrtxTrackID[i];
det_VrtxParticleID[det_n]= det_VrtxParticleID[i];
det_VvvKine[det_n] = det_VvvKine[i];
det_VvvX[det_n] = det_VvvX[i];
det_VvvY[det_n] = det_VvvY[i];
det_VvvZ[det_n] = det_VvvZ[i];
det_VvvVolID[det_n] = det_VvvVolID[i];
det_VvvProcID[det_n] = det_VvvProcID[i];
det_VvvTrackID[det_n] = det_VvvTrackID[i];
det_VvvParticleID[det_n] = det_VvvParticleID[i];
det_n++;
}
}
}
}
// 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()) {
// 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;
Double_t t1,t2;
if (bool_muDecayTimeTransformation) {
// THIS OPTION WAS SUPOSED TO ALLOW THE USER TO SOMEHOW (IN A TRICKY AND NOT 100% RELIABLE WAY)
// SHIFT THE POSITRON COUNTS TO A RESTRICTED TIME WINDOW, AS IF THE MUON DECAYED IN SOME RESTRICTED
// TIME INTERVAL. THIS, HOWEVER, DOES NOT WORK, IN THE WAY IT IS IMPLEMENTED HERE, BECAUSE THE
// MUON SPIN ROTATION AT REST IS IGNORED HERE - IT WOULD NEED SOME FURTHER WORK TO DO, AND
// BECOMES DIFFICULT TO INTERPRET LATER ==> WORK ON THIS WAS STOPPED.
Double_t diff = muDecayTime_t_max - muDecayTime_t_min;
Double_t ttt1 = det_time_start[i];
if ((ttt1 > muDecayTime_Transformation_min) && (ttt1 < muDecayTime_Transformation_max)) {
if (muDecayTime < muDecayTime_t_min) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)+1) * diff + det_time_start[i];
}
else if (muDecayTime > muDecayTime_t_max) {
ttt1 = ( Long64_t((muDecayTime_t_min-muDecayTime)/diff)) * diff + det_time_start[i];
}
}
t1 = (globTime+ttt1 )/tdcresolution;
t2 = (globTime+ttt1+dPhaseShift)/tdcresolution;
// std::cout<<"DEBUG: det_time_start[i]="<<det_time_start[i]<<" ttt1="<<ttt1<<" t1="<<t1<<" t2="<<t2<<std::endl;
}
else {
t1 = (globTime+det_time_start[i] )/tdcresolution;
t2 = (globTime+det_time_start[i]+dPhaseShift)/tdcresolution;
}
Long64_t timeBin = Long64_t(t1);
Long64_t timeBin2 = Long64_t(t2);
(*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 (bool_debugingRequired) {
if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\": pCounterMap.size()="<<pCounterMap.size()<<std::endl;}
}
if (pCounterMap.empty()) return false;
Bool_t positronHitFound = false;
Bool_t goodPositronFound = 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);
Int_t positronQuality = (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 (positronQuality>0) {
if (positronHitFound) {
if (bool_debugingRequired) {
if (debugEventMap[eventID]>2) {std::cout<<"DEBUGEVENT:"<<eventID<<"\"PositronCounterHit\": Coincidence with other positron candidate in other counter."<<std::endl;}
}
doubleHit = true;
return false;
} // There were two hits in two different positron counters
positronHitFound = true;
if (positronQuality>1) goodPositronFound = true;
}
}
}
return goodPositronFound;
}
//================================================================
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);
}
Reference in New Issue View Git Blame Copy Permalink