merge the 2 MCHistogram methods, simplify the filling of pure MC histograms, introduc...
authorgconesab <gustavo.conesa.balbastre@cern.ch>
Thu, 7 Aug 2014 09:53:21 +0000 (11:53 +0200)
committergconesab <gustavo.conesa.balbastre@cern.ch>
Thu, 7 Aug 2014 13:06:16 +0000 (15:06 +0200)
PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.cxx
PWGGA/CaloTrackCorrelations/AliAnaCalorimeterQA.h

index a069c70..34af940 100755 (executable)
@@ -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 ; i<GetMCStack()->GetNprimary(); 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);  
-    }
-  }    
 }
 
 //_________________________________________________________________________________
index de72e79..d95b2f5 100755 (executable)
@@ -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)
 } ;