From 902aa95c7fd8a11de13d5b38b2dd46bea705d47d Mon Sep 17 00:00:00 2001 From: gconesab Date: Tue, 17 Nov 2009 12:11:56 +0000 Subject: [PATCH] Changes in frame required to work without AODs output in case of usage of AliAnaCalorimeterQA - AliAnaCalorimeterQA: More independent of the reader, no gets ESDs or AODs directly from input event. - AliMCAnalysisUtils: New check on the overlap of photons in the clusters. Input for CheckOrigin methods is now either a single label (more contributing to the cluster) or the list of labels stored in calo clusters. --- PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx | 20 +- PWG4/PartCorrBase/AliCaloTrackReader.cxx | 20 +- PWG4/PartCorrBase/AliMCAnalysisUtils.cxx | 435 ++++++++- PWG4/PartCorrBase/AliMCAnalysisUtils.h | 13 +- PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx | 1051 +++++++++++++-------- PWG4/PartCorrDep/AliAnaCalorimeterQA.h | 21 +- PWG4/macros/AddTaskCalorimeterQA.C | 34 +- PWG4/macros/ConfigAnalysisCalorimeterQA.C | 43 +- 8 files changed, 1097 insertions(+), 540 deletions(-) diff --git a/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx b/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx index f2a245cd003..38384e7046a 100755 --- a/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx +++ b/PWG4/PartCorrBase/AliAnaPartCorrMaker.cxx @@ -228,21 +228,25 @@ void AliAnaPartCorrMaker::ProcessEvent(const Int_t iEntry, const char * currentF abort(); } - if(fAnaDebug >= 0 ){ + if(fAnaDebug >= 0 ){ printf("*** Event %d *** \n",iEntry); - if(fAnaDebug > 1 ) - printf("AliAnaPartCorrMaker::Current File Name : %s\n", currentFileName); + if(fAnaDebug > 1 ) { + printf("AliAnaPartCorrMaker::ProcessEvent() - Current File Name : %s\n", currentFileName); + //printf("fAODBranchList %p, entries %d\n",fAODBranchList,fAODBranchList->GetEntries()); + } } - //Each event needs an empty branch - for(Int_t iaod = 0; iaod < fAODBranchList->GetEntries(); iaod++) - fAODBranchList->At(iaod)->Clear(); - + //Each event needs an empty branch + Int_t nAODBranches = fAODBranchList->GetEntries(); + for(Int_t iaod = 0; iaod < nAODBranches; iaod++) + fAODBranchList->At(iaod)->Clear(); + //Tell the reader to fill the data in the 3 detector lists Bool_t ok = fReader->FillInputEvent(iEntry, currentFileName); if(!ok){ printf("*** Skip event *** %d \n",iEntry); return ; } + //printf(">>>>>>>>>> BEFORE >>>>>>>>>>>\n"); //gObjectTable->Print(); //Loop on analysis algorithms @@ -250,7 +254,7 @@ void AliAnaPartCorrMaker::ProcessEvent(const Int_t iEntry, const char * currentF Int_t nana = fAnalysisContainer->GetEntries() ; for(Int_t iana = 0; iana < nana; iana++){ AliAnaPartCorrBaseClass * ana = ((AliAnaPartCorrBaseClass *) fAnalysisContainer->At(iana)) ; - ana->ConnectInputOutputAODBranches(); //Sets branches for each analysis + if(nAODBranches) ana->ConnectInputOutputAODBranches(); //Sets branches for each analysis, if there is any branch. //Make analysis, create aods in aod branch or AODCaloClusters if(fMakeAOD) ana->MakeAnalysisFillAOD() ; //Make further analysis with aod branch and fill histograms diff --git a/PWG4/PartCorrBase/AliCaloTrackReader.cxx b/PWG4/PartCorrBase/AliCaloTrackReader.cxx index 67e894e72ea..5a1b273d278 100755 --- a/PWG4/PartCorrBase/AliCaloTrackReader.cxx +++ b/PWG4/PartCorrBase/AliCaloTrackReader.cxx @@ -344,9 +344,10 @@ void AliCaloTrackReader::InitParameters() fEMCALPtMin = 0.2 ; fPHOSPtMin = 0.2 ; - fFillEMCAL = kTRUE; - fFillPHOS = kTRUE; - fFillCTS = kTRUE; + //Do not filter the detectors input by default. + fFillEMCAL = kFALSE; + fFillPHOS = kFALSE; + fFillCTS = kFALSE; fFillEMCALCells = kFALSE; fFillPHOSCells = kFALSE; @@ -397,15 +398,15 @@ void AliCaloTrackReader::Print(const Option_t * opt) const //___________________________________________________ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * currentFileName) { //Fill the event counter and input lists that are needed, called by the analysis maker. - + fEventNumber = iEntry; fCurrentFileName = TString(currentFileName); - if((fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){ + if(fOutputEvent && (fDataType != kAOD) && ((fOutputEvent->GetCaloClusters())->GetEntriesFast()!=0 ||(fOutputEvent->GetTracks())->GetEntriesFast()!=0)){ printf("AliCaloTrackReader::AODCaloClusters or AODTracks already filled by the filter, do not use the ESD reader, use the AOD reader, STOP\n"); abort(); } - + //In case of analysis of events with jets, skip those with jet pt > 5 pt hard if(fComparePtHardAndJetPt && GetStack()) { if(!ComparePtHardAndJetPt()) return kFALSE ; @@ -421,7 +422,8 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * curre printf("AliCaloTrackReader::FillInputEvent() - Skip events from event %d, no more events in second AOD file \n", iEntry); return kFALSE; } - + printf("Reader 4 \n"); + //Get the Event Int_t nbytes = fSecondInputAODTree->GetEvent(iEntry+fSecondInputFirstEvent); if ( nbytes == 0 ) {//If nothing in AOD @@ -430,7 +432,7 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * curre } } - //if(iEntry > 10) return kFALSE; + if(fFillCTS) FillInputCTS(); if(fFillEMCAL) FillInputEMCAL(); if(fFillPHOS) FillInputPHOS(); @@ -449,7 +451,7 @@ void AliCaloTrackReader::ResetLists() { if(fAODPHOS) fAODPHOS -> Clear(); if(fEMCALCells) fEMCALCells -> Clear(); if(fPHOSCells) fPHOSCells -> Clear(); - if(fCleanOutputStdAOD){ + if(fCleanOutputStdAOD && fOutputEvent ){ //Only keep copied tracks and clusters if requested fOutputEvent->GetTracks() ->Clear(); fOutputEvent->GetCaloClusters()->Clear(); diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx index 52f755286f4..24415c97069 100755 --- a/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.cxx @@ -83,24 +83,41 @@ AliMCAnalysisUtils::~AliMCAnalysisUtils() } +//_________________________________________________________________________ +Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, const Int_t nlabels, AliCaloTrackReader* reader, const Int_t input = 0) { + //Play with the montecarlo particles if available + Int_t tag = 0; + + //Select where the information is, ESD-galice stack or AOD mcparticles branch + if(reader->ReadStack()){ + tag = CheckOriginInStack(label, nlabels, reader->GetStack()); + } + else if(reader->ReadAODMCParticles()){ + tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(input)); + } + + return tag ; +} + //_________________________________________________________________________ Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, AliCaloTrackReader* reader, const Int_t input = 0) { //Play with the montecarlo particles if available Int_t tag = 0; + Int_t labels[]={label}; //Select where the information is, ESD-galice stack or AOD mcparticles branch if(reader->ReadStack()){ - tag = CheckOriginInStack(label, reader->GetStack()); + tag = CheckOriginInStack(labels, 1,reader->GetStack()); } else if(reader->ReadAODMCParticles()){ - tag = CheckOriginInAOD(label, reader->GetAODMCParticles(input)); + tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(input)); } return tag ; } //_________________________________________________________________________ -Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) { +Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack* stack) { // Play with the MC stack if available. Tag particles depending on their origin. // Do same things as in CheckOriginInAOD but different input. @@ -115,9 +132,12 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) } Int_t tag = 0; + Int_t label=labels[0];//Most significant particle contributing to the cluster + if(label >= 0 && label < stack->GetNtrack()){ //MC particle of interest is the "mom" of the entity TParticle * mom = stack->Particle(label); + Int_t iMom = label; Int_t mPdg = TMath::Abs(mom->GetPdgCode()); Int_t mStatus = mom->GetStatusCode() ; Int_t iParent = mom->GetFirstMother() ; @@ -134,33 +154,54 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) } else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent); + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n"); + printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); + printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus); + } + //Check if "mother" of entity is converted, if not, get the first non converted mother if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){ SetTagBit(tag,kMCConversion); //Check if the mother is photon or electron with status not stable while ((pPdg == 22 || pPdg == 11) && mStatus != 1) { //Mother - mom = stack->Particle(mom->GetFirstMother()); + iMom = mom->GetFirstMother(); + mom = stack->Particle(iMom); mPdg = TMath::Abs(mom->GetPdgCode()); mStatus = mom->GetStatusCode() ; iParent = mom->GetFirstMother() ; if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent); //GrandParent - if(iParent > 0){ + if(iParent >= 0){ parent = stack->Particle(iParent); pPdg = TMath::Abs(parent->GetPdgCode()); pStatus = parent->GetStatusCode(); } }//while + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n"); + printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); + printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus); + } + }//mother and parent are electron or photon and have status 0 - else if(mStatus == 0){ - //Still a conversion but only one electron/photon generated. Just from hadrons - if(pPdg == 2112 || pPdg == 211 || - pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 ) - SetTagBit(tag,kMCConversion); - mom = stack->Particle(mom->GetFirstMother()); - mPdg = TMath::Abs(mom->GetPdgCode()); + else if((mPdg == 22 || mPdg == 11) && mStatus == 0){ + //Still a conversion but only one electron/photon generated. Just from hadrons but not decays. + if(pPdg == 2112 || pPdg == 211 || pPdg == 321 || + pPdg == 2212 || pPdg == 130 || pPdg == 13 ) { + SetTagBit(tag,kMCConversion); + iMom = mom->GetFirstMother(); + mom = stack->Particle(iMom); + mPdg = TMath::Abs(mom->GetPdgCode()); + + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n"); + printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); + } + }//hadron converted + //Comment for the next lines, we do not check the parent of the hadron for the moment. //iParent = mom->GetFirstMother() ; //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent); @@ -171,7 +212,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) // pPdg = TMath::Abs(parent->GetPdgCode()); //} } - // conversion fo electrons/photons checked + // conversion into electrons/photons checked //first check for typical charged particles if(mPdg == 13) SetTagBit(tag,kMCMuon); @@ -179,8 +220,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) else if(mPdg == 321) SetTagBit(tag,kMCKaon); else if(mPdg == 2212) SetTagBit(tag,kMCProton); //check for pi0 and eta (shouldn't happen unless their decays were turned off) - else if(mPdg == 111) SetTagBit(tag,kMCPi0); - else if(mPdg == 221) SetTagBit(tag,kMCEta); + else if(mPdg == 111) { + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if(mPdg == 221) { + SetTagBit(tag,kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster + } //Photons else if(mPdg == 22){ SetTagBit(tag,kMCPhoton); @@ -194,9 +243,17 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) SetTagBit(tag, kMCISR); //Initial state radiation } else if(pStatus == 11){//Decay - if(pPdg == 111) SetTagBit(tag,kMCPi0Decay); - else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay); - else SetTagBit(tag,kMCOtherDecay); + if(pPdg == 111) { + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if (pPdg == 221) { + SetTagBit(tag, kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster + } + else SetTagBit(tag,kMCOtherDecay); }//Decay else { printf("AliMCAnalysisUtils::ChecOrigingInAOD() - what is it? Mother mPdg %d, status %d \n Parent iParent %d, pPdg %d %s, status %d\n", @@ -212,7 +269,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) iParent = parent->GetFirstMother(); parent=stack->Particle(iParent); pStatus= parent->GetStatusCode(); - pPdg = parent->GetPdgCode(); + pPdg = TMath::Abs(parent->GetPdgCode()); }//Look for the parton if(iParent < 8 && iParent > 5) { @@ -222,8 +279,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) else SetTagBit(tag,kMCISR);//Initial state radiation }//Not decay else{//Decay - if(pPdg == 111) SetTagBit(tag,kMCPi0Decay); - else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay); + if(pPdg == 111) { + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if (pPdg == 221) { + SetTagBit(tag,kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster + } else SetTagBit(tag,kMCOtherDecay); }//Decay }//HERWIG @@ -233,8 +298,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) }//Status 1 : created by event generator else if(mStatus == 0){ // geant - if(pPdg == 111) SetTagBit(tag,kMCPi0Decay); - else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay); + if(pPdg == 111) { + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if (pPdg == 221) { + SetTagBit(tag,kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster + } else SetTagBit(tag,kMCOtherDecay); }//status 0 : geant generated @@ -253,7 +326,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) } } SetTagBit(tag,kMCElectron); - if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons"); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n"); if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay @@ -269,18 +342,21 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) //if it is not from any of the above, where is it from? if(pPdg > 10000) SetTagBit(tag,kMCUnknown); else SetTagBit(tag,kMCOtherDecay); - if(fDebug > 0) printf("STACK Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg); } }//electron check //Cluster was made by something else else { - if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg); SetTagBit(tag,kMCUnknown); } }//Good label value - else{ - if(label < 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label); - if(label >= stack->GetNtrack()) printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); + else{// Bad label + + if(label < 0 && (fDebug >= 0)) + printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label); + if(label >= stack->GetNtrack() && (fDebug >= 0)) + printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); SetTagBit(tag,kMCUnknown); }//Bad label @@ -290,7 +366,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t label, AliStack* stack) //_________________________________________________________________________ -Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcparticles) { +Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, const Int_t nlabels, TClonesArray *mcparticles) { // Play with the MCParticles in AOD if available. Tag particles depending on their origin. // Do same things as in CheckOriginInStack but different input. if(!mcparticles) { @@ -299,10 +375,13 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa } Int_t tag = 0; + Int_t label=labels[0];//Most significant particle contributing to the cluster + Int_t nprimaries = mcparticles->GetEntriesFast(); if(label >= 0 && label < nprimaries){ //Mother AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label); + Int_t iMom = label; Int_t mPdg = TMath::Abs(mom->GetPdgCode()); Int_t iParent = mom->GetMother() ; if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); @@ -316,13 +395,20 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa } else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent); + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n"); + printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); + printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); + } + //Check if mother is converted, if not, get the first non converted mother if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){ SetTagBit(tag,kMCConversion); //Check if the mother is photon or electron with status not stable while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) { //Mother - mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother()); + iMom = mom->GetMother(); + mom = (AliAODMCParticle *) mcparticles->At(iMom); mPdg = TMath::Abs(mom->GetPdgCode()); iParent = mom->GetMother() ; if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); @@ -332,15 +418,33 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa parent = (AliAODMCParticle *) mcparticles->At(iParent); pPdg = TMath::Abs(parent->GetPdgCode()); } - }//while + // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); + // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); + + }//while + + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n"); + printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); + printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); + } + }//mother and parent are electron or photon and have status 0 and parent is photon or electron - else if(!mom->IsPrimary()){ + else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){ //Still a conversion but only one electron/photon generated. Just from hadrons - if(pPdg == 2112 || pPdg == 211 || - pPdg == 321 || pPdg == 2212 || pPdg == 130 || pPdg == 13 ) - SetTagBit(tag,kMCConversion); - mom = (AliAODMCParticle *) mcparticles->At(mom->GetMother()); - mPdg = TMath::Abs(mom->GetPdgCode()); + if(pPdg == 2112 || pPdg == 211 || pPdg == 321 || + pPdg == 2212 || pPdg == 130 || pPdg == 13 ) { + SetTagBit(tag,kMCConversion); + iMom = mom->GetMother(); + mom = (AliAODMCParticle *) mcparticles->At(iMom); + mPdg = TMath::Abs(mom->GetPdgCode()); + + if(fDebug > 2 ) { + printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n"); + printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); + } + }//hadron converted + //Comment for next lines, we do not check the parent of the hadron for the moment. //iParent = mom->GetMother() ; //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); @@ -360,8 +464,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa else if(mPdg == 321) SetTagBit(tag,kMCKaon); else if(mPdg == 2212) SetTagBit(tag,kMCProton); //check for pi0 and eta (shouldn't happen unless their decays were turned off) - else if(mPdg == 111) SetTagBit(tag,kMCPi0); - else if(mPdg == 221) SetTagBit(tag,kMCEta); + else if(mPdg == 111) { + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if(mPdg == 221) { + SetTagBit(tag,kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster + } //Photons else if(mPdg == 22){ SetTagBit(tag,kMCPhoton); @@ -374,8 +486,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa SetTagBit(tag, kMCISR); //Initial state radiation } else if(parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay - if(pPdg == 111) SetTagBit(tag,kMCPi0Decay); - else if (pPdg == 221) SetTagBit(tag, kMCEtaDecay); + if(pPdg == 111){ + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if (pPdg == 221) { + SetTagBit(tag, kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster + } else SetTagBit(tag,kMCOtherDecay); }//Decay else { @@ -385,8 +505,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa } }// Pythia generated else if(!mom->IsPrimary()){ //Decays - if(pPdg == 111) SetTagBit(tag,kMCPi0Decay); - else if (pPdg == 221) SetTagBit(tag,kMCEtaDecay); + if(pPdg == 111){ + SetTagBit(tag,kMCPi0Decay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster + } + else if (pPdg == 221) { + SetTagBit(tag,kMCEtaDecay); + if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n"); + CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster + } else SetTagBit(tag,kMCOtherDecay); }//not primary : geant generated, decays else { @@ -424,28 +552,237 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t label, TClonesArray *mcpa } else { //prompt or other decay TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg); TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg); - if(fDebug > 0) printf("AOD Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(),pPdg,foo1->GetName(),mPdg); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg); if(pPdg > 10000) SetTagBit(tag,kMCUnknown); else SetTagBit(tag,kMCOtherDecay); } }//electron check //cluster was made by something else else { - if(fDebug > 0) printf("\tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg); + if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg); SetTagBit(tag,kMCUnknown); } }//Good label value - else{ - if(label < 0 && (fDebug >= 0) ) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no stack ***: label %d \n", label); - if(label >= mcparticles->GetEntriesFast() && (fDebug > 0) ) printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast()); + else{//Bad label + + if(label < 0 && (fDebug >= 0) ) + printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label); + if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) ) + printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast()); SetTagBit(tag,kMCUnknown); + }//Bad label return tag; } +//_________________________________________________________________________ +void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex, + AliStack *stack, Int_t &tag) { + //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack + + if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) { + if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n", + labels[0],stack->GetNtrack(), nlabels); + return; + } + + TParticle * meson = stack->Particle(mesonIndex); + Int_t mesonPdg = meson->GetPdgCode(); + if(mesonPdg!=111 && mesonPdg!=221){ + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg); + return; + } + + if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex); + //Check if meson decayed into 2 daughters or if both were kept. + if(meson->GetNDaughters() != 2){ + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters()); + return; + } + + //Get the daughters + Int_t iPhoton0 = meson->GetDaughter(0); + Int_t iPhoton1 = meson->GetDaughter(1); + TParticle *photon0 = stack->Particle(iPhoton0); + TParticle *photon1 = stack->Particle(iPhoton1); + + //Check if both daughters are photons + if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){ + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode()); + return; + } + + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1); + + //Check if both photons contribute to the cluster + Bool_t okPhoton0 = kFALSE; + Bool_t okPhoton1 = kFALSE; + + if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n"); + + for(Int_t i = 0; i < nlabels; i++){ + if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1); + + //If we already found both, break the loop + if(okPhoton0 && okPhoton1) break; + + Int_t index = labels[i]; + if (iPhoton0 == index) { + okPhoton0 = kTRUE; + continue; + } + else if (iPhoton1 == index) { + okPhoton1 = kTRUE; + continue; + } + + //Trace back the mother in case it was a conversion + TParticle * daught = stack->Particle(index); + Int_t tmpindex = daught->GetFirstMother(); + if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex); + while(tmpindex>=0){ + //MC particle of interest is the mother + if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex); + daught = stack->Particle(tmpindex); + if (iPhoton0 == tmpindex) { + okPhoton0 = kTRUE; + break; + } + else if (iPhoton1 == tmpindex) { + okPhoton1 = kTRUE; + break; + } + tmpindex = daught->GetFirstMother(); + }//While to check if pi0/eta daughter was one of these contributors to the cluster + + if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n"); + + }//loop on list of labels + + //If both photons contribute tag as the corresponding meson. + if(okPhoton0 && okPhoton1){ + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName()); + + if(mesonPdg == 111) SetTagBit(tag,kMCPi0); + else SetTagBit(tag,kMCEta); + } + +} + +//_________________________________________________________________________ +void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex, + TClonesArray *mcparticles, Int_t &tag) { + //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles + + if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) { + if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n", + labels[0],mcparticles->GetEntriesFast(), nlabels); + return; + } + + + AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex); + Int_t mesonPdg = meson->GetPdgCode(); + if(mesonPdg != 111 && mesonPdg != 221) { + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg); + return; + } + + if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters()); + + + //Get the daughters + if(meson->GetNDaughters() != 2){ + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters()); + return; + } + Int_t iPhoton0 = meson->GetDaughter(0); + Int_t iPhoton1 = meson->GetDaughter(1); + //if((iPhoton0 == -1) || (iPhoton1 == -1)){ + // if(fDebug > 2) + // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overalapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1); + // return; + //} + AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0); + AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1); + + //Check if both daughters are photons + if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22){ + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode()); + return; + } + + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1); + + //Check if both photons contribute to the cluster + Bool_t okPhoton0 = kFALSE; + Bool_t okPhoton1 = kFALSE; + + if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n"); + + for(Int_t i = 0; i < nlabels; i++){ + if(fDebug > 3) printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1); + + //If we already found both, break the loop + if(okPhoton0 && okPhoton1) break; + + Int_t index = labels[i]; + if (iPhoton0 == index) { + okPhoton0 = kTRUE; + continue; + } + else if (iPhoton1 == index) { + okPhoton1 = kTRUE; + continue; + } + + //Trace back the mother in case it was a conversion + AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index); + Int_t tmpindex = daught->GetMother(); + if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex); + + while(tmpindex>=0){ + + //MC particle of interest is the mother + if(fDebug > 3) printf("\t parent index %d\n",tmpindex); + daught = (AliAODMCParticle*) mcparticles->At(tmpindex); + //printf("tmpindex %d\n",tmpindex); + if (iPhoton0 == tmpindex) { + okPhoton0 = kTRUE; + break; + } + else if (iPhoton1 == tmpindex) { + okPhoton1 = kTRUE; + break; + } + tmpindex = daught->GetMother(); + }//While to check if pi0/eta daughter was one of these contributors to the cluster + + if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 ) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n"); + + }//loop on list of labels + + //If both photons contribute tag as the corresponding meson. + if(okPhoton0 && okPhoton1){ + if(fDebug > 2) + printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()); + + if(mesonPdg == 111) SetTagBit(tag,kMCPi0); + else SetTagBit(tag,kMCEta); + } + +} //_________________________________________________________________________ TList * AliMCAnalysisUtils::GetJets(AliCaloTrackReader * const reader){ diff --git a/PWG4/PartCorrBase/AliMCAnalysisUtils.h b/PWG4/PartCorrBase/AliMCAnalysisUtils.h index 5284a574f9c..980d613a233 100755 --- a/PWG4/PartCorrBase/AliMCAnalysisUtils.h +++ b/PWG4/PartCorrBase/AliMCAnalysisUtils.h @@ -41,11 +41,18 @@ public: kMCElectron, kMCEFromCFromB, kMCEFromC, kMCEFromB, kMCZDecay, kMCWDecay, kMCMuon, kMCPion, kMCPi0, kMCKaon, kMCEta, kMCProton, kMCOther, kMCUnknown}; - + + //Check only the label of the most significant particle Int_t CheckOrigin(const Int_t label, AliCaloTrackReader * reader, const Int_t input) ; - Int_t CheckOriginInStack(const Int_t label, AliStack * stack) ; - Int_t CheckOriginInAOD(const Int_t label, TClonesArray* mcparticles) ; + //Check the label of the most significant particle but do checks on the rest of the contributing labels + Int_t CheckOrigin(const Int_t *label, const Int_t nlabels, AliCaloTrackReader * reader, const Int_t input) ; + Int_t CheckOriginInStack(const Int_t *labels, const Int_t nlabels, AliStack * stack) ; + Int_t CheckOriginInAOD (const Int_t *labels, const Int_t nlabels, TClonesArray* mcparticles) ; + + void CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex, AliStack * stack, Int_t & tag); + void CheckOverlapped2GammaDecay(const Int_t *labels, const Int_t nlabels, const Int_t mesonIndex, TClonesArray* mcparticles, Int_t & tag); + TList * GetJets(AliCaloTrackReader * const reader) ; void SetTagBit(Int_t &tag, const UInt_t set) const { diff --git a/PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx b/PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx index b4467a8d3fb..4703e87e3cf 100755 --- a/PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx +++ b/PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx @@ -39,11 +39,18 @@ #include "AliCaloTrackReader.h" #include "AliStack.h" #include "AliAODCaloCells.h" +#include "AliESDCaloCells.h" #include "AliAODCaloCluster.h" #include "AliFidutialCut.h" #include "AliESDtrack.h" +#include "AliAODTrack.h" #include "AliESDCaloCluster.h" #include "AliESDEvent.h" +#include "AliAODEvent.h" +#include "AliVEventHandler.h" +#include "AliAnalysisManager.h" +#include "AliAODMCParticle.h" +#include "AliMCAnalysisUtils.h" ClassImp(AliAnaCalorimeterQA) @@ -131,10 +138,10 @@ AliAnaCalorimeterQA & AliAnaCalorimeterQA::operator = (const AliAnaCalorimeterQA if(this == &qa)return *this; ((AliAnaPartCorrBaseClass *)this)->operator=(qa); - fCalorimeter = qa.fCalorimeter; - fStyleMacro = qa.fStyleMacro; - fMakePlots = qa.fMakePlots; - + fCalorimeter = qa.fCalorimeter; + fStyleMacro = qa.fStyleMacro; + fMakePlots = qa.fMakePlots; + fhE = qa.fhE; fhPt = qa.fhPt; fhPhi = qa.fhPhi; @@ -856,460 +863,661 @@ void AliAnaCalorimeterQA::Print(const Option_t * opt) const //__________________________________________________________________ void AliAnaCalorimeterQA::MakeAnalysisFillHistograms() { - //Do analysis and fill histograms + //Fill Calorimeter QA histograms + TLorentzVector mom ; TLorentzVector mom2 ; - //Play with the MC stack if available - AliStack * stack = 0x0; - TParticle * primary = 0x0; - if(IsDataMC()) { - stack = GetMCStack() ; - if(!stack) { - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Stack not available, have you switched on the MC data?\n"); - return; + TRefArray * caloClusters = new TRefArray(); + Int_t nLabel = 0; + Int_t *labels=0x0; + Int_t nCaloClusters = 0; + Int_t nCaloCellsPerCluster = 0; + Int_t nTracksMatched = 0; + Int_t trackIndex = 0; + + //Play with the MC stack if available + //Get the MC arrays and do some checks + if(IsDataMC()){ + if(GetReader()->ReadStack()){ + + if(!GetMCStack()) { + printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n"); + abort(); + } + //Fill some pure MC histograms, only primaries. + for(Int_t i=0 ; iGetNprimary(); i++){//Only primary particles, for all MC transport put GetNtrack() + TParticle *primary = GetMCStack()->Particle(i) ; + //printf("i %d, %s: status = %d, primary? %d\n",i, primary->GetName(), primary->GetStatusCode(), primary->IsPrimary()); + if (primary->GetStatusCode() > 11) continue; //Working for PYTHIA and simple generators, check for HERWIG + primary->Momentum(mom); + MCHistograms(mom,TMath::Abs(primary->GetPdgCode())); + } //primary loop + } + else if(GetReader()->ReadAODMCParticles()){ + + if(!GetReader()->GetAODMCParticles(0)) { + printf("AliAnaPhoton::MakeAnalysisFillHistograms() - AODMCParticles not available!\n"); + abort(); + } + //Fill some pure MC histograms, only primaries. + for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles(0))->GetEntriesFast(); i++){ + AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(i) ; + //printf("i %d, %s: primary? %d physical primary? %d, flag %d\n", + // i,(TDatabasePDG::Instance()->GetParticle(aodprimary->GetPdgCode()))->GetName(), + // aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), aodprimary->GetFlag()); + if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons + //aodprimary->Momentum(mom); + mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E()); + MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode())); + } //primary loop + + } + }// is data and MC + + + //Get List with CaloClusters + + if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) { + if (fCalorimeter == "EMCAL") ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL(); + else if(fCalorimeter == "PHOS") ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS(); + else { + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()); + abort(); + } + } + else if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD) { + if (fCalorimeter == "EMCAL") ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALClusters(caloClusters);//GetAODEMCAL(); + else if(fCalorimeter == "PHOS") ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSClusters (caloClusters);//GetAODPHOS(); + else { + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()); + abort(); } } - //Get List with clusters - TObjArray * partList = new TObjArray(); - if(fCalorimeter == "EMCAL") partList = GetAODEMCAL(); - else if(fCalorimeter == "PHOS") partList = GetAODPHOS(); - else { - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Wrong calorimeter name <%s>, END\n", fCalorimeter.Data()); + if(!caloClusters) { + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - No CaloClusters available\n"); abort(); } - if(!partList || partList->GetEntriesFast() == 0) return ; - - Int_t nclusters = partList->GetEntriesFast() ; - fhNClusters->Fill(nclusters); - + nCaloClusters = caloClusters->GetEntriesFast() ; + fhNClusters->Fill(nCaloClusters); + if(GetDebug() > 0) - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nclusters); + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In %s there are %d clusters \n", fCalorimeter.Data(), nCaloClusters); //Get vertex for photon momentum calculation Double_t v[3] ; //vertex ; GetReader()->GetVertex(v); + TObject * track = 0x0; + //Loop over CaloClusters + for(Int_t iclus = 0; iclus < nCaloClusters; iclus++){ + + if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster: %d/%d, data %d \n", + iclus+1,nCaloClusters,GetReader()->GetDataType()); + + if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){ + AliESDCaloCluster* clus = (AliESDCaloCluster*) (caloClusters->At(iclus)); + //Get cluster kinematics + clus->GetMomentum(mom,v); + //MC labels + nLabel = clus->GetNLabels(); + if(clus->GetLabels()) labels = (clus->GetLabels())->GetArray(); + //Cells per cluster + nCaloCellsPerCluster = clus->GetNCells(); + //matched cluster with tracks + nTracksMatched = clus->GetNTracksMatched(); + trackIndex = clus->GetTrackMatched(); + if(trackIndex >= 0){ + track = (AliESDtrack*) ((AliESDEvent*)GetReader()->GetInputEvent())->GetTrack(trackIndex); + } + else{ + if(nTracksMatched == 1) nTracksMatched = 0; + track = 0; + } + } + else{ + AliAODCaloCluster* clus = (AliAODCaloCluster*) (caloClusters->At(iclus)); - for(Int_t iclus = 0; iclus < partList->GetEntriesFast(); iclus++){ - //printf(" cluster %d\n",iclus); - AliAODCaloCluster * calo = (AliAODCaloCluster*) (partList->At(iclus)); - //if(calo->GetNCells() <= 2) continue; - //Get cluster kinematics - calo->GetMomentum(mom,v); - Float_t e = mom.E(); - //if(e < 0.5) continue; - //printf("e %2.2f, n %f\n",e, calo->GetNCells()); - Float_t pt = mom.Pt(); - Float_t eta = mom.Eta(); - Float_t phi = mom.Phi(); - if(phi < 0) phi +=TMath::TwoPi(); - if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - cluster %d: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",iclus,e,pt,eta,phi*TMath::RadToDeg()); - fhE ->Fill(e); - fhPt ->Fill(pt); - fhPhi ->Fill(phi); - fhEta ->Fill(eta); - fhEtaPhi->Fill(eta,phi); - - //matched cluster with tracks - Int_t ntracksmatched = calo->GetNTracksMatched(); - - //Fill histograms only possible when simulation - if(IsDataMC()){ - //Play with the MC stack if available - Int_t label = calo->GetLabel(0); - - if(label < 0 || !stack) { - if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** bad label or no stack ***: label %d \n", label); - continue; - } - - if(label >= stack->GetNtrack()) { - if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); - continue ; - } - - //cout<<"LABEL "<Particle(label); + //Get cluster kinematics + clus->GetMomentum(mom,v); + //MC labels + nLabel = clus->GetNLabel(); + if(clus->GetLabels()) labels = clus->GetLabels(); + //Cells per cluster + nCaloCellsPerCluster = clus->GetNCells(); + //matched cluster with tracks + nTracksMatched = clus->GetNTracksMatched(); + if(nTracksMatched > 0) + track = (AliAODTrack*)clus->GetTrackMatched(0); + } - Int_t pdg = primary->GetPdgCode(); - Float_t vx = primary->Vx(); - Float_t vy = primary->Vy(); - //Float_t vz = primary->Vz(); - Float_t r = TMath::Sqrt(vx*vx + vy*vy); - if((pdg == 22 || TMath::Abs(pdg)==11) && primary->GetStatusCode()!=1) { - fhEMVxyz ->Fill(vx,vy);//,vz); - fhEMR ->Fill(e,r); - } - - //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta); - //printf("prim e %f, pt %f, phi %f, eta %f \n", primary->Energy(),primary->Pt() ,primary->Phi() ,primary->Eta() ); - //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vx, vy, vz, r); - - //Get final particle, no conversion products - Int_t status = primary->GetStatusCode(); - Int_t mother= primary->GetFirstMother(); - if(status == 0){ - while(mother >= 0){ - primary = GetMCStack()->Particle(mother); - status = primary->GetStatusCode(); - mother= primary->GetFirstMother(); - if(status == 1) break; - //printf("mother %d\n",mother); - } - } - - fh2E ->Fill(e, primary->Energy()); - fh2Pt ->Fill(pt, primary->Pt()); - fh2Phi ->Fill(phi, primary->Phi()); - fh2Eta ->Fill(eta, primary->Eta()); - fhDeltaE ->Fill(primary->Energy()-e); - fhDeltaPt ->Fill(primary->Pt()-pt); - fhDeltaPhi->Fill(primary->Phi()-phi); - fhDeltaEta->Fill(primary->Eta()-eta); - if(primary->Energy() > 0) fhRatioE ->Fill(e/primary->Energy()); - if(primary->Pt() > 0) fhRatioPt ->Fill(pt/primary->Pt()); - if(primary->Phi() > 0) fhRatioPhi->Fill(phi/primary->Phi()); - if(primary->Eta() > 0) fhRatioEta->Fill(eta/primary->Eta()); + //Fill histograms related to single cluster or track matching + ClusterHistograms(mom, nCaloCellsPerCluster, nTracksMatched, track, labels, nLabel); - //cout<<"Final Label "<GetPdgCode(); - Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); - - if(pdg == 22){ - //cout<<"pdg "<GetStatusCode()<GetNLabel() > 1){ - if(pdg !=111){ - while(mother >= 0){ - pi0 = GetMCStack()->Particle(mother); - pdgpi0 = pi0->GetPdgCode(); - if(pdgpi0 == 111) { - //cout<<"pi0!!!"<GetFirstMother(); - //printf("mother %d\n",mother); - } - } - else pi0 = primary; - - if(!pi0 || mother < 0 ) continue ; - //cout<<"pi0 pointer "<GetPdgCode()<Pt()<<" status "<GetStatusCode()<GetNDaughters() == 2){ - //cout<<"pi0, 2 daughters "<GetFirstDaughter(); - Int_t id2 = pi0->GetFirstDaughter()+1; - p1=GetMCStack()->Particle(id1); - p2=GetMCStack()->Particle(id2); - - //if(p1->GetFirstMother()!=p2->GetFirstMother()) cout <<"Decay photon mothers are not the same!!"<GetPdgCode()==22 && p2->GetPdgCode()==22){ - // cout<<"2 photons, labels "<< id1<<" "<GetNLabel(); ilabel++){ - Int_t iprim = calo->GetLabel(ilabel); - //cout<<"iprim "<= 0){ - tmp = GetMCStack()->Particle(mothertmp); - mothertmp= tmp->GetFirstMother(); - // cout<<"mothertmp "<GetName()<< " pt "<Pt()<Energy()<Fill(e,primary->Energy()); - fhPi0Pt ->Fill(pt,primary->Pt()); - fhPi0Eta ->Fill(eta,primary->Eta()); - fhPi0Phi ->Fill(phi,primary->Phi()); - if( ntracksmatched > 0){ - fhPi0ECharged ->Fill(e,primary->Energy()); - fhPi0PtCharged ->Fill(pt,primary->Pt()); - fhPi0PhiCharged ->Fill(phi,primary->Phi()); - fhPi0EtaCharged ->Fill(eta,primary->Eta()); - } - } - else{ - fhGamE ->Fill(e,primary->Energy()); - fhGamPt ->Fill(pt,primary->Pt()); - fhGamEta ->Fill(eta,primary->Eta()); - fhGamPhi ->Fill(phi,primary->Phi()); - fhGamDeltaE ->Fill(primary->Energy()-e); - fhGamDeltaPt ->Fill(primary->Pt()-pt); - fhGamDeltaPhi->Fill(primary->Phi()-phi); - fhGamDeltaEta->Fill(primary->Eta()-eta); - if(primary->Energy() > 0) fhGamRatioE ->Fill(e/primary->Energy()); - if(primary->Pt() > 0) fhGamRatioPt ->Fill(pt/primary->Pt()); - if(primary->Phi() > 0) fhGamRatioPhi->Fill(phi/primary->Phi()); - if(primary->Eta() > 0) fhGamRatioEta->Fill(eta/primary->Eta()); - if( ntracksmatched > 0){ - fhGamECharged ->Fill(e,primary->Energy()); - fhGamPtCharged ->Fill(pt,primary->Pt()); - fhGamPhiCharged ->Fill(phi,primary->Phi()); - fhGamEtaCharged ->Fill(eta,primary->Eta()); - } - } - }//pdg == 22 - else if(TMath::Abs(pdg) == 11) { - fhEleE ->Fill(e,primary->Energy()); - fhElePt ->Fill(pt,primary->Pt()); - fhEleEta ->Fill(eta,primary->Eta()); - fhElePhi ->Fill(phi,primary->Phi()); - fhEMVxyz ->Fill(vx,vy);//,vz); - fhEMR ->Fill(e,r); - if( ntracksmatched > 0){ - fhEleECharged ->Fill(e,primary->Energy()); - fhElePtCharged ->Fill(pt,primary->Pt()); - fhElePhiCharged ->Fill(phi,primary->Phi()); - fhEleEtaCharged ->Fill(eta,primary->Eta()); - } - } - else if(pdg == 111) { - fhPi0E ->Fill(e,primary->Energy()); - fhPi0Pt ->Fill(pt,primary->Pt()); - fhPi0Eta ->Fill(eta,primary->Eta()); - fhPi0Phi ->Fill(phi,primary->Phi()); - if( ntracksmatched > 0){ - fhPi0ECharged ->Fill(e,primary->Energy()); - fhPi0PtCharged ->Fill(pt,primary->Pt()); - fhPi0PhiCharged ->Fill(phi,primary->Phi()); - fhPi0EtaCharged ->Fill(eta,primary->Eta()); - } - } - else if(charge == 0){ - fhNeHadE ->Fill(e,primary->Energy()); - fhNeHadPt ->Fill(pt,primary->Pt()); - fhNeHadEta ->Fill(eta,primary->Eta()); - fhNeHadPhi ->Fill(phi,primary->Phi()); - fhHaVxyz ->Fill(vx,vy);//,vz); - fhHaR ->Fill(e,r); - if( ntracksmatched > 0){ - fhNeHadECharged ->Fill(e,primary->Energy()); - fhNeHadPtCharged ->Fill(pt,primary->Pt()); - fhNeHadPhiCharged ->Fill(phi,primary->Phi()); - fhNeHadEtaCharged ->Fill(eta,primary->Eta()); + //Invariant mass + if (nCaloClusters > 1 ) { + for(Int_t jclus = iclus + 1 ; jclus < nCaloClusters-1 ; jclus++) { + if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){ + AliESDCaloCluster* clus2 = (AliESDCaloCluster*) (caloClusters->At(jclus)); + //Get cluster kinematics + clus2->GetMomentum(mom2,v); } - } - else if(charge!=0){ - fhChHadE ->Fill(e,primary->Energy()); - fhChHadPt ->Fill(pt,primary->Pt()); - fhChHadEta ->Fill(eta,primary->Eta()); - fhChHadPhi ->Fill(phi,primary->Phi()); - fhHaVxyz ->Fill(vx,vy);//,vz); - fhHaR ->Fill(e,r); - if( ntracksmatched > 0){ - fhChHadECharged ->Fill(e,primary->Energy()); - fhChHadPtCharged ->Fill(pt,primary->Pt()); - fhChHadPhiCharged ->Fill(phi,primary->Phi()); - fhChHadEtaCharged ->Fill(eta,primary->Eta()); + if(GetReader()->GetDataType()==AliCaloTrackReader::kAOD){ + AliAODCaloCluster* clus2 = (AliAODCaloCluster*) (caloClusters->At(jclus)); + //Get cluster kinematics + clus2->GetMomentum(mom2,v); } - } - }//Work with stack also - - //matched cluster with tracks - if( ntracksmatched > 0){ - fhECharged ->Fill(e); - fhPtCharged ->Fill(pt); - fhPhiCharged ->Fill(phi); - fhEtaCharged ->Fill(eta); - fhEtaPhiCharged->Fill(eta,phi); - if((!strcmp(GetReader()->GetInputEvent()->GetName(),"AliESDEvent"))) { - AliESDEvent *esd = (AliESDEvent*) GetReader()->GetInputEvent(); - AliESDCaloCluster * esdcalo = (AliESDCaloCluster*) esd->GetCaloCluster(calo->GetID()); - Int_t trackIndex = esdcalo->GetTrackMatched(); - //printf("track index %d ntracks %d\n", trackIndex, esd->GetNumberOfTracks()); - if(trackIndex >= 0){ - AliESDtrack* track = (AliESDtrack*) esd->GetTrack(trackIndex); - if (track && track->GetOuterParam() ) { - - Double_t tphi = track->GetOuterParam()->Phi(); - Double_t teta = track->GetOuterParam()->Eta(); - Double_t tmom = track->GetOuterParam()->P(); - - Double_t deta = teta - eta; - Double_t dphi = tphi - phi; - if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); - if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); - Double_t dR = sqrt(dphi*dphi + deta*deta); - - Double_t pOverE = tmom/e; - - fh1pOverE->Fill(pOverE); - if(dR < 0.02) fh1pOverER02->Fill(pOverE); - - fh1dR->Fill(dR); - fh2MatchdEdx->Fill(track->P(),track->GetTPCsignal()); - - if(IsDataMC() && primary){ - Int_t pdg = primary->GetPdgCode(); - Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); - - if(TMath::Abs(pdg) == 11){ - fhMCEle1pOverE->Fill(pOverE); - fhMCEle1dR->Fill(dR); - fhMCEle2MatchdEdx->Fill(track->P(),track->GetTPCsignal()); - if(dR < 0.02) fhMCEle1pOverER02->Fill(pOverE); - } - else if(charge!=0){ - fhMCChHad1pOverE->Fill(pOverE); - fhMCChHad1dR->Fill(dR); - fhMCChHad2MatchdEdx->Fill(track->P(),track->GetTPCsignal()); - if(dR < 0.02) fhMCChHad1pOverER02->Fill(pOverE); - } - else if(charge == 0){ - fhMCNeutral1pOverE->Fill(pOverE); - fhMCNeutral1dR->Fill(dR); - fhMCNeutral2MatchdEdx->Fill(track->P(),track->GetTPCsignal()); - if(dR < 0.02) fhMCNeutral1pOverER02->Fill(pOverE); - } - } - int nITS = track->GetNcls(0); - int nTPC = track->GetNcls(1); - if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5 - && calo->GetNCells() > 1 && nITS > 3 && nTPC > 20) { - fh2EledEdx->Fill(track->P(),track->GetTPCsignal()); - } - } - else if(!track->GetOuterParam()){ - ULong_t status=AliESDtrack::kTPCrefit; - status|=AliESDtrack::kITSrefit; - //printf("track status %d\n", track->GetStatus() ); - fhEChargedNoOut ->Fill(e); - fhPtChargedNoOut ->Fill(pt); - fhPhiChargedNoOut ->Fill(phi); - fhEtaChargedNoOut ->Fill(eta); - fhEtaPhiChargedNoOut ->Fill(eta,phi); - if(GetDebug() >= 0 && ((track->GetStatus() & status) == status)) printf("ITS+TPC\n"); - } - else { - if(GetDebug() >= 0) printf("ERROR: Could not receive track %d\n", trackIndex); - } - }// non negative track index - }//do only if input are ESDs - }// at least one track matched - - //Invariant mass - if (nclusters > 1 ) { - for(Int_t jclus = iclus + 1 ; jclus < nclusters ; jclus++) { - AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (partList->At(jclus)); - //Get cluster kinematics - calo2->GetMomentum(mom2,v); + fhIM ->Fill((mom+mom2).E(),(mom+mom2).M()); - fhAsym->Fill((mom+mom2).E(),TMath::Abs((e-mom2.E())/(e+mom2.E()))); + fhAsym->Fill((mom+mom2).E(),TMath::Abs((mom.E()-mom2.E())/(mom.E()+mom2.E()))); + }// 2nd cluster loop - }////more than 1 cluster in calorimeter + }////more than 1 cluster in calorimeter + }//cluster loop + + //CaloCells + if(GetReader()->GetDataType()==AliCaloTrackReader::kESD){ + AliESDCaloCells * cell = 0x0; + Int_t ncells = 0; + if(fCalorimeter == "PHOS") cell = ((AliESDEvent*)GetReader()->GetInputEvent())->GetPHOSCells(); + else cell = ((AliESDEvent*)GetReader()->GetInputEvent())->GetEMCALCells(); + + if(!cell) { + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis"); + abort(); + } + + ncells = cell->GetNumberOfCells() ; + fhNCells->Fill(ncells) ; + if(GetDebug() > 0) + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells); + + for (Int_t iCell = 0; iCell < ncells; iCell++) { + if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell)); + fhAmplitude->Fill(cell->GetAmplitude(iCell)); + } + + }//ESD + else{//AOD + AliESDCaloCells * cell = 0x0; + Int_t ncells = 0; - //Cells per cluster - fhNCellsPerCluster->Fill(e,calo->GetNCells()); + if(fCalorimeter == "PHOS") cell = (AliESDCaloCells *) ((AliAODEvent*)GetReader()->GetInputEvent())->GetPHOSCells(); + else cell = (AliESDCaloCells *) ((AliAODEvent*)GetReader()->GetInputEvent())->GetEMCALCells(); - }// 1st cluster loop + if(!cell) { + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis"); + abort(); + } + + ncells = cell->GetNumberOfCells() ; + fhNCells->Fill(ncells) ; + if(GetDebug() > 0) + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), ncells); - //CaloCells - AliAODCaloCells * cell = new AliAODCaloCells ; - if(fCalorimeter == "PHOS") - cell = (AliAODCaloCells *) GetPHOSCells(); - else - cell = (AliAODCaloCells *) GetEMCALCells(); - - if(!cell) { - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - STOP: No CELLS available for analysis"); - abort(); - } + for (Int_t iCell = 0; iCell < ncells; iCell++) { + if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell)); + fhAmplitude->Fill(cell->GetAmplitude(iCell)); + } - //Some prints - if(GetDebug() > 0 && cell ) - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - In ESD %s cell entries %d\n", fCalorimeter.Data(), cell->GetNumberOfCells()); + }//AOD - Int_t ncells = cell->GetNumberOfCells() ; - fhNCells->Fill(ncells) ; + if(GetDebug() > 0) + printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n"); +} + +//__________________________________ +void AliAnaCalorimeterQA::ClusterHistograms(const TLorentzVector mom, const Int_t nCaloCellsPerCluster, + const Int_t nTracksMatched, const TObject * track, + const Int_t * labels, const Int_t nLabels){ + //Fill CaloCluster related histograms - for (Int_t iCell = 0; iCell < ncells; iCell++) { - if(GetDebug() > 0) printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", cell->GetAmplitude(iCell), cell->GetCellNumber(iCell)); - fhAmplitude->Fill(cell->GetAmplitude(iCell)); + AliAODMCParticle * aodprimary = 0x0; + TParticle * primary = 0x0; + Int_t tag = 0; + + Float_t e = mom.E(); + Float_t pt = mom.Pt(); + Float_t eta = mom.Eta(); + Float_t phi = mom.Phi(); + if(phi < 0) phi +=TMath::TwoPi(); + if(GetDebug() > 0) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - cluster: E %2.3f, pT %2.3f, eta %2.3f, phi %2.3f \n",e,pt,eta,phi*TMath::RadToDeg()); + if(IsDataMC()) { + printf("\t Primaries: nlabels %d, labels pointer %p\n",nLabels,labels); + if(!nLabels || !labels) printf("\t Strange, no labels!!!\n"); + } } - - //Monte Carlo - if(IsDataMC()){ + + fhE ->Fill(e); + fhPt ->Fill(pt); + fhPhi ->Fill(phi); + fhEta ->Fill(eta); + fhEtaPhi->Fill(eta,phi); + //Cells per cluster + fhNCellsPerCluster->Fill(e, nCaloCellsPerCluster); + + //Fill histograms only possible when simulation + if(IsDataMC() && nLabels > 0 && labels){ + //Play with the MC stack if available - for(Int_t i=8 ; iGetNprimary(); i++){ - primary = stack->Particle(i) ; - - //if (!primary->IsPrimary()) continue; - if (TMath::Abs(primary->Eta()) > 1) continue; + Int_t label = labels[0]; + + if(label < 0) { + if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** bad label ***: label %d \n", label); + return; + } + + Int_t pdg =-1; Int_t pdg0 =-1;Int_t status = -1; Int_t iMother = -1; Int_t iParent = -1; + Float_t vxMC= 0; Float_t vyMC = 0; + Float_t eMC = 0; Float_t ptMC= 0; Float_t phiMC =0; Float_t etaMC = 0; + Int_t charge = 0; + + //Check the origin. + tag = GetMCAnalysisUtils()->CheckOrigin(labels,nLabels, GetReader(),0); + + if(GetReader()->ReadStack() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){ //it MC stack and known tag + + if( label >= GetMCStack()->GetNtrack()) { + if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() *** large label ***: label %d, n tracks %d \n", label, GetMCStack()->GetNtrack()); + return ; + } - Int_t kf = primary->GetPdgCode(); - //printf("kf %d\n",kf); + primary = GetMCStack()->Particle(label); + iMother = label; + pdg0 = TMath::Abs(primary->GetPdgCode()); + pdg = pdg0; + status = primary->GetStatusCode(); + vxMC = primary->Vx(); + vyMC = primary->Vy(); + iParent = primary->GetFirstMother(); + + if(GetDebug() > 1 ) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n"); + printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg0, primary->GetName(),status, iParent); + } + + //Get final particle, no conversion products + if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){ + //Get the parent + primary = GetMCStack()->Particle(iParent); + pdg = TMath::Abs(primary->GetPdgCode()); + if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n"); + while((pdg == 22 || pdg == 11) && status != 1){ + iMother = iParent; + primary = GetMCStack()->Particle(iMother); + status = primary->GetStatusCode(); + iParent = primary->GetFirstMother(); + pdg = TMath::Abs(primary->GetPdgCode()); + if(GetDebug() > 1 )printf("\t pdg %d, index %d, %s, status %d \n",pdg, iMother, primary->GetName(),status); + } + + if(GetDebug() > 1 ) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n"); + printf("\t Mother label %d, pdg %d, %s, status %d, parent %d \n",iMother, pdg, primary->GetName(), status, iParent); + } + + } + + //Overlapped pi0 (or eta, there will be very few), get the meson + if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || + GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){ + if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: \n"); + while(pdg != 111 && pdg != 221){ + iMother = iParent; + primary = GetMCStack()->Particle(iMother); + status = primary->GetStatusCode(); + iParent = primary->GetFirstMother(); + pdg = TMath::Abs(primary->GetPdgCode()); + if(GetDebug() > 1 ) printf("\t pdg %d, %s, index %d\n",pdg, primary->GetName(),iMother); + if(iMother==-1) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n"); + //break; + } + } + + if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", + primary->GetName(),iMother); + } + + eMC = primary->Energy(); + ptMC = primary->Pt(); + phiMC = primary->Phi(); + etaMC = primary->Eta(); + pdg = TMath::Abs(primary->GetPdgCode()); + charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); + + } + else if(GetReader()->ReadAODMCParticles() && !GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCUnknown)){//it MC AOD and known tag + //Get the list of MC particles + if(!GetReader()->GetAODMCParticles(0)) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - MCParticles not available!\n"); + abort(); + } - Bool_t in = kTRUE; - primary->Momentum(mom); - if(IsFidutialCutOn()) in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ; + aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(0))->At(label); + iMother = label; + pdg0 = TMath::Abs(aodprimary->GetPdgCode()); + pdg = pdg0; + status = aodprimary->IsPrimary(); + vxMC = aodprimary->Xv(); + vyMC = aodprimary->Yv(); + iParent = aodprimary->GetMother(); + + if(GetDebug() > 1 ) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Cluster most contributing mother: \n"); + printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d, parent %d \n", + iMother, pdg0, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary(), iParent); + } - if (kf==22) { - fhGenGamPt ->Fill(primary->Pt()); - fhGenGamEta->Fill(primary->Eta()); - fhGenGamPhi->Fill(primary->Phi()); - if(in){ - fhGenGamAccE ->Fill(primary->Energy()); - fhGenGamAccPt ->Fill(primary->Pt()); - fhGenGamAccEta->Fill(primary->Eta()); - fhGenGamAccPhi->Fill(primary->Phi()); + //Get final particle, no conversion products + if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCConversion)){ + if(GetDebug() > 1 ) + printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted cluster!. Find before conversion: \n"); + //Get the parent + aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iParent); + pdg = TMath::Abs(aodprimary->GetPdgCode()); + while ((pdg == 22 || pdg == 11) && !aodprimary->IsPhysicalPrimary()) { + iMother = iParent; + aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother); + status = aodprimary->IsPrimary(); + iParent = aodprimary->GetMother(); + pdg = TMath::Abs(aodprimary->GetPdgCode()); + if(GetDebug() > 1 ) + printf("\t pdg %d, index %d, Primary? %d, Physical Primary? %d \n", + pdg, iMother, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()); + } + + if(GetDebug() > 1 ) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Converted Cluster mother before conversion: \n"); + printf("\t Mother label %d, pdg %d, parent %d, Primary? %d, Physical Primary? %d \n", + iMother, pdg, iParent, aodprimary->IsPrimary(), aodprimary->IsPhysicalPrimary()); } + } - else if (kf==111) { - fhGenPi0Pt ->Fill(primary->Pt()); - fhGenPi0Eta->Fill(primary->Eta()); - fhGenPi0Phi->Fill(primary->Phi()); - if(in){ - fhGenPi0AccE ->Fill(primary->Energy()); - fhGenPi0AccPt ->Fill(primary->Pt()); - fhGenPi0AccEta->Fill(primary->Eta()); - fhGenPi0AccPhi->Fill(primary->Phi()); + + //Overlapped pi0 (or eta, there will be very few), get the meson + if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || + GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){ + if(GetDebug() > 1 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped Meson decay!, Find it: PDG %d, mom %d \n",pdg, iMother); + while(pdg != 111 && pdg != 221){ + + iMother = iParent; + aodprimary = (AliAODMCParticle*)(GetReader()->GetAODMCParticles(0))->At(iMother); + status = aodprimary->IsPrimary(); + iParent = aodprimary->GetMother(); + pdg = TMath::Abs(aodprimary->GetPdgCode()); + + if(GetDebug() > 1 ) printf("\t pdg %d, index %d\n",pdg, iMother); + + if(iMother==-1) { + printf("AliAnaCalorimeterQA::ClusterHistograms() - Tagged as Overlapped photon but meson not found, why?\n"); + //break; + } + } + + if(GetDebug() > 2 ) printf("AliAnaCalorimeterQA::ClusterHistograms() - Overlapped %s decay, label %d \n", + aodprimary->GetName(),iMother); + } + + status = aodprimary->IsPrimary(); + eMC = aodprimary->E(); + ptMC = aodprimary->Pt(); + phiMC = aodprimary->Phi(); + etaMC = aodprimary->Eta(); + pdg = TMath::Abs(aodprimary->GetPdgCode()); + charge = aodprimary->Charge(); + + } + + //Float_t vz = primary->Vz(); + Float_t r = TMath::Sqrt(vxMC*vxMC + vyMC*vyMC); + if((pdg == 22 || TMath::Abs(pdg)==11) && status!=1) { + fhEMVxyz ->Fill(vxMC,vyMC);//,vz); + fhEMR ->Fill(e,r); + } + + //printf("reco e %f, pt %f, phi %f, eta %f \n", e, pt, phi, eta); + //printf("prim e %f, pt %f, phi %f, eta %f \n", eMC,ptMC,phiMC ,etaMC ); + //printf("vertex: vx %f, vy %f, vz %f, r %f \n", vxMC, vyMC, vz, r); + + + fh2E ->Fill(e, eMC); + fh2Pt ->Fill(pt, ptMC); + fh2Phi ->Fill(phi, phiMC); + fh2Eta ->Fill(eta, etaMC); + fhDeltaE ->Fill(eMC-e); + fhDeltaPt ->Fill(ptMC-pt); + fhDeltaPhi->Fill(phiMC-phi); + fhDeltaEta->Fill(etaMC-eta); + if(eMC > 0) fhRatioE ->Fill(e/eMC); + if(ptMC > 0) fhRatioPt ->Fill(pt/ptMC); + if(phiMC > 0) fhRatioPhi->Fill(phi/phiMC); + if(etaMC > 0) fhRatioEta->Fill(eta/etaMC); + + + //Overlapped pi0 (or eta, there will be very few) + if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPi0) || + GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCEta)){ + //cout<<"Fill pi0"<< "E "<< e <<" prim E "<Fill(e,eMC); + fhPi0Pt ->Fill(pt,ptMC); + fhPi0Eta ->Fill(eta,etaMC); + fhPi0Phi ->Fill(phi,phiMC); + if( nTracksMatched > 0){ + fhPi0ECharged ->Fill(e,eMC); + fhPi0PtCharged ->Fill(pt,ptMC); + fhPi0PhiCharged ->Fill(phi,phiMC); + fhPi0EtaCharged ->Fill(eta,etaMC); } + }//Overlapped pizero decay + else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCPhoton)){ + fhGamE ->Fill(e,eMC); + fhGamPt ->Fill(pt,ptMC); + fhGamEta ->Fill(eta,etaMC); + fhGamPhi ->Fill(phi,phiMC); + fhGamDeltaE ->Fill(eMC-e); + fhGamDeltaPt ->Fill(ptMC-pt); + fhGamDeltaPhi->Fill(phiMC-phi); + fhGamDeltaEta->Fill(etaMC-eta); + if(eMC > 0) fhGamRatioE ->Fill(e/eMC); + if(ptMC > 0) fhGamRatioPt ->Fill(pt/ptMC); + if(phiMC > 0) fhGamRatioPhi->Fill(phi/phiMC); + if(etaMC > 0) fhGamRatioEta->Fill(eta/etaMC); + if( nTracksMatched > 0){ + fhGamECharged ->Fill(e,eMC); + fhGamPtCharged ->Fill(pt,ptMC); + fhGamPhiCharged ->Fill(phi,phiMC); + fhGamEtaCharged ->Fill(eta,etaMC); + } + }//photon + else if(GetMCAnalysisUtils()->CheckTagBit(tag, AliMCAnalysisUtils::kMCElectron)){ + fhEleE ->Fill(e,eMC); + fhElePt ->Fill(pt,ptMC); + fhEleEta ->Fill(eta,etaMC); + fhElePhi ->Fill(phi,phiMC); + fhEMVxyz ->Fill(vxMC,vyMC);//,vz); + fhEMR ->Fill(e,r); + if( nTracksMatched > 0){ + fhEleECharged ->Fill(e,eMC); + fhElePtCharged ->Fill(pt,ptMC); + fhElePhiCharged ->Fill(phi,phiMC); + fhEleEtaCharged ->Fill(eta,etaMC); + } + } + else if(charge == 0){ + fhNeHadE ->Fill(e,eMC); + fhNeHadPt ->Fill(pt,ptMC); + fhNeHadEta ->Fill(eta,etaMC); + fhNeHadPhi ->Fill(phi,phiMC); + fhHaVxyz ->Fill(vxMC,vyMC);//,vz); + fhHaR ->Fill(e,r); + if( nTracksMatched > 0){ + fhNeHadECharged ->Fill(e,eMC); + fhNeHadPtCharged ->Fill(pt,ptMC); + fhNeHadPhiCharged ->Fill(phi,phiMC); + fhNeHadEtaCharged ->Fill(eta,etaMC); } - else if (kf==221) { - fhGenEtaPt ->Fill(primary->Pt()); - fhGenEtaEta->Fill(primary->Eta()); - fhGenEtaPhi->Fill(primary->Phi()); + } + else if(charge!=0){ + fhChHadE ->Fill(e,eMC); + fhChHadPt ->Fill(pt,ptMC); + fhChHadEta ->Fill(eta,etaMC); + fhChHadPhi ->Fill(phi,phiMC); + fhHaVxyz ->Fill(vxMC,vyMC);//,vz); + fhHaR ->Fill(e,r); + if( nTracksMatched > 0){ + fhChHadECharged ->Fill(e,eMC); + fhChHadPtCharged ->Fill(pt,ptMC); + fhChHadPhiCharged ->Fill(phi,phiMC); + fhChHadEtaCharged ->Fill(eta,etaMC); } - else if (kf==223) { - fhGenOmegaPt ->Fill(primary->Pt()); - fhGenOmegaEta->Fill(primary->Eta()); - fhGenOmegaPhi->Fill(primary->Phi()); + } + }//Work with MC + + + //Match tracks and clusters + //To be Modified in case of AODs + + //if(ntracksmatched==1 && trackIndex==-1) ntracksmatched=0; + + if( nTracksMatched > 0){ + fhECharged ->Fill(e); + fhPtCharged ->Fill(pt); + fhPhiCharged ->Fill(phi); + fhEtaCharged ->Fill(eta); + fhEtaPhiCharged->Fill(eta,phi); + + //printf("track index %d ntracks %d\n", esd->GetNumberOfTracks()); + if(GetReader()->GetDataType()==AliCaloTrackReader::kESD) { + if (track && ((AliESDtrack*)track)->GetOuterParam() ) { + + Double_t tphi = ((AliESDtrack*)track)->GetOuterParam()->Phi(); + Double_t teta = ((AliESDtrack*)track)->GetOuterParam()->Eta(); + Double_t tmom = ((AliESDtrack*)track)->GetOuterParam()->P(); + + Double_t tmom2 = ((AliESDtrack*)track)->P(); + Double_t tpcSignal = ((AliESDtrack*)track)->GetTPCsignal(); + + Double_t deta = teta - eta; + Double_t dphi = tphi - phi; + if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi(); + if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi(); + Double_t dR = sqrt(dphi*dphi + deta*deta); + + Double_t pOverE = tmom/e; + + fh1pOverE->Fill(pOverE); + if(dR < 0.02) fh1pOverER02->Fill(pOverE); + + fh1dR->Fill(dR); + fh2MatchdEdx->Fill(tmom2,tpcSignal); + + if(IsDataMC() && primary){ + Int_t pdg = primary->GetPdgCode(); + Double_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge(); + + if(TMath::Abs(pdg) == 11){ + fhMCEle1pOverE->Fill(pOverE); + fhMCEle1dR->Fill(dR); + fhMCEle2MatchdEdx->Fill(tmom2,tpcSignal); + if(dR < 0.02) fhMCEle1pOverER02->Fill(pOverE); + } + else if(charge!=0){ + fhMCChHad1pOverE->Fill(pOverE); + fhMCChHad1dR->Fill(dR); + fhMCChHad2MatchdEdx->Fill(tmom2,tpcSignal); + if(dR < 0.02) fhMCChHad1pOverER02->Fill(pOverE); + } + else if(charge == 0){ + fhMCNeutral1pOverE->Fill(pOverE); + fhMCNeutral1dR->Fill(dR); + fhMCNeutral2MatchdEdx->Fill(tmom2,tpcSignal); + if(dR < 0.02) fhMCNeutral1pOverER02->Fill(pOverE); + } + }//DataMC + int nITS = ((AliESDtrack*)track)->GetNcls(0); + int nTPC = ((AliESDtrack*)track)->GetNcls(1); + if(dR < 0.02 && pOverE > 0.5 && pOverE < 1.5 + && nCaloCellsPerCluster > 1 && nITS > 3 && nTPC > 20) { + fh2EledEdx->Fill(tmom2,tpcSignal); + } + }//Outer param available + else if(!((AliESDtrack*)track)->GetOuterParam()){ + ULong_t status=AliESDtrack::kTPCrefit; + status|=AliESDtrack::kITSrefit; + //printf("track status %d\n", track->GetStatus() ); + fhEChargedNoOut ->Fill(e); + fhPtChargedNoOut ->Fill(pt); + fhPhiChargedNoOut ->Fill(phi); + fhEtaChargedNoOut ->Fill(eta); + fhEtaPhiChargedNoOut ->Fill(eta,phi); + if(GetDebug() >= 0 && ((((AliESDtrack*)track)->GetStatus() & status) == status)) printf("ITS+TPC\n"); } - else if (TMath::Abs(kf)==11) { - fhGenElePt ->Fill(primary->Pt()); - fhGenEleEta->Fill(primary->Eta()); - fhGenElePhi->Fill(primary->Phi()); + else { + if(GetDebug() >= 0) printf("AliAnaCalorimeterQA::ClusterHistograms() ERROR: Could not receive track %d\n", ((AliESDtrack*)track)->GetID()); } - } //primary loop - } //Is data MC + }//do only if input are ESDs + } - if(GetDebug() > 0) - printf("AliAnaCalorimeterQA::MakeAnalysisFillHistograms() - End \n"); +}// Clusters +//______________________________________________________________________________ +void AliAnaCalorimeterQA::MCHistograms(const TLorentzVector mom, const Int_t pdg){ + //Fill pure monte carlo related histograms -} + Float_t eMC = mom.E(); + Float_t ptMC = mom.Pt(); + Float_t phiMC = mom.Phi(); + if(phiMC < 0) + phiMC += TMath::TwoPi(); + Float_t etaMC = mom.Eta(); + + if (TMath::Abs(etaMC) > 1) return; + Bool_t in = kTRUE; + if(IsFidutialCutOn()) in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ; + + if (pdg==22) { + fhGenGamPt ->Fill(ptMC); + fhGenGamEta->Fill(etaMC); + fhGenGamPhi->Fill(phiMC); + if(in){ + fhGenGamAccE ->Fill(eMC); + fhGenGamAccPt ->Fill(ptMC); + fhGenGamAccEta->Fill(etaMC); + fhGenGamAccPhi->Fill(phiMC); + } + } + else if (pdg==111) { + fhGenPi0Pt ->Fill(ptMC); + fhGenPi0Eta->Fill(etaMC); + fhGenPi0Phi->Fill(phiMC); + if(in){ + fhGenPi0AccE ->Fill(eMC); + fhGenPi0AccPt ->Fill(ptMC); + fhGenPi0AccEta->Fill(etaMC); + fhGenPi0AccPhi->Fill(phiMC); + } + } + else if (pdg==221) { + fhGenEtaPt ->Fill(ptMC); + fhGenEtaEta->Fill(etaMC); + fhGenEtaPhi->Fill(phiMC); + } + else if (pdg==223) { + fhGenOmegaPt ->Fill(ptMC); + fhGenOmegaEta->Fill(etaMC); + fhGenOmegaPhi->Fill(phiMC); + } + else if (TMath::Abs(pdg)==11) { + fhGenElePt ->Fill(ptMC); + fhGenEleEta->Fill(etaMC); + fhGenElePhi->Fill(phiMC); + } + +} + //________________________________________________________________________ void AliAnaCalorimeterQA::ReadHistograms(TList* outputList) { @@ -1492,7 +1700,14 @@ void AliAnaCalorimeterQA::ReadHistograms(TList* outputList) //__________________________________________________________________ void AliAnaCalorimeterQA::Terminate(TList* outputList) { - //Do plots if requested + //Do plots if requested + char line[1024] ; + + //if(fRemoveOutputAOD){ + // sprintf(line, ".!rm -fR %s",((AliVEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()))->GetOutputFileName()); + // gROOT->ProcessLine(line); + //} + if(!fMakePlots) return; //Do some plots to end @@ -2801,7 +3016,7 @@ void AliAnaCalorimeterQA::Terminate(TList* outputList) } - char line[1024] ; + sprintf(line, ".!tar -zcf QA_%s_%s.tar.gz *%s*.eps", fCalorimeter.Data(), GetName(),fCalorimeter.Data()) ; gROOT->ProcessLine(line); sprintf(line, ".!rm -fR *.eps"); diff --git a/PWG4/PartCorrDep/AliAnaCalorimeterQA.h b/PWG4/PartCorrDep/AliAnaCalorimeterQA.h index 675bb2ff28c..766e41f2db0 100755 --- a/PWG4/PartCorrDep/AliAnaCalorimeterQA.h +++ b/PWG4/PartCorrDep/AliAnaCalorimeterQA.h @@ -13,7 +13,6 @@ // --- Root system --- class TH2F; class TH1F; -//class TH3F; // --- Analysis system --- #include "AliAnaPartCorrBaseClass.h" @@ -27,15 +26,21 @@ class AliAnaCalorimeterQA : public AliAnaPartCorrBaseClass { AliAnaCalorimeterQA & operator = (const AliAnaCalorimeterQA & g) ;//cpy assignment virtual ~AliAnaCalorimeterQA() {;} //virtual dtor + void ClusterHistograms(const TLorentzVector mom, const Int_t nCaloCellsPerCluster, + const Int_t nTracksMatched, const TObject* track, + const Int_t * labels, const Int_t nLabels); + TList * GetCreateOutputObjects(); - + void Init(); void InitParameters(); void Print(const Option_t * opt) const; void MakeAnalysisFillHistograms() ; - + + void MCHistograms(const TLorentzVector mom, const Int_t pdg); + TString GetCalorimeter() const {return fCalorimeter ;} void SetCalorimeter( TString calo ) {fCalorimeter = calo; } TString GetStyleMacro() const {return fStyleMacro ;} @@ -43,16 +48,16 @@ class AliAnaCalorimeterQA : public AliAnaPartCorrBaseClass { void SwithOnPlotsMaking() {fMakePlots = kTRUE;} void SwithOffPlotsMaking() {fMakePlots = kFALSE;} - + void Terminate(TList * outputList); void ReadHistograms(TList * outputList); //Fill histograms with histograms in ouput list, needed in Terminate. private: - TString fCalorimeter ; //Calorimeter selection - TString fStyleMacro ; //Location of macro for plots style - Bool_t fMakePlots ; //Print plots - + TString fCalorimeter ; //Calorimeter selection + TString fStyleMacro ; //Location of macro for plots style + Bool_t fMakePlots ; //Print plots + //Histograms //CaloClusters TH1F * fhE ; //! E distribution, Reco diff --git a/PWG4/macros/AddTaskCalorimeterQA.C b/PWG4/macros/AddTaskCalorimeterQA.C index 42af4390092..633a3950414 100644 --- a/PWG4/macros/AddTaskCalorimeterQA.C +++ b/PWG4/macros/AddTaskCalorimeterQA.C @@ -22,17 +22,13 @@ AliAnalysisTaskParticleCorrelation *AddTaskCalorimeterQA(TString data, Bool_t kU //=========================================================================== //Reader + //For this particular analysis few things done by the reader. + //Nothing else needs to be set. AliCaloTrackReader * reader = 0x0; - if(data=="AOD") reader = new AliCaloTrackAODReader(); + if(data=="AOD") reader = new AliCaloTrackAODReader(); else if(data=="ESD") reader = new AliCaloTrackESDReader(); - reader->SetDebug(-1);//10 for lots of messages - reader->SwitchOnEMCALCells(); - reader->SwitchOnPHOSCells(); - //Min particle pT - reader->SetEMCALPtMin(0.2); - reader->SetPHOSPtMin(0.2); - reader->SetCTSPtMin(0.2); - if(kPrintSettings) reader->Print(""); + //reader->SetDebug(10);//10 for lots of messages + if(kPrintSettings) reader->Print(""); if(kUseKinematics){ if(inputDataType == "ESD"){ @@ -44,10 +40,9 @@ AliAnalysisTaskParticleCorrelation *AddTaskCalorimeterQA(TString data, Bool_t kU reader->SwitchOnAODMCParticles(); } } - reader->SetDeltaAODFileName(""); //Do not create deltaAOD file, this analysis do not create branches. - //Do not keep the temporal AOD in file. - //reader->SwitchOnCleanStdAOD(); + reader->SetDeltaAODFileName(""); //Do not create deltaAOD file, this analysis do not create branches. + // ##### Analysis algorithm settings #### @@ -57,24 +52,27 @@ AliAnalysisTaskParticleCorrelation *AddTaskCalorimeterQA(TString data, Bool_t kU fidCut->DoPHOSFidutialCut(kTRUE) ; AliAnaCalorimeterQA *emcalQA = new AliAnaCalorimeterQA(); - emcalQA->SetDebug(-1); //10 for lots of messages + //emcalQA->SetDebug(2); //10 for lots of messages emcalQA->SetCalorimeter("EMCAL"); - if(kUseKinematics && inputDataType!="AOD") emcalQA->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet. + if(kUseKinematics) emcalQA->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet. else emcalQA->SwitchOffDataMC() ; emcalQA->AddToHistogramsName("EMCAL_"); //Begining of histograms name emcalQA->SetFidutialCut(fidCut); emcalQA->SwitchOnFidutialCut(); if(kPrintSettings) emcalQA->Print(""); + //emcalQA->GetMCAnalysisUtils()->SetDebug(10); AliAnaCalorimeterQA *phosQA = new AliAnaCalorimeterQA(); - phosQA->SetDebug(-1); //10 for lots of messages + //phosQA->SetDebug(2); //10 for lots of messages phosQA->SetCalorimeter("PHOS"); - if(kUseKinematics && inputDataType!="AOD") phosQA->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet. + if(kUseKinematics) phosQA->SwitchOnDataMC() ;//Access MC stack and fill more histograms, AOD MC not implemented yet. else phosQA->SwitchOffDataMC() ; phosQA->AddToHistogramsName("PHOS_");//Begining of histograms name phosQA->SetFidutialCut(fidCut); phosQA->SwitchOnFidutialCut(); if(kPrintSettings)phosQA->Print(""); + //phosQA->GetMCAnalysisUtils()->SetDebug(10); + // #### Configure Maker #### AliAnaPartCorrMaker * maker = new AliAnaPartCorrMaker(); @@ -94,7 +92,7 @@ AliAnalysisTaskParticleCorrelation *AddTaskCalorimeterQA(TString data, Bool_t kU //=========================================================================== AliAnalysisTaskParticleCorrelation * task = new AliAnalysisTaskParticleCorrelation ("CalorimeterPerformance"); task->SetConfigFileName(""); //Don't configure the analysis via configuration file. - //task->SetDebugLevel(-1); + task->SetDebugLevel(0); task->SetAnalysisMaker(maker); mgr->AddTask(task); @@ -106,7 +104,7 @@ AliAnalysisTaskParticleCorrelation *AddTaskCalorimeterQA(TString data, Bool_t kU //============================================================================== mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer()); // AOD output slot will be used in a different way in future - mgr->ConnectOutput (task, 0, mgr->GetCommonOutputContainer()); + //mgr->ConnectOutput (task, 0, mgr->GetCommonOutputContainer()); mgr->ConnectOutput (task, 1, cout_pc); return task; diff --git a/PWG4/macros/ConfigAnalysisCalorimeterQA.C b/PWG4/macros/ConfigAnalysisCalorimeterQA.C index b229f36114a..e07b1f60697 100644 --- a/PWG4/macros/ConfigAnalysisCalorimeterQA.C +++ b/PWG4/macros/ConfigAnalysisCalorimeterQA.C @@ -3,7 +3,8 @@ //------------------------------------ // Configuration macro example: // -// Calorimeters QA +// Calorimeters QA: Validation of data (MC) +// Valid for ESDs, comment/uncomment below in the reader part for AODs // // Author : Gustavo Conesa Balbastre (INFN-LNF) //------------------------------------ @@ -22,35 +23,23 @@ AliAnaPartCorrMaker* ConfigAnalysis() //----------------------------------------------------------- // Reader //----------------------------------------------------------- + //For this particular analysis few things done by the reader. + //Nothing else needs to be set. AliCaloTrackESDReader *reader = new AliCaloTrackESDReader(); reader->SetDebug(-1); - - //Switch on or off the detectors information that you want - reader->SwitchOnEMCAL(); - reader->SwitchOffCTS(); - reader->SwitchOnPHOS(); - reader->SwitchOnEMCALCells(); - reader->SwitchOnPHOSCells(); - - //Min particle pT - reader->SetEMCALPtMin(0.); - reader->SetPHOSPtMin(0.); - - // //We want tracks fitted in the detectors: - // ULong_t status=AliAODTrack::kTPCrefit; - // status|=AliAODTrack::kITSrefit; //(default settings) - - // We want tracks whose PID bit is set: - // ULong_t status =AliAODTrack::kITSpid; - // status|=AliAODTrack::kTPCpid; - - // reader->SetTrackStatus(status); - - //Remove the temporal AODs we create. - reader->SwitchOnCleanStdAOD(); - + reader->SwitchOnStack(); + reader->SwitchOffAODMCParticles(); + reader->SetDeltaAODFileName(""); //Do not create deltaAOD file, this analysis do not create branches. reader->Print(""); + //For AODs: +// AliCaloTrackAODReader *reader = new AliCaloTrackAODReader(); +// reader->SetDebug(-1); +// reader->SwitchOffStack(); +// reader->SwitchOnAODMCParticles(); +// reader->SetDeltaAODFileName(""); //Do not create deltaAOD file, this analysis do not create branches. +// reader->Print(""); + //--------------------------------------------------------------------- // Analysis algorithm @@ -92,7 +81,7 @@ AliAnaPartCorrMaker* ConfigAnalysis() maker->AddAnalysis(anaPHOS,0); //maker->SetAnaDebug(0) ; maker->SwitchOnHistogramsMaker() ; - maker->SwitchOffAODsMaker() ; + maker->SwitchOffAODsMaker() ;//No AODs created in this task. maker->Print(""); // -- 2.43.5