From 95aee5e1cbb9d15d85f47d3185631b5dd81f5e77 Mon Sep 17 00:00:00 2001 From: gconesab Date: Thu, 7 Aug 2014 11:53:21 +0200 Subject: [PATCH] merge the 2 MCHistogram methods, simplify the filling of pure MC histograms, introduce the method to check the acceptance of the calorimeter in CaloUtils; move the Correlate method a bit before, check inside the existence of the calorimeter arrays --- .../AliAnaCalorimeterQA.cxx | 292 ++++++++++-------- .../AliAnaCalorimeterQA.h | 12 +- 2 files changed, 174 insertions(+), 130 deletions(-) diff --git a/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx b/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx index a069c705fb7..34af9409eb0 100755 --- a/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx +++ b/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx @@ -157,8 +157,8 @@ fhRecoMCDeltaE(), fhRecoMCRatioE(), fhRecoMCDeltaPhi(), fhRecoMCDeltaEta(), // MC only -fhGenMCE(), fhGenMCEtaPhi(), -fhGenMCAccE(), fhGenMCAccEtaPhi(), +fhGenMCE(), fhGenMCPt(), fhGenMCEtaPhi(), +fhGenMCAccE(), fhGenMCAccPt(), fhGenMCAccEtaPhi(), //matched MC fhEMVxyz(0), fhEMR(0), @@ -226,6 +226,16 @@ fhTrackMatchedDEtaPos(0), fhTrackMatchedDPhiPos(0), f fhRecoMCDeltaEta[i][0] = 0; fhRecoMCDeltaEta[i][1] = 0; } + for(Int_t i = 0; i < 4; i++) + { + fhGenMCE[i] = 0; + fhGenMCPt[i] = 0; + fhGenMCEtaPhi[i] = 0; + fhGenMCAccE[i] = 0; + fhGenMCAccPt[i] = 0; + fhGenMCAccEtaPhi[i] = 0; + } + //Initialize parameters InitParameters(); } @@ -756,7 +766,6 @@ void AliAnaCalorimeterQA::ClusterAsymmetryHistograms(AliVCluster* clus, Int_t ab if(goodCluster) { - // Was cluster matched? Bool_t matched = GetCaloPID()->IsTrackMatched(clus,GetCaloUtils(),GetReader()->GetInputEvent()); @@ -943,7 +952,6 @@ void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters, AliVCaloCells* cells) { // Fill clusters related histograms - TLorentzVector mom ; Int_t nLabel = 0 ; Int_t *labels = 0x0; @@ -956,7 +964,7 @@ void AliAnaCalorimeterQA::ClusterLoopHistograms(const TObjArray *caloClusters, // Get vertex for photon momentum calculation and event selection Double_t v[3] = {0,0,0}; //vertex ; //GetReader()->GetVertex(v); - + Int_t *nClustersInModule = new Int_t[fNModules]; for(Int_t imod = 0; imod < fNModules; imod++ ) nClustersInModule[imod] = 0; @@ -1542,10 +1550,29 @@ void AliAnaCalorimeterQA::Correlate() { // Correlate information from PHOS and EMCAL and with V0 and track multiplicity - //Clusters + //Clusters arrays TObjArray * caloClustersEMCAL = GetEMCALClusters(); TObjArray * caloClustersPHOS = GetPHOSClusters(); + if(!caloClustersEMCAL || !caloClustersPHOS) + { + if( GetDebug() > 0 ) printf("AliAnaCalorimeterQA::Correlate() - PHOS (%p) or EMCAL (%p) clusters array not available, do not correlate\n", + caloClustersPHOS,caloClustersEMCAL); + return ; + } + + //Cells arrays + AliVCaloCells * cellsEMCAL = GetEMCALCells(); + AliVCaloCells * cellsPHOS = GetPHOSCells(); + + if(!cellsEMCAL || !cellsPHOS) + { + if( GetDebug() > 0 ) printf("AliAnaCalorimeterQA::Correlate() - PHOS (%p) or EMCAL (%p) cells array ot available, do not correlate\n", + cellsPHOS,cellsEMCAL); + return ; + } + + // Clusters parameters Int_t nclEMCAL = caloClustersEMCAL->GetEntriesFast(); Int_t nclPHOS = caloClustersPHOS ->GetEntriesFast(); @@ -1560,11 +1587,7 @@ void AliAnaCalorimeterQA::Correlate() for(iclus = 0 ; iclus < caloClustersPHOS->GetEntriesFast(); iclus++) sumClusterEnergyPHOS += ((AliVCluster*)caloClustersPHOS->At(iclus))->E(); - //Cells - - AliVCaloCells * cellsEMCAL = GetEMCALCells(); - AliVCaloCells * cellsPHOS = GetPHOSCells(); - + //Cells parameters Int_t ncellsEMCAL = cellsEMCAL->GetNumberOfCells(); Int_t ncellsPHOS = cellsPHOS ->GetNumberOfCells(); @@ -2928,33 +2951,45 @@ TList * AliAnaCalorimeterQA::GetCreateOutputObjects() //Pure MC for(Int_t iPart = 0; iPart < 4; iPart++) { - fhGenMCE[iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) , + fhGenMCE [iPart] = new TH1F(Form("hGenMCE_%s",particleName[iPart].Data()) , + Form("#it{E} of generated %s",particleName[iPart].Data()), + nptbins,ptmin,ptmax); + + fhGenMCPt[iPart] = new TH1F(Form("hGenMCPt_%s",particleName[iPart].Data()) , Form("#it{p}_{T} of generated %s",particleName[iPart].Data()), nptbins,ptmin,ptmax); + fhGenMCEtaPhi[iPart] = new TH2F(Form("hGenMCEtaPhi_%s",particleName[iPart].Data()), Form("Y vs #phi of generated %s",particleName[iPart].Data()), - netabins,etamin,etamax,nphibins,phimin,phimax); + 200,-1,1,360,0,TMath::TwoPi()); - fhGenMCE[iPart] ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); + fhGenMCE [iPart] ->SetXTitle("#it{E} (GeV)"); + fhGenMCPt[iPart] ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhGenMCEtaPhi[iPart]->SetXTitle("#eta"); fhGenMCEtaPhi[iPart]->SetYTitle("#phi (rad)"); - - outputContainer->Add(fhGenMCE[iPart]); + + outputContainer->Add(fhGenMCE [iPart]); + outputContainer->Add(fhGenMCPt [iPart]); outputContainer->Add(fhGenMCEtaPhi[iPart]); - fhGenMCAccE[iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) , + fhGenMCAccE [iPart] = new TH1F(Form("hGenMCAccE_%s",particleName[iPart].Data()) , Form("#it{E} of generated %s",particleName[iPart].Data()), nptbins,ptmin,ptmax); + fhGenMCAccPt[iPart] = new TH1F(Form("hGenMCAccPt_%s",particleName[iPart].Data()) , + Form("#it{p}_{T} of generated %s",particleName[iPart].Data()), + nptbins,ptmin,ptmax); fhGenMCAccEtaPhi[iPart] = new TH2F(Form("hGenMCAccEtaPhi_%s",particleName[iPart].Data()), Form("Y vs #phi of generated %s",particleName[iPart].Data()), netabins,etamin,etamax,nphibins,phimin,phimax); - fhGenMCAccE[iPart] ->SetXTitle("#it{E} (GeV)"); + fhGenMCAccE [iPart] ->SetXTitle("#it{E} (GeV)"); + fhGenMCAccPt[iPart] ->SetXTitle("#it{p}_{T} (GeV/#it{c})"); fhGenMCAccEtaPhi[iPart]->SetXTitle("#eta"); fhGenMCAccEtaPhi[iPart]->SetYTitle("#phi (rad)"); - outputContainer->Add(fhGenMCAccE[iPart]); + outputContainer->Add(fhGenMCAccE [iPart]); + outputContainer->Add(fhGenMCAccPt [iPart]); outputContainer->Add(fhGenMCAccEtaPhi[iPart]); } @@ -3335,6 +3370,9 @@ void AliAnaCalorimeterQA::MakeAnalysisFillHistograms() //Play with the MC stack if available if(IsDataMC()) MCHistograms(); + // Correlate Calorimeters and V0 and track Multiplicity + if(fCorrelate) Correlate(); + //Get List with CaloClusters , calo Cells, init min amplitude TObjArray * caloClusters = NULL; AliVCaloCells * cells = 0x0; @@ -3363,10 +3401,7 @@ void AliAnaCalorimeterQA::MakeAnalysisFillHistograms() //printf("QA: N cells %d, N clusters %d \n",cells->GetNumberOfCells(),caloClusters->GetEntriesFast()); - // Correlate Calorimeters and V0 and track Multiplicity - if(fCorrelate) Correlate(); - - // Clusters + // Clusters ClusterLoopHistograms(caloClusters,cells); // Cells @@ -3382,124 +3417,133 @@ void AliAnaCalorimeterQA::MCHistograms() { //Get the MC arrays and do some checks before filling MC histograms + Int_t pdg = 0 ; + Int_t status = 0 ; + Int_t nprim = 0 ; + + TParticle * primStack = 0; + AliAODMCParticle * primAOD = 0; TLorentzVector mom ; - if(GetReader()->ReadStack()) - { - if(!GetMCStack()) - AliFatal("Stack not available, is the MC handler called?\n"); - - //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) ; - - 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 - } // ESD - else if(GetReader()->ReadAODMCParticles()) + // Get the ESD MC particles container + AliStack * stack = 0; + if( GetReader()->ReadStack() ) { - if(!GetReader()->GetAODMCParticles()) - AliFatal("AODMCParticles not available!"); - - //Fill some pure MC histograms, only primaries. - for(Int_t i=0 ; i < (GetReader()->GetAODMCParticles())->GetEntriesFast(); i++) - { - AliAODMCParticle *aodprimary = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(i) ; - - if (!aodprimary->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons - - mom.SetPxPyPzE(aodprimary->Px(),aodprimary->Py(),aodprimary->Pz(),aodprimary->E()); - MCHistograms(mom,TMath::Abs(aodprimary->GetPdgCode())); - } //primary loop - } // AOD -} - -//_______________________________________________________________________________ -void AliAnaCalorimeterQA::MCHistograms(TLorentzVector mom, Int_t pdg) -{ - //Fill pure monte carlo related histograms - - Float_t eMC = mom.E(); - Float_t phiMC = mom.Phi(); - if(phiMC < 0) - phiMC += TMath::TwoPi(); - Float_t etaMC = mom.Eta(); - - if (TMath::Abs(etaMC) > 1) return; + stack = GetMCStack(); + if(!stack ) return; + nprim = stack->GetNtrack(); + } - Bool_t in = kFALSE; + // Get the AOD MC particles container + TClonesArray * mcparticles = 0; + if( GetReader()->ReadAODMCParticles() ) + { + mcparticles = GetReader()->GetAODMCParticles(); + if( !mcparticles ) return; + nprim = mcparticles->GetEntriesFast(); + } - //Rough stimate of acceptance for pi0, Eta and electrons - if(fCalorimeter == "PHOS") + //printf("N primaries %d\n",nprim); + for(Int_t i=0 ; i < nprim; i++) { - if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter)) - in = kTRUE ; - if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in); + if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ; - } - else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()) - { - if(GetEMCALGeometry()) + // Get the generated particles, check that it is primary (not parton, resonance) + // and get its momentum. Different way to recover from ESD or AOD + if(GetReader()->ReadStack()) { - Int_t absID=0; - GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(mom.Eta(),mom.Phi(),absID); + primStack = stack->Particle(i) ; + if(!primStack) + { + printf("AliAnaCalorimeterQA::MCHistograms() - ESD primaries pointer not available!!\n"); + continue; + } + + pdg = primStack->GetPdgCode(); + status = primStack->GetStatusCode(); - if( absID >= 0) - in = kTRUE; + //printf("Input: i %d, %s, pdg %d, status %d \n",i, primStack->GetName(), pdg, status); - if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s Real acceptance? %d\n",fCalorimeter.Data(),in); + if ( status > 11 ) continue; //Working for PYTHIA and simple generators, check for HERWIG, HIJING? + + if ( primStack->Energy() == TMath::Abs(primStack->Pz()) ) continue ; //Protection against floating point exception + + //printf("Take : i %d, %s, pdg %d, status %d \n",i, primStack->GetName(), pdg, status); + + //Photon kinematics + primStack->Momentum(mom); } else { - if(GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter)) - in = kTRUE ; - if(GetDebug() > 2) printf("AliAnaCalorimeterQA::MCHistograms() - In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),in); - } - } - - if (pdg==22) - { - fhGenMCE[kmcPhoton] ->Fill(eMC); - if(eMC > 0.5) fhGenMCEtaPhi[kmcPhoton]->Fill(etaMC,phiMC); - if(in) - { - fhGenMCAccE[kmcPhoton] ->Fill(eMC); - if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPhoton]->Fill(etaMC,phiMC); + primAOD = (AliAODMCParticle *) mcparticles->At(i); + if(!primAOD) + { + printf("AliAnaCalorimeterQA::MCHistograms() - AOD primaries pointer not available!!\n"); + continue; + } + + pdg = primAOD->GetPdgCode(); + status = primAOD->GetStatus(); + + //printf("Input: i %d, %s, pdg %d, status %d \n",i, primAOD->GetName(), pdg, status); + + if (!primAOD->IsPrimary()) continue; //accept all which is not MC transport generated. Don't know how to avoid partons + + if ( status > 11 ) continue; //Working for PYTHIA and simple generators, check for HERWIG + + if ( primAOD->E() == TMath::Abs(primAOD->Pz()) ) continue ; //Protection against floating point exception + + //printf("Take : i %d, %s, pdg %d, status %d \n",i, primAOD->GetName(), pdg, status); + + //kinematics + mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E()); } - } - else if (pdg==111) - { - fhGenMCE[kmcPi0] ->Fill(eMC); - if(eMC > 0.5) fhGenMCEtaPhi[kmcPi0]->Fill(etaMC,phiMC); - if(in) + + Float_t eMC = mom.E(); + if(eMC < 0.2) continue; + Float_t ptMC = mom.E(); + + Float_t etaMC = mom.Eta(); + // Only particles in |eta| < 1 + if (TMath::Abs(etaMC) > 1) continue; + + Float_t phiMC = mom.Phi(); + if(phiMC < 0) + phiMC += TMath::TwoPi(); + + Int_t mcIndex = -1; + if (pdg==22) mcIndex = kmcPhoton; + else if (pdg==111) mcIndex = kmcPi0; + else if (pdg==221) mcIndex = kmcEta; + else if (TMath::Abs(pdg)==11) mcIndex = kmcElectron; + + if( mcIndex >=0 ) { - fhGenMCAccE[kmcPi0] ->Fill(eMC); - if(eMC > 0.5) fhGenMCAccEtaPhi[kmcPi0]->Fill(etaMC,phiMC); + fhGenMCE [mcIndex]->Fill( eMC); + fhGenMCPt[mcIndex]->Fill(ptMC); + if(eMC > 0.5) fhGenMCEtaPhi[mcIndex]->Fill(etaMC,phiMC); + + Bool_t inacceptance = kTRUE; + // Check same fidutial borders as in data analysis on top of real acceptance if real was requested. + if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ) inacceptance = kFALSE ; + + if(IsRealCaloAcceptanceOn()) // defined on base class + { + if(GetReader()->ReadStack() && + !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fCalorimeter, primStack)) inacceptance = kFALSE ; + if(GetReader()->ReadAODMCParticles() && + !GetCaloUtils()->IsMCParticleInCalorimeterAcceptance(fCalorimeter, primAOD )) inacceptance = kFALSE ; + } + + if(!inacceptance) continue; + + fhGenMCAccE [mcIndex]->Fill( eMC); + fhGenMCAccPt[mcIndex]->Fill(ptMC); + if(eMC > 0.5) fhGenMCAccEtaPhi[mcIndex]->Fill(etaMC,phiMC); + } + } - else if (pdg==221) - { - fhGenMCE[kmcEta] ->Fill(eMC); - if(eMC > 0.5) fhGenMCEtaPhi[kmcEta]->Fill(etaMC,phiMC); - if(in) - { - fhGenMCAccE[kmcEta] ->Fill(eMC); - if(eMC > 0.5) fhGenMCAccEtaPhi[kmcEta]->Fill(etaMC,phiMC); - } - } - else if (TMath::Abs(pdg)==11) - { - fhGenMCE[kmcElectron] ->Fill(eMC); - if(eMC > 0.5) fhGenMCEtaPhi[kmcElectron]->Fill(etaMC,phiMC); - if(in) - { - fhGenMCAccE[kmcElectron] ->Fill(eMC); - if(eMC > 0.5) fhGenMCAccEtaPhi[kmcElectron]->Fill(etaMC,phiMC); - } - } } //_________________________________________________________________________________ diff --git a/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.h b/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.h index de72e7920f1..d95b2f5f6dd 100755 --- a/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.h +++ b/PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.h @@ -54,7 +54,7 @@ public: void CellInClusterPositionHistograms(AliVCluster* cluster); - void ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax, const Bool_t goodCluster ); + void ClusterAsymmetryHistograms(AliVCluster* clus, Int_t absIdMax, Bool_t goodCluster ); void ClusterHistograms(AliVCluster* cluster, const TObjArray *caloClusters, AliVCaloCells * cells, Int_t absIdMax, Double_t maxCellFraction, Float_t eCrossFrac, Double_t tmax); @@ -81,8 +81,6 @@ public: void MCHistograms(); - void MCHistograms(const TLorentzVector mom, Int_t pdg); - void WeightHistograms(AliVCluster *clus, AliVCaloCells* cells); // Setters and Getters @@ -390,9 +388,11 @@ public: TH2F * fhRecoMCDeltaPhi[6][2]; //! Gen-Reco phi generated particle vs reconstructed E TH2F * fhRecoMCDeltaEta[6][2]; //! Gen-Reco eta generated particle vs reconstructed E - TH1F * fhGenMCE[4] ; //! pt of primary particle + TH1F * fhGenMCE [4] ; //! pt of primary particle + TH1F * fhGenMCPt[4] ; //! pt of primary particle TH2F * fhGenMCEtaPhi[4] ; //! eta vs phi of primary particle - TH1F * fhGenMCAccE[4] ; //! pt of primary particle, in acceptance + TH1F * fhGenMCAccE [4] ; //! pt of primary particle, in acceptance + TH1F * fhGenMCAccPt[4] ; //! pt of primary particle, in acceptance TH2F * fhGenMCAccEtaPhi[4] ; //! eta vs phi of primary particle, in acceptance TH2F * fhEMVxyz ; //! Electromagnetic particle production vertex @@ -439,7 +439,7 @@ public: AliAnaCalorimeterQA & operator = (const AliAnaCalorimeterQA & qa) ;//cpy assignment AliAnaCalorimeterQA( const AliAnaCalorimeterQA & qa) ; // cpy ctor - ClassDef(AliAnaCalorimeterQA,28) + ClassDef(AliAnaCalorimeterQA,29) } ; -- 2.43.0