20 #include "FairRootManager.h"
21 #include "FairRunAna.h"
22 #include "FairRuntimeDb.h"
27 #include "TLorentzVector.h"
29 #include "TClonesArray.h"
30 #include "TDatabasePDG.h"
31 #include "TParticlePDG.h"
64 FairTask(
"Panda Analysis Task")
71 FairTask(
"Panda Analysis Task")
95 FairRootManager* ioman = FairRootManager::Instance();
97 cout <<
"-E- PndSimpleAnalysis::Init: "
98 <<
"RootManager not instantiated!" << endl;
112 fpInit.SetXYZT(0.,0.,0.,0.);
115 fChargedArray = (TClonesArray*) ioman->GetObject(
"PidChargedCand");
116 fNeutralArray = (TClonesArray*) ioman->GetObject(
"PidNeutralCand");
121 fMcArray = (TClonesArray*) ioman->GetObject(
"PndMcTracks");
125 cout <<
"-W- PndSimpleAnalysis::Init: "
126 <<
"None of PidChargedCand, PndNeutralCand, PidMcTracks, PndPidCandidates available!" << endl;
131 cout <<
"-E- PndSimpleAnalysis::Init: "
132 <<
"Error reading config file "<<
fCfgFileName<<
"." << endl;
137 cout <<
"-I- PndSimpleAnalysis: Intialization successfull" << endl;
226 int pcodes[10] = {-11,11,-13,13,211,-211,321,-321,2212,-2212};
228 for (i=0; i<40; i++) {
232 ldef->
fCharge = (float)(1-(i%2*2));
333 std::vector<std::string> shortcut;
337 shortcut.push_back(
"px");
338 shortcut.push_back(
"py");
339 shortcut.push_back(
"pz");
340 shortcut.push_back(
"en");
345 shortcut.push_back(
"pxcm");
346 shortcut.push_back(
"pycm");
347 shortcut.push_back(
"pzcm");
348 shortcut.push_back(
"encm");
353 shortcut.push_back(
"phi");
354 shortcut.push_back(
"tht");
355 shortcut.push_back(
"cth");
356 shortcut.push_back(
"p");
361 shortcut.push_back(
"phicm");
362 shortcut.push_back(
"thtcm");
363 shortcut.push_back(
"cthcm");
364 shortcut.push_back(
"pcm");
369 shortcut.push_back(
"vx");
370 shortcut.push_back(
"vy");
371 shortcut.push_back(
"vz");
376 shortcut.push_back(
"sel");
377 shortcut.push_back(
"des");
378 shortcut.push_back(
"dem");
379 shortcut.push_back(
"det");
380 shortcut.push_back(
"tcb");
381 shortcut.push_back(
"tcd");
382 shortcut.push_back(
"tcr");
383 shortcut.push_back(
"tm2");
388 shortcut.push_back(
"elh");
389 shortcut.push_back(
"mulh");
390 shortcut.push_back(
"pilh");
391 shortcut.push_back(
"klh");
392 shortcut.push_back(
"plh");
402 FairRun*
run = FairRun::Instance();
403 if ( ! run ) { Fatal(
"SetParContainers",
"No analysis run"); }
417 TVector3 pInitBoost(0,0,0);
418 if (
fpInit.Z()!=0.0) { pInitBoost=-
fpInit.BoostVector(); }
429 bool dumpAList=
false;
433 std::vector<TLorentzVector> lvcache;
438 if (!cur->
fIsUsed) {
continue; }
497 if (cur->
fHisto.size()>0) {
524 for (j=0; j<(int)curdump->
fNtpFNames.size(); j++) {
530 bool cachefilled=(lvcache.size()>0);
543 theArF[off+k]=cur->
fList[k]->Px();
546 theArF[off+k]=cur->
fList[k]->Py();
549 theArF[off+k]=cur->
fList[k]->Pz();
552 theArF[off+k]=cur->
fList[k]->E();
556 theArF[off+k]=cur->
fList[k]->M();
559 theArF[off+k]=cur->
fList[k]->Charge();
563 theArF[off+k]=cur->
fList[k]->Pos().X();
566 theArF[off+k]=cur->
fList[k]->Pos().Y();
569 theArF[off+k]=cur->
fList[k]->Pos().Z();
573 theArF[off+k]=cur->
fList[k]->GetVect().Phi();
576 theArF[off+k]=cur->
fList[k]->GetVect().Theta();
579 theArF[off+k]=cur->
fList[k]->GetVect().CosTheta();
582 theArF[off+k]=cur->
fList[k]->GetVect().Mag();
585 theArF[off+k]=(
fpInit - (cur->
fList[k]->P4()) ).M();
590 TLorentzVector lv=cur->
fList[k]->P4();
591 lv.
Boost(pInitBoost);
592 lvcache.push_back(lv);
594 theArF[off+k]=lvcache[k].Px();
599 TLorentzVector lv=cur->
fList[k]->P4();
600 lv.
Boost(pInitBoost);
601 lvcache.push_back(lv);
603 theArF[off+k]=lvcache[k].Py();
608 TLorentzVector lv=cur->
fList[k]->P4();
609 lv.
Boost(pInitBoost);
610 lvcache.push_back(lv);
612 theArF[off+k]=lvcache[k].Pz();
617 TLorentzVector lv=cur->
fList[k]->P4();
618 lv.
Boost(pInitBoost);
619 lvcache.push_back(lv);
621 theArF[off+k]=lvcache[k].E();
626 TLorentzVector lv=cur->
fList[k]->P4();
627 lv.
Boost(pInitBoost);
628 lvcache.push_back(lv);
630 theArF[off+k]=lvcache[k].Vect().Phi();
635 TLorentzVector lv=cur->
fList[k]->P4();
636 lv.
Boost(pInitBoost);
637 lvcache.push_back(lv);
639 theArF[off+k]=lvcache[k].Vect().Theta();
644 TLorentzVector lv=cur->
fList[k]->P4();
645 lv.
Boost(pInitBoost);
646 lvcache.push_back(lv);
648 theArF[off+k]=lvcache[k].Vect().CosTheta();
653 TLorentzVector lv=cur->
fList[k]->P4();
654 lv.
Boost(pInitBoost);
655 lvcache.push_back(lv);
657 theArF[off+k]=lvcache[k].Vect().Mag();
661 theArF[off+k]=cur->
fList[k]->Daughter(0)->M();
664 theArF[off+k]=cur->
fList[k]->Daughter(1)->M();
667 theArF[off+k]=cur->
fList[k]->Daughter(2)->M();
670 theArF[off+k]=cur->
fList[k]->Daughter(3)->M();
673 theArF[off+k]=cur->
fList[k]->Daughter(4)->M();
687 if (mic) { theArF[off+k]=mic->
GetTofM2(); }
697 if (mic) { theArF[off+k]=cur->
fList[k]->GetPidInfo(0); }
700 if (mic) { theArF[off+k]=cur->
fList[k]->GetPidInfo(1); }
703 if (mic) { theArF[off+k]=cur->
fList[k]->GetPidInfo(2); }
706 if (mic) { theArF[off+k]=cur->
fList[k]->GetPidInfo(3); }
709 if (mic) { theArF[off+k]=cur->
fList[k]->GetPidInfo(4); }
735 for (j=0; j<(int)curdump->
fNtpINames.size(); j++) {
761 if (0==dauC) {
break; }
774 theArI[k+off]=dauindex;
779 Error(
"PndSimpleAnalysis::Exec",
"case 510 requested and MC truth handling changed! Please fix it.");
808 if (i>=30) { pidl=4; }
809 else if (i>=20) { pidl=3; }
810 else if (i>=10) { pidl=2; }
830 theArI[k+off]=cur->
fList[k]->PdgCode();
843 if (dumpAList) {
ntp->Fill(); }
856 if (nd==0) {
return; }
861 if (level==0) { cout <<endl; }
875 if (n==
"Charged" || n==
"Neutral") {
return true; }
877 if (n==
"ElectronVeryLoose" || n==
"MuonVeryLoose" || n==
"PionVeryLoose"
878 || n==
"KaonVeryLoose" || n==
"ProtonVeryLoose") {
return true; }
880 if (n==
"ElectronLoose" || n==
"MuonLoose" || n==
"PionLoose"
881 || n==
"KaonLoose" || n==
"ProtonLoose") {
return true; }
883 if (n==
"ElectronTight" || n==
"MuonTight" || n==
"PionTight"
884 || n==
"KaonTight" || n==
"ProtonTight") {
return true; }
886 if (n==
"ElectronVeryTight" || n==
"MuonVeryTight" || n==
"PionVeryTight"
887 || n==
"KaonVeryTight" || n==
"ProtonVeryTight") {
return true; }
888 else {
return false; }
966 for (i=0; i<
fMcArray->GetEntriesFast(); i++) {
971 for (i=0; i<40; i++) {
974 int pidhint=(i%10)/2;
1031 cout<<
"setupAnalysis"<<endl;
1036 unsigned int i=0,j=0;
1038 std::vector<int> daupdgs;
1039 unsigned int daucnt=0;
1040 unsigned int daulistcnt=0;
1043 unsigned int linecnt=0;
1048 bool defineOpen=
false;
1049 bool decaySet=
false;
1058 while (!cfgFile.eof()) {
1059 cfgFile.getline(line,200);
1065 char* token = tokenizer.
GetFirst(line,
" \t");
1068 tokenVec.push_back(token);
1069 token=tokenizer.
GetNext(
" \t");
1072 if (tokenVec.size()==0) {
continue; }
1079 if (tokenVec[0]==
"DefineList") {
1106 if (tokenVec[0]==
"decayMode") {
1115 for (i=3; i<tokenVec.size(); i++) {
1119 daucnt=daupdgs.size();
1132 if (tokenVec[0]==
"daughterList") {
1140 if (daulistcnt>daucnt) {
return ErrorMessage(302,linecnt); }
1143 if (tokenVec[1]==
"McTruth") {
return ErrorMessage(305,linecnt); }
1145 std::string dName=tokenVec[1];
1146 bool isgeneric=
false;
1154 if (dName==
"Charged")
1155 switch (abs(daupdgs[daulistcnt])) {
1157 dName=
"ElectronVeryLoose";
1160 dName=
"MuonVeryLoose";
1163 dName=
"PionVeryLoose";
1166 dName=
"KaonVeryLoose";
1169 dName=
"ProtonVeryLoose";
1173 float charge=TDatabasePDG::Instance()->GetParticle(daupdgs[daulistcnt])->Charge();
1174 if (charge<-0.1) { dName+=
"M"; }
1175 else if (charge>0.1) { dName+=
"P"; }
1187 int decpdg=daupdgs[daulistcnt];
1192 if (decpdg!=listpdg && decpdg!=alistpdg) {
return ErrorMessage(304,linecnt); }
1197 if (decpdg==listpdg) {
1198 currentList->
fDauIdx.push_back(idx);
1199 }
else if (decpdg==alistpdg) {
1213 if (tokenVec[0]==
"End" || tokenVec[0]==
"EndDefine") {
1217 if (!decaySet || daucnt<2 || daulistcnt!=daucnt) {
return ErrorMessage(400,linecnt); }
1228 std::vector<int> finstate;
1229 std::vector<int> ccstate;
1233 for (i=0; i<daulistcnt; i++) {
1237 finstate.push_back(idx);
1238 ccstate.push_back(ccidx);
1246 std::sort(finstate.begin(),finstate.end());
1247 std::sort(ccstate.begin(),ccstate.end());
1251 cout <<currentList->
fName<<
" : list indices of list and cc : ";
1252 for (i=0; i<finstate.size(); i++) {
1253 cout <<
"( "<<finstate[
i]<<
" , "<<ccstate[
i]<<
" ) ";
1254 if (finstate[i]!=ccstate[i]) { hasAnti=
true; }
1262 std::string aName=currentList->
fName+
"_A";
1272 for (i=0; i<ccstate.size(); i++) {
1276 currentAntiList->
fDauIdx.push_back(ccidx);
1287 currentAntiList->
fCharge=charge;
1289 currentAntiList->
fIsUsed=
true;
1308 if (tokenVec[0]==
"pbarMom") {
1311 double p = atof(tokenVec[1].c_str());
1312 double mp = 0.938272;
1313 double E =
sqrt(mp*mp+p*p)+
mp;
1315 fpInit.SetXYZT(0.0,0.0,p,E);
1318 if (tokenVec[0]==
"Ecms") {
1321 double M = atof(tokenVec[1].c_str());
1322 double mp = 0.938272;
1324 double X = (M*M-2*mp*
mp)/(2*mp);
1325 double p =
sqrt(X*X-mp*mp);
1326 double E =
sqrt(mp*mp+p*p)+
mp;
1328 fpInit.SetXYZT(0.0,0.0,p,E);
1350 if (tokenVec[0]==
"dumpList") {
1352 if (defineOpen && tokenVec.size()<2) {
return ErrorMessage(501,linecnt); }
1355 if (!defineOpen && tokenVec.size()<3) {
return ErrorMessage(502,linecnt); }
1359 unsigned int startcols=2;
1361 std::string dName=tokenVec[1];
1366 if (dName==
"Charged") { dName=
"PionVeryLoose"; }
1386 currentList->
fColName=tokenVec[startcols-1];
1388 if (tokenVec.size()==startcols) {
1389 tokenVec.push_back(
"p4");
1390 tokenVec.push_back(
"ch");
1391 tokenVec.push_back(
"m");
1394 for (i=startcols; i<tokenVec.size(); i++) {
1395 std::string
col=tokenVec[
i];
1398 cout <<
"-W- PndSimpleAnalysis::SetupAnalysis: In '"<<
fCfgFileName
1399 <<
"', line "<<linecnt<<
": ";
1400 cout <<
"Unknown column name '"<<col<<
"'."<<endl;
1408 if (key==520 || key==521) {
continue; }
1411 if (key>=150 && key<=154) {
continue; }
1414 if (key>=501 && key<=505) {
continue; }
1419 for (j=0; j<keylist.size(); j++) { tokenVec.push_back(keylist[j]); }
1421 }
else if (key<500) { currentList->
fNtpFNames.push_back(col); }
1422 else { currentList->
fNtpINames.push_back(col); }
1426 if (dName==
"McTruth") {
1436 if (tokenVec[0]==
"histogram") {
1444 mass=TDatabasePDG::Instance()->GetParticle(currentList->
fPdgCode)->Mass();
1447 if (tokenVec.size()==2) {
1448 width=atof(tokenVec[1].c_str());
1451 double min=mass - width;
1452 double max=mass + width;
1454 if (tokenVec.size()==3) {
1455 double num1=atof(tokenVec[1].c_str());
1456 double num2=atof(tokenVec[2].c_str());
1457 if (num1<num2) {min=num1; max=num2;}
1458 else {min=num2; max=num1;}
1461 std::string hname=currentList->
fName+
"_M";
1462 std::string htitle=currentList->
fName+
" mass";
1464 TH1F*
h=
new TH1F(hname.c_str(),htitle.c_str(),100,
min,
max);
1465 currentList->
fHisto.push_back(h);
1471 if (tokenVec[0]==
"select") {
1474 if (tokenVec[1]==
"Mass") {
1476 if (currentList->
fPdgCode==0 && tokenVec.size()<4) {
1480 double mean=TDatabasePDG::Instance()->GetParticle(currentList->
fPdgCode)->Mass();
1484 if (tokenVec.size()==3) {
1485 width=atof(tokenVec[2].c_str());
1488 if (tokenVec.size()>3) {
1489 mean=atof(tokenVec[2].c_str());
1490 width=atof(tokenVec[3].c_str());
1504 std::map<std::string, int> colMap;
1510 if (0==
ntp) {
ntp=
new TTree(
"ntp",
"PndSimpleAnalysis NTuple"); }
1511 std::string pre=currentList->
fColName;
1513 int nd=currentList->
GetNDau();
1516 if (nd>0) { currentList->
fNtpINames.push_back(
"d1"); }
1517 if (nd>1) { currentList->
fNtpINames.push_back(
"d2"); }
1518 if (nd>2) { currentList->
fNtpINames.push_back(
"d3"); }
1519 if (nd>3) { currentList->
fNtpINames.push_back(
"d4"); }
1520 if (nd>4) { currentList->
fNtpINames.push_back(
"d5"); }
1522 for (i=0; i<currentList->
fNtpFNames.size(); i++) {
1526 if (nd>0) { currentList->
fNtpFNames.push_back(
"d1m"); }
1527 if (nd>1) { currentList->
fNtpFNames.push_back(
"d2m"); }
1528 if (nd>2) { currentList->
fNtpFNames.push_back(
"d3m"); }
1529 if (nd>3) { currentList->
fNtpFNames.push_back(
"d4m"); }
1530 if (nd>4) { currentList->
fNtpFNames.push_back(
"d5m"); }
1536 std::string brname=
"n"+pre;
1537 if (colMap.find(brname)!=colMap.end()) {
return ErrorMessage(504,0,brname); }
1541 ntp->Branch(brname.c_str(),&(currentList->
fNEntries),(brname+
"/I").c_str());
1544 for (i=0; i<currentList->
fNtpFNames.size(); i++) {
1547 if (colMap.find(brname)!=colMap.end()) {
return ErrorMessage(504,0,brname); }
1553 ntp->Branch(brname.c_str(),currentList->
fNtpFArrays[
i],(brname+
"[n"+pre+
"]/F").c_str());
1557 for (i=0; i<currentList->
fNtpINames.size(); i++) {
1560 if (colMap.find(brname)!=colMap.end()) {
return ErrorMessage(504,0,brname); }
1565 ntp->Branch(brname.c_str(), currentList->
fNtpIArrays[
i], (brname+
"[n"+pre+
"]/I").c_str());
1607 if (TDatabasePDG::Instance()->GetParticle(name.c_str())) {
1608 code=TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
1620 int code=0, pdgcode=0;
1621 if (TDatabasePDG::Instance()->GetParticle(name.c_str())) {
1622 pdgcode=TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
1623 if (TDatabasePDG::Instance()->GetParticle(pdgcode)->AntiParticle()) {
1640 if (TDatabasePDG::Instance()->GetParticle(pdgcode)) {
1641 if (TDatabasePDG::Instance()->GetParticle(pdgcode)->AntiParticle()) {
1657 cout <<
"-E- PndSimpleAnalysis::SetupAnalysis: In '"<<
fCfgFileName <<
"', line "<<line<<
": ";
1661 cout <<
"Missing 'End' statement!";
1665 cout <<
"List with name '"<< arg<<
"' already defined";
1669 cout <<
"Missing 'DefineList' statement!";
1673 cout <<
"DecayMode already defined!";
1677 cout <<
"More than 5 daughters not allowed!";
1681 cout <<
"Missing 'decayMode' statement!";
1685 cout <<
"Too many daughterLists defined!";
1689 cout <<
"Unknown List '"<<arg<<
"'!";
1693 cout <<
"Particle type mismatch between daughterList and decayMode!";
1697 cout <<
"McTruth cannot be used as daughterList!";
1701 cout <<
"Wrong number of daughters or daughterLists!";
1705 cout <<
"Mass of particle '"<<arg <<
"' unknown; must specify selector mean value!";
1709 cout <<
"Must specify name of tuple column!";
1713 cout <<
"Must specify name of list and tuple column!";
1717 cout <<
"Warning: Initial 4 vector already set before!";
1721 cout <<
"Branch '"<<arg<<
"' already exists in NTuple!";
1725 cout <<
"Invalid number for parameter 'MaxEntries'.";
1729 cout <<
"Unexpected error!";
int GetPdgCode(std::string name)
void Add(const RhoCandidate *c)
RhoMinusParticleSelector * minusSel
Double_t GetProtonPidProb(PndPidProbability *flux=NULL) const
void PrintTree(RhoCandidate *tc, int level=0)
int GetAntiPdgCode(std::string name)
Double_t GetKaonPidProb(PndPidProbability *flux=NULL) const
std::vector< TH1F * > fHisto
std::vector< std::string > fNtpINames
friend F32vec4 sqrt(const F32vec4 &a)
Float_t GetSttMeanDEDX() const
std::vector< std::string > fGenericListNames
void Boost(const TVector3 &)
TClonesArray * fNeutralProbability
RhoSimpleElectronSelector * eSel
TClonesArray * fChargedArray
virtual Bool_t Accept(RhoCandidate *)=0
RhoCandidate * Daughter(Int_t n)
RhoSimplePionSelector * piSel
cout<< "POINTs for new FwEndCap == "<< tsim-> GetEntriesFast()
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
TClonesArray * fMicroArray
Float_t GetRichThetaC() const
bool IsGenericListName(std::string n)
std::vector< float * > fNtpFArrays
Float_t GetDrcThetaC() const
Float_t GetDiscThetaC() const
int uid(int lev, int lrun, int lmode)
void Combine(RhoCandList &l1, RhoCandList &l2)
void SetPidInfo(double *pidinfo=0)
Double_t GetMuonPidProb(PndPidProbability *flux=NULL) const
TClonesArray * fChargedProbability
std::map< int, std::vector< std::string > > fColShortKeyMap
char * GetFirst(char *lpsz, const char *lpcszDelimiters)
void Select(RhoParticleSelectorBase *pidmgr)
Float_t GetMvdDEDX() const
std::vector< int > fDauIdx
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
static RhoFactory * Instance()
void SetConfigFile(std::string filename="analysis.cfg")
bool ErrorMessage(int mid, int line=0, std::string arg="")
std::map< std::string, int > fColKeyMap
RhoSimpleMuonSelector * muSel
virtual void SetCriterion(const char *crit)
std::vector< int * > fNtpIArrays
RhoPlusParticleSelector * plusSel
RhoSimpleKaonSelector * kSel
RhoSimpleProtonSelector * pSel
Bool_t IsCloneOf(const RhoCandidate &, Bool_t checkType=kFALSE) const
virtual InitStatus Init()
char * GetNext(const char *lpcszDelimiters)
std::vector< std::string > ArgVector
virtual void Exec(Option_t *opt)
std::vector< std::string > fNtpFNames
virtual void SetParContainers()
std::vector< RhoParticleSelectorBase * > fSelector
Double_t GetElectronPidProb(PndPidProbability *flux=NULL) const
Double_t GetPionPidProb(PndPidProbability *flux=NULL) const
double GetPidInfo(int hypo)
TClonesArray * fNeutralArray
std::map< std::string, int > fListMap
std::vector< PndListDefiner * > fListDefiners