declare only once different TLorentzVectors and TVector3s
authorgconesab <gustavo.conesa.balbastre@cern.ch>
Sun, 19 Oct 2014 17:45:08 +0000 (19:45 +0200)
committergconesab <gustavo.conesa.balbastre@cern.ch>
Sun, 19 Oct 2014 20:47:35 +0000 (22:47 +0200)
14 files changed:
PWGGA/CaloTrackCorrelations/AliAnaElectron.cxx
PWGGA/CaloTrackCorrelations/AliAnaElectron.h
PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.h
PWGGA/CaloTrackCorrelations/AliAnaParticleJetFinderCorrelation.cxx
PWGGA/CaloTrackCorrelations/AliAnaParticleJetFinderCorrelation.h
PWGGA/CaloTrackCorrelations/AliAnaPhoton.cxx
PWGGA/CaloTrackCorrelations/AliAnaPhoton.h
PWGGA/CaloTrackCorrelations/AliAnaPhotonConvInCalo.cxx
PWGGA/CaloTrackCorrelations/AliAnaPhotonConvInCalo.h
PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx
PWGGA/CaloTrackCorrelations/AliAnaPi0.h
PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.h

index d8cfa9e..9e9bbd5 100755 (executable)
@@ -59,6 +59,7 @@ AliAnaElectron::AliAnaElectron() :
     fdEdxMin(0.),                         fdEdxMax (200.), 
     fEOverPMin(0),                        fEOverPMax (2),
     fAODParticle(0),
+    fMomentum(),                          fMomentumMC(),                         fProdVertex(),
     // Histograms
     fhdEdxvsE(0),                         fhdEdxvsP(0),                 
     fhEOverPvsE(0),                       fhEOverPvsP(0),
@@ -158,17 +159,17 @@ AliAnaElectron::AliAnaElectron() :
 }
 
 //____________________________________________________________________________
-Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
+Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, Int_t nMaxima)
 {
   //Select clusters if they pass different cuts
   if(GetDebug() > 2) 
     printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
            GetReader()->GetEventNumber(),
-           mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+           fMomentum.E(), fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta());
   
   //.......................................
   //If too small or big energy, skip it
-  if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ; 
+  if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) return kFALSE ; 
   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
   
   //.......................................
@@ -184,7 +185,7 @@ Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, I
   //.......................................
   //Check acceptance selection
   if(IsFiducialCutOn()){
-    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom.Eta(),mom.Phi(),GetCalorimeter()) ;
+    Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
     if(! in ) return kFALSE ;
   }
   if(GetDebug() > 2) printf("Fiducial cut passed \n");
@@ -215,7 +216,7 @@ Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, I
   if(GetDebug() > 0) 
     printf("AliAnaElectron::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
            GetReader()->GetEventNumber(), 
-           mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+           fMomentum.E(), fMomentum.Pt(),calo->E(),fMomentum.Phi()*TMath::RadToDeg(),fMomentum.Eta());
   
   //All checks passed, cluster selected
   return kTRUE;
@@ -243,19 +244,10 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTa
   
   Float_t l0   = 0., l1   = 0.;
   Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
-  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;    
+  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
   
-  TLorentzVector mom;
-  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
-  {
-    cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
-  else{
-    Double_t vertex[]={0,0,0};
-    cluster->GetMomentum(mom,vertex) ;
-  }
-  
-  Float_t eta = mom.Eta();
-  Float_t phi = mom.Phi();
+  Float_t eta = fMomentum.Eta();
+  Float_t phi = fMomentum.Phi();
   if(phi < 0) phi+=TMath::TwoPi();
   
   fhLam0E[pidIndex] ->Fill(energy,lambda0);
@@ -355,11 +347,12 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTa
         
         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
         Int_t ancPDG = 0, ancStatus = -1;
-        TLorentzVector momentum; TVector3 prodVertex;
         Int_t ancLabel = 0;
         Int_t noverlaps = 1;      
-        for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
-          ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
+        for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ )
+        {
+          ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),
+                                                               ancPDG,ancStatus,fMomentumMC,fProdVertex);
           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
         }
         
@@ -424,7 +417,7 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTa
       fhMCEDispPhi        [pidIndex][index]-> Fill(energy,dPhi);
       fhMCESumEtaPhi      [pidIndex][index]-> Fill(energy,sEtaPhi);
       fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
-      if(dEta+dPhi>0)fhMCESphericity     [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
+      if(dEta+dPhi>0) fhMCESphericity[pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
     }
     
   }//MC data
@@ -1119,7 +1112,6 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
   }
   
   //Init arrays, variables, get number of clusters
-  TLorentzVector mom, mom2 ;
   Int_t nCaloClusters = pl->GetEntriesFast();
   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
   
@@ -1146,10 +1138,10 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     //Cluster selection, not charged, with photon id and in fiducial cut         
     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
     {
-      calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
+      calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
     else{
       Double_t vertex[]={0,0,0};
-      calo->GetMomentum(mom,vertex) ;
+      calo->GetMomentum(fMomentum,vertex) ;
     }
     
     //--------------------------------------
@@ -1160,7 +1152,7 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     else                           cells = GetPHOSCells();
     
     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
-    if(!ClusterSelected(calo,mom,nMaxima)) continue;
+    if(!ClusterSelected(calo,nMaxima)) continue;
     
     //-------------------------------------
     // PID selection via dE/dx
@@ -1339,15 +1331,15 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     }
         
     if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
-                              mom.Pt(), pid);
+                              fMomentum.Pt(), pid);
     
     Float_t maxCellFraction = 0;
     Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
-    if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(mom.E(),maxCellFraction);
+    if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(fMomentum.E(),maxCellFraction);
     
-    fhNCellsE[pidIndex] ->Fill(mom.E(),calo->GetNCells());
-    fhNLME   [pidIndex] ->Fill(mom.E(),nMaxima);
-    fhTimeE  [pidIndex] ->Fill(mom.E(),calo->GetTOF()*1.e9);
+    fhNCellsE[pidIndex] ->Fill(fMomentum.E(),calo->GetNCells());
+    fhNLME   [pidIndex] ->Fill(fMomentum.E(),nMaxima);
+    fhTimeE  [pidIndex] ->Fill(fMomentum.E(),calo->GetTOF()*1.e9);
 
     
     //----------------------------
@@ -1357,7 +1349,7 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     //Add AOD with electron/hadron object to aod branch
     if ( pid == fAODParticle || fAODParticle == 0 ) 
     {
-      AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
+      AliAODPWG4Particle aodpart = AliAODPWG4Particle(fMomentum);
       
       //...............................................
       //Set the indeces of the original caloclusters (MC, ID), and calorimeter
index 3f007e2..5c33ae0 100755 (executable)
@@ -55,7 +55,7 @@ class AliAnaElectron : public AliAnaCaloTrackCorrBaseClass {
   
   // Analysis methods
   
-  Bool_t       ClusterSelected(AliVCluster* cl, TLorentzVector mom, Int_t nMaxima) ;
+  Bool_t       ClusterSelected(AliVCluster* cl, Int_t nMaxima) ;
   
   void         FillShowerShapeHistograms( AliVCluster* cluster, Int_t mcTag , Int_t pidTag) ;
   
@@ -138,6 +138,10 @@ class AliAnaElectron : public AliAnaCaloTrackCorrBaseClass {
 
   Int_t    fAODParticle;                       // Select the type of particle to put in AODs for other analysis
   
+  TLorentzVector fMomentum;                    //! cluster momentum
+  TLorentzVector fMomentumMC;                  //! mc particle momentum
+  TVector3       fProdVertex;                  //! mc particle production vertex
+
   //Histograms
   TH2F * fhdEdxvsE;                            //! matched track dEdx vs cluster E 
   TH2F * fhdEdxvsP;                            //! matched track dEdx vs track P
index 310abd4..795fcb5 100755 (executable)
@@ -63,6 +63,11 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
   fFillArmenterosHisto(0),                   fFillThetaStarHisto(0),
   fSSWeightN(0),                             fSSECellCutN(0),
   fNLMSettingN(0),                           fWSimu(),
+  fClusterMomentum(),                        fSubClusterMom1(),                         fSubClusterMom2(),
+  fSubClusterMomSum(),                       fSubClusterMomBoost(),
+  fPrimaryMom(),                             fMCDaughMom1(),                            fMCDaughMom2(),
+  fProdVertex(),
+// Histograms
   fhMassAsyCutNLocMax1(0),                   fhMassAsyCutNLocMax2(0),                   fhMassAsyCutNLocMaxN(0),
   fhM02AsyCutNLocMax1(0),                    fhM02AsyCutNLocMax2(0),                    fhM02AsyCutNLocMaxN(0),
   fhMassM02CutNLocMax1(0),                   fhMassM02CutNLocMax2(0),                   fhMassM02CutNLocMaxN(0),
@@ -582,10 +587,10 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
 //      Int_t   mpdg = -999999;
 //      Int_t   mstatus = -1;
 //      Int_t   grandLabel = -1;
-//      TLorentzVector mother = GetMCAnalysisUtils()->GetMother(mclabel,GetReader(),mpdg,mstatus,mOK,grandLabel);
+//      fPrimaryMom = GetMCAnalysisUtils()->GetMother(mclabel,GetReader(),mpdg,mstatus,mOK,grandLabel);
 //      
 //      printf("******** mother %d : Label %d, pdg %d; status %d, E %2.2f, Eta %2.2f, Phi %2.2f, ok %d, mother label %d\n",
-//             ilab, mclabel, mpdg, mstatus,mother.E(), mother.Eta(),mother.Phi()*TMath::RadToDeg(),mOK,grandLabel);
+//             ilab, mclabel, mpdg, mstatus,fPrimaryMom.E(), fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg(),mOK,grandLabel);
 //      
 //      if( ( mpdg == 22 || TMath::Abs(mpdg)==11 ) && grandLabel >=0 )
 //      {
@@ -709,7 +714,6 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
   
   // Compare the common ancestors of the 2 highest energy local maxima
   Int_t ancPDG = 0, ancStatus = -1;
-  TLorentzVector momentum; TVector3 prodVertex;
   Int_t ancLabel = 0;
   Bool_t high = kFALSE;
   Bool_t low  = kFALSE;
@@ -749,7 +753,7 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       }
       
       ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(mcLabel1,mcLabel2,
-                                                           GetReader(),ancPDG,ancStatus,momentum,prodVertex);
+                                                           GetReader(),ancPDG,ancStatus,fPrimaryMom,fProdVertex);
       if(ancPDG==111)
       {
         if((i==imax && j==imax2) ||  (j==imax && i==imax2))
@@ -770,7 +774,7 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
      
       Bool_t ok  =kFALSE;
       Int_t pdg = -22222, status = -1;
-      TLorentzVector primary  =GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
+      fPrimaryMom = GetMCAnalysisUtils()->GetMother(ancLabel,GetReader(), pdg, status, ok);
       //printf("\t i %d label %d - j %d label %d; ancestor label %d, PDG %d-%d; E %2.2f; high %d, any %d \n",i,mcLabel1,j,mcLabel2, ancLabel, ancPDG,pdg, primary.E(), high, low);
 
     }
@@ -829,11 +833,10 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
   Int_t gLabel = -1;
   
   Int_t label = cluster->GetLabel();
-  TLorentzVector pi0Kine;
-    
+  
   while( pdg!=111 && label>=0 )
   {
-    pi0Kine = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),pdg,status,ok, label,gLabel);
+    fPrimaryMom = GetMCAnalysisUtils()->GetGrandMother(label,GetReader(),pdg,status,ok, label,gLabel);
   }
   
   if(pdg!=111 || label < 0)
@@ -852,9 +855,9 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
   
   // Get daughter photon kinematics
   Int_t pdg0 = -22222, status0   = -1; Int_t label0 = -1;
-  TLorentzVector photon0Kine = GetMCAnalysisUtils()->GetDaughter(0,label,GetReader(),pdg0,status0,ok,label0);
+  fMCDaughMom1 = GetMCAnalysisUtils()->GetDaughter(0,label,GetReader(),pdg0,status0,ok,label0);
   Int_t pdg1 = -22222, status1   = -1; Int_t label1 = -1;
-  TLorentzVector photon1Kine = GetMCAnalysisUtils()->GetDaughter(1,label,GetReader(),pdg1,status1,ok,label1);
+  fMCDaughMom2 = GetMCAnalysisUtils()->GetDaughter(1,label,GetReader(),pdg1,status1,ok,label1);
 
   if(pdg1!=22 || pdg0 != 22)
   {
@@ -863,18 +866,18 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
   }
   
   // In what cells did the photons hit
-  Float_t eta0 = photon0Kine.Eta();
-  Float_t eta1 = photon1Kine.Eta();
+  Float_t eta0 = fMCDaughMom1.Eta();
+  Float_t eta1 = fMCDaughMom2.Eta();
   
-  Float_t phi0 = photon0Kine.Phi();
-  Float_t phi1 = photon1Kine.Phi();
+  Float_t phi0 = fMCDaughMom1.Phi();
+  Float_t phi1 = fMCDaughMom2.Phi();
 
 // if((mass < 0.06 || mass > 1.8) && mcindex==kmcPi0 && noverlaps == 0)
 //  {
 //    printf("MC pi0 label %d E  %2.2f, eta %2.2f, phi %2.2f, mass (ph1, ph2) %2.2f: \n \t photon0 label %d E %2.2f, eta %2.2f, phi %2.2f \n \t photon1 label %d E %2.2f eta %2.2f, phi %2.2f\n",
-//           label , pi0Kine.E()    , pi0Kine.Eta(),pi0Kine.Phi()*TMath::RadToDeg(), (photon0Kine+photon1Kine).M(),
-//           label0, photon0Kine.E(),          eta0,         phi0*TMath::RadToDeg(),
-//           label1, photon1Kine.E(),          eta1,         phi1*TMath::RadToDeg());
+//           label , fPrimaryMom.E()    , fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg(), (fMCDaughMom1+fMCDaughMom2).M(),
+//           label0, fMCDaughMom1.E(),          eta0,         phi0*TMath::RadToDeg(),
+//           label1, fMCDaughMom2.E(),          eta1,         phi1*TMath::RadToDeg());
 //    
 //    TLorentzVector momclus;
 //    cluster->GetMomentum(momclus,GetVertex(0));
@@ -978,7 +981,7 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       Int_t tmplabel   = mclabel;
       while((secLabel0 < 0 || secLabel1 < 0) && tmplabel > 0 )
       {
-        TLorentzVector mother = GetMCAnalysisUtils()->GetMother(tmplabel,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+        fPrimaryMom = GetMCAnalysisUtils()->GetMother(tmplabel,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
         
         //printf("\t \t while secLabel %d, mom %d, granmom %d\n",mclabel,tmplabel,secgrandLabel);
         
@@ -998,14 +1001,14 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
     // Get the position of the found secondaries mother
     if(!match0 && secLabel0 > 0)
     {
-      TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel0,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+      fPrimaryMom = GetMCAnalysisUtils()->GetMother(secLabel0,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
       
-      //Float_t eta = mother.Eta();
-      //Float_t phi = mother.Phi();
+      //Float_t eta = fPrimaryMom.Eta();
+      //Float_t phi = fPrimaryMom.Phi();
       //if(phi < 0 ) phi+=TMath::TwoPi();
       //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId0second);
       
-      //printf("Secondary MC0 label %d, absId %d E %2.2F eta %2.2f, phi %f\n", secLabel0,absId0second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
+      //printf("Secondary MC0 label %d, absId %d E %2.2F eta %2.2f, phi %f\n", secLabel0,absId0second, fPrimaryMom.E(),fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg());
       
       if(absId0second == list[imax] ) { match0 = kTRUE ; imatch0 = imax  ; }
       if(absId0second == list[imax2]) { match0 = kTRUE ; imatch0 = imax2 ; }
@@ -1013,14 +1016,14 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
 
     if(!match1 && secLabel1 > 0)
     {
-      TLorentzVector mother = GetMCAnalysisUtils()->GetMother(secLabel1,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
+      fPrimaryMom = GetMCAnalysisUtils()->GetMother(secLabel1,GetReader(),secpdg,secstatus,secOK,secgrandLabel);
       
-      //Float_t eta = mother.Eta();
-      //Float_t phi = mother.Phi();
+      //Float_t eta = fPrimaryMom.Eta();
+      //Float_t phi = fPrimaryMom.Phi();
       //if(phi < 0 ) phi+=TMath::TwoPi();
       //GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta, phi, absId1second);
       
-      //printf("Secondary MC1 label %d absId %d E %2.2F eta %2.2f, phi %f\n",secLabel1, absId1second, mother.E(),mother.Eta(),mother.Phi()*TMath::RadToDeg());
+      //printf("Secondary MC1 label %d absId %d E %2.2F eta %2.2f, phi %f\n",secLabel1, absId1second, fPrimaryMom.E(),fPrimaryMom.Eta(),fPrimaryMom.Phi()*TMath::RadToDeg());
       
       if(absId1second == list[imax] ) { match1 = kTRUE ; imatch1 = imax  ; }
       if(absId1second == list[imax2]) { match1 = kTRUE ; imatch1 = imax2 ; }
@@ -1042,17 +1045,17 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonHitHighLMMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     else
@@ -1061,17 +1064,17 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonHitHighLMOverlapMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax )
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
 
     }
@@ -1117,17 +1120,17 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
 
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     else
@@ -1136,17 +1139,17 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonAdjHighLMOverlapMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM1vsELM1[inlm]->Fill(e1,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjHighLMOverlapDiffELM2vsELM2[inlm]->Fill(e2,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     
@@ -1205,13 +1208,13 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitOtherLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     else
@@ -1220,13 +1223,13 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonHitOtherLMMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonHitOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     
@@ -1279,13 +1282,13 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonAdjOtherLMMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjOtherLMDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     else
@@ -1294,13 +1297,13 @@ void AliAnaInsideClusterInvariantMass::CheckLocalMaximaMCOrigin(AliVCluster* clu
       fhMCPi0DecayPhotonAdjOtherLMOverlapMass[inlm]->Fill(en,mass);
       if(match0 && imatch0 == imax)
       {
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon0Kine.E())/photon0Kine.E());
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon1Kine.E())/photon1Kine.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom1.E())/fMCDaughMom1.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom2.E())/fMCDaughMom2.E());
       }
       else
       {
-        if(photon1Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-photon1Kine.E())/photon1Kine.E());
-        if(photon0Kine.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-photon0Kine.E())/photon0Kine.E());
+        if(fMCDaughMom2.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM1[inlm]->Fill(en,(e1-fMCDaughMom2.E())/fMCDaughMom2.E());
+        if(fMCDaughMom1.E()>0)fhMCPi0DecayPhotonAdjOtherLMOverlapDiffELM2[inlm]->Fill(en,(e2-fMCDaughMom1.E())/fMCDaughMom1.E());
       }
     }
     
@@ -1496,24 +1499,23 @@ void AliAnaInsideClusterInvariantMass::FillAngleHistograms(Int_t   nMax,      Bo
 
 //______________________________________________________________________________________________________________________
 void AliAnaInsideClusterInvariantMass::FillArmenterosHistograms(Int_t nMax, Int_t ebin, Int_t mcIndex,
-                                                                Float_t en, TLorentzVector g1, TLorentzVector g2,
-                                                                Float_t m02, Int_t pid)
+                                                                Float_t en, Float_t m02, Int_t pid)
 {
   // Fill Armeteros type histograms
   
   // Get pTArm and AlphaArm
-  TLorentzVector pi0 = g1+g2;
-  Float_t momentumSquaredMother = pi0.P()*pi0.P();
+  fSubClusterMomSum = fSubClusterMom1+fSubClusterMom2;
+  Float_t momentumSquaredMother = fSubClusterMomSum.P()*fSubClusterMomSum.P();
   Float_t momentumDaughter1AlongMother = 0.;
   Float_t momentumDaughter2AlongMother = 0.;
 
   if (momentumSquaredMother > 0.)
   {
-    momentumDaughter1AlongMother = (g1.Px()*pi0.Px() + g1.Py()*pi0.Py()+ g1.Pz()*pi0.Pz()) / sqrt(momentumSquaredMother);
-    momentumDaughter2AlongMother = (g2.Px()*pi0.Px() + g2.Py()*pi0.Py()+ g2.Pz()*pi0.Pz()) / sqrt(momentumSquaredMother);
+    momentumDaughter1AlongMother = (fSubClusterMom1.Px()*fSubClusterMomSum.Px() + fSubClusterMom1.Py()*fSubClusterMomSum.Py()+ fSubClusterMom1.Pz()*fSubClusterMomSum.Pz()) / sqrt(momentumSquaredMother);
+    momentumDaughter2AlongMother = (fSubClusterMom2.Px()*fSubClusterMomSum.Px() + fSubClusterMom2.Py()*fSubClusterMomSum.Py()+ fSubClusterMom2.Pz()*fSubClusterMomSum.Pz()) / sqrt(momentumSquaredMother);
   }
 
-  Float_t momentumSquaredDaughter1 = g1.P()*g1.P();
+  Float_t momentumSquaredDaughter1 = fSubClusterMom1.P()*fSubClusterMom1.P();
   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
   
   Float_t pTArm = 0.;
@@ -1524,7 +1526,7 @@ void AliAnaInsideClusterInvariantMass::FillArmenterosHistograms(Int_t nMax, Int_
   if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
   
-  Float_t asym = TMath::Abs( g1.Energy()-g2.Energy() )/( g1.Energy()+g2.Energy() ) ;
+  Float_t asym = TMath::Abs( fSubClusterMom1.Energy()-fSubClusterMom2.Energy() )/( fSubClusterMom1.Energy()+fSubClusterMom2.Energy() ) ;
   
    if(GetDebug() > 2 ) Info("FillArmenterosHistograms()","E %f, alphaArm %f, pTArm %f\n",en,alphaArm,pTArm);
   
@@ -1537,7 +1539,7 @@ void AliAnaInsideClusterInvariantMass::FillArmenterosHistograms(Int_t nMax, Int_
   Int_t inlm = nMax-1;
   if(inlm > 2 ) inlm = 2;
   Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
-  if     (ensubcut > 0.1 && ensubcut < g1.E() && ensubcut < g2.E() ) eCutOK = kTRUE;
+  if     (ensubcut > 0.1 && ensubcut < fSubClusterMom1.E() && ensubcut < fSubClusterMom2.E() ) eCutOK = kTRUE;
   else if(ensubcut < 0.1)                                            eCutOK = kTRUE;
 
   
@@ -1599,19 +1601,18 @@ void AliAnaInsideClusterInvariantMass::FillArmenterosHistograms(Int_t nMax, Int_
 
 //______________________________________________________________________________________________________________
 void AliAnaInsideClusterInvariantMass::FillThetaStarHistograms(Int_t nMax, Bool_t matched, Int_t mcIndex,
-                                                               Float_t en, TLorentzVector g1, TLorentzVector g2,
-                                                               Float_t m02, Int_t pid)
+                                                               Float_t en, Float_t m02, Int_t pid)
 {
   // Fill cos Theta^star histograms
 
   
   // Get cos Theta^star
-  TLorentzVector pi0 = g1+g2;
-  TLorentzVector g1Boost = g1;
-  g1Boost.Boost(-pi0.BoostVector());
-  Float_t  cosThStar=TMath::Cos(g1Boost.Vect().Angle(pi0.Vect()));
+  fSubClusterMomSum = fSubClusterMom1+fSubClusterMom2;
+  fSubClusterMomBoost = fSubClusterMom1;
+  fSubClusterMomBoost.Boost(-fSubClusterMomSum.BoostVector());
+  Float_t  cosThStar=TMath::Cos(fSubClusterMomBoost.Vect().Angle(fSubClusterMomSum.Vect()));
   
-  Float_t asym = TMath::Abs( g1.Energy()-g2.Energy() )/( g1.Energy()+g2.Energy() ) ;
+  Float_t asym = TMath::Abs( fSubClusterMom1.Energy()-fSubClusterMom2.Energy() )/( fSubClusterMom1.Energy()+fSubClusterMom2.Energy() ) ;
 
   Bool_t m02OK = GetCaloPID()->IsInPi0M02Range(en,m02,nMax);
   Bool_t asyOK = GetCaloPID()->IsInPi0SplitAsymmetryRange(en,asym,nMax);
@@ -1622,7 +1623,7 @@ void AliAnaInsideClusterInvariantMass::FillThetaStarHistograms(Int_t nMax, Bool_
   Int_t inlm = nMax-1;
   if(inlm > 2 ) inlm = 2;
   Float_t ensubcut = GetCaloPID()->GetSubClusterEnergyMinimum(inlm);
-  if     (ensubcut > 0.1 && ensubcut < g1.E() && ensubcut < g2.E() ) eCutOK = kTRUE;
+  if     (ensubcut > 0.1 && ensubcut < fSubClusterMom1.E() && ensubcut < fSubClusterMom2.E() ) eCutOK = kTRUE;
   else if(ensubcut < 0.1)                                            eCutOK = kTRUE;
 
   //printf("Reco cos %f, asy %f\n",cosThStar,asym);
@@ -2613,7 +2614,6 @@ void AliAnaInsideClusterInvariantMass::FillNLMDiffCutHistograms(AliVCluster *clu
   
   Int_t    nlm  = 0;
   Double_t mass = 0., angle = 0.;
-  TLorentzVector    lv1, lv2;
   Int_t    absId1   =-1; Int_t   absId2   =-1;
   Float_t  distbad1 =-1; Float_t distbad2 =-1;
   Bool_t   fidcut1  = 0; Bool_t  fidcut2  = 0;
@@ -2633,7 +2633,7 @@ void AliAnaInsideClusterInvariantMass::FillNLMDiffCutHistograms(AliVCluster *clu
 
       pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(clus,cells,GetCaloUtils(),
                                                                                  GetVertex(0), nlm, mass, angle,
-                                                                                 lv1,lv2,absId1,absId2,
+                                                                                 fSubClusterMom1,fSubClusterMom2,absId1,absId2,
                                                                                  distbad1,distbad2,fidcut1,fidcut2);
       if (nlm <= 0)
       {
@@ -6244,8 +6244,8 @@ void AliAnaInsideClusterInvariantMass::GetMCPrimaryKine(AliVCluster* cluster, In
   Bool_t ok      = kFALSE;
   Int_t  mcLabel = cluster->GetLabel();
   
-  TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
-  eprim = primary.E();
+  fPrimaryMom = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
+  eprim = fPrimaryMom.E();
   
   Int_t mesonLabel = -1;
   
@@ -6255,15 +6255,15 @@ void AliAnaInsideClusterInvariantMass::GetMCPrimaryKine(AliVCluster* cluster, In
     {
       GetMCAnalysisUtils()->GetMCDecayAsymmetryAngleForPDG(mcLabel,111,GetReader(),asymGen,angleGen,ok);
       asymGen = TMath::Abs(asymGen);
-      TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
-      if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
+      fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
+      if(fGrandMotherMom.E() > 0 && ok) eprim =  fGrandMotherMom.E();
     }
     else
     {
       GetMCAnalysisUtils()->GetMCDecayAsymmetryAngleForPDG(mcLabel,221,GetReader(),asymGen,angleGen,ok);
       asymGen = TMath::Abs(asymGen);
-      TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
-      if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
+      fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
+      if(fGrandMotherMom.E() > 0 && ok) eprim =  fGrandMotherMom.E();
     }
   }
   
@@ -6407,10 +6407,9 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
 
     // Get cluster angles
     
-    TLorentzVector lv;
-    cluster->GetMomentum(lv, GetVertex(0));
-    Float_t eta = lv.Eta();
-    Float_t phi = lv.Phi();
+    cluster->GetMomentum(fClusterMomentum, GetVertex(0));
+    Float_t eta = fClusterMomentum.Eta();
+    Float_t phi = fClusterMomentum.Phi();
     if(phi<0) phi=+TMath::TwoPi();
     
     //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
@@ -6426,15 +6425,16 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
     
     Int_t    nMax = 0;
     Double_t mass = 0., angle = 0.;
-    TLorentzVector    lv1, lv2;
     Int_t    absId1   =-1; Int_t   absId2   =-1;
     Float_t  distbad1 =-1; Float_t distbad2 =-1;
     Bool_t   fidcut1  = 0; Bool_t  fidcut2  = 0;
     
     Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
                                                                                GetVertex(0), nMax, mass, angle,
-                                                                               lv1,lv2,absId1,absId2,
-                                                                               distbad1,distbad2,fidcut1,fidcut2);
+                                                                               fSubClusterMom1,fSubClusterMom2,
+                                                                               absId1,absId2,
+                                                                               distbad1,distbad2,
+                                                                               fidcut1,fidcut2);
     if (nMax <= 0) 
     {
       if(GetDebug() > 0 )
@@ -6476,8 +6476,8 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
 
     // Get sub-cluster parameters
     
-    Float_t e1 = lv1.Energy();
-    Float_t e2 = lv2.Energy();
+    Float_t e1 = fSubClusterMom1.Energy();
+    Float_t e2 = fSubClusterMom2.Energy();
     
     Double_t tof1  = cells->GetCellTime(absId1);
     GetCaloUtils()->RecalibrateCellTime(tof1, GetCalorimeter(), absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
@@ -6551,10 +6551,10 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
       FillAngleHistograms(nMax,matched,mcindex,en,e1,e2,angle,mass,angleGen,l0, asym,pidTag,noverlaps);
 
     if(fFillArmenterosHisto && ebin >= 0)
-      FillArmenterosHistograms(nMax, ebin, mcindex, en, lv1, lv2, l0, pidTag);
+      FillArmenterosHistograms(nMax, ebin, mcindex, en, l0, pidTag);
 
     if(fFillThetaStarHisto)
-      FillThetaStarHistograms(nMax,matched,mcindex, en, lv1, lv2, l0, pidTag);
+      FillThetaStarHistograms(nMax,matched,mcindex, en, l0, pidTag);
 
     
     //---------------------------------------------------------------------
index 6311f3c..5edb7e7 100755 (executable)
@@ -49,12 +49,10 @@ class AliAnaInsideClusterInvariantMass : public AliAnaCaloTrackCorrBaseClass {
                                    Float_t m02,   Float_t asym,    Int_t   pid,    Int_t   noverlaps);
   
   void         FillArmenterosHistograms(Int_t nMax, Int_t ebin, Int_t mcindex,
-                                        Float_t pi0E, TLorentzVector g1, TLorentzVector g2,
-                                        Float_t m02, Int_t pid);
+                                        Float_t pi0E, Float_t m02, Int_t pid);
 
   void         FillThetaStarHistograms(Int_t nMax, Bool_t matched, Int_t mcindex,
-                                       Float_t pi0E, TLorentzVector g1, TLorentzVector g2,
-                                       Float_t m02, Int_t pid);
+                                       Float_t pi0E, Float_t m02, Int_t pid);
 
   void         FillEBinHistograms(Int_t ebin, Int_t nMax, Int_t mcindex, Float_t splitFrac,
                                   Float_t mass, Float_t asym, Float_t l0);
@@ -235,6 +233,18 @@ class AliAnaInsideClusterInvariantMass : public AliAnaCaloTrackCorrBaseClass {
   Float_t      fWSimu[2];              // Constant and slope of the linear correction factor for the shower
                                        // shape weight in simulation, about 1-0.07*w
   
+  TLorentzVector fClusterMomentum;     //! Cluster momentum
+  TLorentzVector fSubClusterMom1;      //! Sub-Cluster momentum
+  TLorentzVector fSubClusterMom2;      //! Sub-Cluster momentum
+  TLorentzVector fSubClusterMomSum;    //! Sub-Cluster momentum sum, armenteros
+  TLorentzVector fSubClusterMomBoost;  //! Sub-Cluster momentum boosted, armenteros
+  
+  TLorentzVector fPrimaryMom;          //! Primary momentum
+  TLorentzVector fGrandMotherMom;      //! Primary momentum
+  TLorentzVector fMCDaughMom1;         //! Primary momentum
+  TLorentzVector fMCDaughMom2;         //! Primary momentum
+  TVector3       fProdVertex;          //! primary production vertex
+
   //Histograms
   
   TH2F       * fhMassNLocMax1[7][2]  ;                  //! Split Inv Mass vs cluster E, NLM=1,  different MC particle types, track matching on/off
index deced35..7531f12 100644 (file)
@@ -61,6 +61,7 @@ AliAnaCaloTrackCorrBaseClass(),
   fGammaConeSize(0.3),fUseBackgroundSubtractionGamma(kFALSE),fSaveGJTree(kTRUE),
   fMostEnergetic(kFALSE),fMostOpposite(kTRUE), fUseHistogramJetBkg(kTRUE),
   fUseHistogramTracks(kTRUE),fUseHistogramJetTracks(kTRUE),fMCStudies(kFALSE),fGenerator(0),
+  fMomentum(),
   fhDeltaEta(0), /*fhDeltaPhi(0),*/fhDeltaPhiCorrect(0),fhDeltaPhi0PiCorrect(0), fhDeltaPt(0), fhPtRatio(0), fhPt(0),
   fhFFz(0),fhFFxi(0),fhFFpt(0),fhNTracksInCone(0),
   fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetFFzCor(0),fhJetFFxiCor(0),
@@ -1715,25 +1716,24 @@ void  AliAnaParticleJetFinderCorrelation::MakeAnalysisFillHistograms()
       
       fGamSumPtNeu=0;
       fGamNclusters=0;
-      TLorentzVector mom ;
       //TVector3 p3Tmp;
       //Double_t scalarProduct=0;
       //Double_t vectorLength=particlecorr->P();
       for(Int_t icalo=0; icalo <clusters->GetEntriesFast(); icalo++){
         AliVCluster* calo = (AliVCluster *) clusters->At(icalo);
         if(clusterID==calo->GetID()) continue;//the same cluster as trigger
-        calo->GetMomentum(mom,vertex) ;//Assume that come from vertex in straight line
+        calo->GetMomentum(fMomentum,vertex) ;//Assume that come from vertex in straight line
         //printf("min pt %f\n",GetMinPt());
-        if(mom.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
-        p3Tmp.SetXYZ(mom.Px(),mom.Py(),mom.Pz());
+        if(fMomentum.Pt()<GetMinPt()) continue; //<<hardcoded here //FIXME 0.5 check if correct
+        p3Tmp.SetXYZ(fMomentum.Px(),fMomentum.Py(),fMomentum.Pz());
         //calculate sum pt in the cone
         if( TMath::Sqrt((particlecorr->Eta()-p3Tmp.Eta())*(particlecorr->Eta()-p3Tmp.Eta()) +
                         (particlecorr->Phi()-p3Tmp.Phi())*(particlecorr->Phi()-p3Tmp.Phi()) )<fGammaConeSize ){
-          //scalarProduct = particlecorr->Px()*mom.Px() + particlecorr->Py()*mom.Py() + particlecorr->Pz()*mom.Pz();
-          //scalarProduct/=mom.P();
+          //scalarProduct = particlecorr->Px()*fMomentum.Px() + particlecorr->Py()*fMomentum.Py() + particlecorr->Pz()*fMomentum.Pz();
+          //scalarProduct/=fMomentum.P();
           //scalarProduct/=vectorLength;
           //if(scalarProduct>TMath::Cos(0.3)) {//FIXME photon radius
-          fGamSumPtNeu+=mom.Pt();
+          fGamSumPtNeu+=fMomentum.Pt();
           fGamNclusters++;
         }
       }
@@ -2380,8 +2380,7 @@ void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
          fhMCPhotonEtaPhi->Fill(photonPhi,photonEta);
          
          //Check if photons hit the Calorimeter
-         TLorentzVector lv;
-         lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
+         fMomentum.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
          inacceptance = kFALSE;
          if(GetCaloUtils()->IsEMCALGeoMatrixSet()){
            fhMCPhotonCuts->Fill(4);
@@ -2392,7 +2391,7 @@ void AliAnaParticleJetFinderCorrelation::FindMCgenInfo(){
              if(GetDebug() > 3) printf("In EMCAL Real acceptance? %d\n",inacceptance);
            }
            else{
-             if(GetFiducialCut()->IsInFiducialCut(lv.Eta(),lv.Phi(),kEMCAL)) inacceptance = kTRUE ;
+             if(GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kEMCAL)) inacceptance = kTRUE ;
              if(GetDebug() > 3) printf("In EMCAL fiducial cut acceptance? %d\n",inacceptance);
            }
          }else{//no EMCAL nor EMCALGeoMatrixSet
index 943bf98..dccff8e 100644 (file)
@@ -156,6 +156,8 @@ private:
 
   TRandom2 * fGenerator;//! pointer to random generator object
 
+  TLorentzVector fMomentum; //! momentum
+  
   // Histograms
   TH2F *     fhDeltaEta;           //! Difference of jet eta and trigger particle eta as function of trigger particle pT
   //TH2F *     fhDeltaPhi;         //! Difference of jet phi and trigger particle phi as function of trigger particle pT
index f9176ed..5c1a09c 100755 (executable)
@@ -61,7 +61,7 @@ fNCellsCut(0),
 fNLMCutMin(-1),               fNLMCutMax(10),
 fFillSSHistograms(kFALSE),    fFillOnlySimpleSSHisto(1),
 fNOriginHistograms(8),        fNPrimaryHistograms(4),
-fMomentum(),
+fMomentum(),                  fPrimaryMom(),
 // Histograms
 
 // Control histograms
@@ -344,7 +344,6 @@ void AliAnaPhoton::FillAcceptanceHistograms()
   
   TParticle        * primStack = 0;
   AliAODMCParticle * primAOD   = 0;
-  TLorentzVector lv;
   
   // Get the ESD MC particles container
   AliStack * stack = 0;
@@ -386,7 +385,7 @@ void AliAnaPhoton::FillAcceptanceHistograms()
       //       prim->GetName(), prim->GetPdgCode());
       
       //Photon kinematics
-      primStack->Momentum(lv);
+      primStack->Momentum(fMomentum);
       
       photonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
     }
@@ -405,7 +404,7 @@ void AliAnaPhoton::FillAcceptanceHistograms()
       if(primAOD->E() == TMath::Abs(primAOD->Pz()))  continue ; //Protection against floating point exception
       
       //Photon kinematics
-      lv.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
+      fMomentum.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
 
       photonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
     }
@@ -414,13 +413,13 @@ void AliAnaPhoton::FillAcceptanceHistograms()
     if(pdg != 22 ) continue ;
     
     // If too small or too large pt, skip, same cut as for data analysis
-    photonPt  = lv.Pt () ;
+    photonPt  = fMomentum.Pt () ;
     
     if(photonPt < GetMinPt() || photonPt > GetMaxPt() ) continue ;
     
-    photonE   = lv.E  () ;
-    photonEta = lv.Eta() ;
-    photonPhi = lv.Phi() ;
+    photonE   = fMomentum.E  () ;
+    photonEta = fMomentum.Eta() ;
+    photonPhi = fMomentum.Phi() ;
     
     if(photonPhi < 0) photonPhi+=TMath::TwoPi();
     
@@ -428,7 +427,7 @@ void AliAnaPhoton::FillAcceptanceHistograms()
     inacceptance = kTRUE;
     
     // Check same fidutial borders as in data analysis on top of real acceptance if real was requested.
-    if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(lv.Eta(),lv.Phi(),GetCalorimeter())) inacceptance = kFALSE ;
+    if( IsFiducialCutOn() && !GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter())) inacceptance = kFALSE ;
     
     // Check if photons hit the Calorimeter acceptance
     if(IsRealCaloAcceptanceOn()) // defined on base class
@@ -2567,11 +2566,11 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
       Float_t eprim   = 0;
       Float_t ptprim  = 0;
       Bool_t ok = kFALSE;
-      TLorentzVector primary = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
+      fPrimaryMom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
       if(ok)
       {
-        eprim   = primary.Energy();
-        ptprim  = primary.Pt();
+        eprim   = fPrimaryMom.Energy();
+        ptprim  = fPrimaryMom.Pt();
       }
       
       Int_t tag =ph->GetTag();
index f3745a0..cd81d8d 100755 (executable)
@@ -129,6 +129,7 @@ class AliAnaPhoton : public AliAnaCaloTrackCorrBaseClass {
   Int_t    fNPrimaryHistograms;                     // Fill only NPrimaryHistograms of the 7 defined types
   
   TLorentzVector fMomentum;                         //! Cluster momentum
+  TLorentzVector fPrimaryMom;                       //! Primary MC momentum
   
   //Histograms 
   TH1F * fhClusterCutsE [10];                       //! control histogram on the different photon selection cuts, E
index fbaeadd..bb78fd1 100755 (executable)
@@ -53,7 +53,7 @@ fAddConvertedPairsToAOD(kFALSE),
 fMassCut(0),                  
 fConvAsymCut(1.),                  fConvDEtaCut(2.),
 fConvDPhiMinCut(-1.),              fConvDPhiMaxCut(7.), 
-
+fMomentum(),                       fProdVertex(),
 // Histograms
 fhPtPhotonConv(0),                 fhEtaPhiPhotonConv(0),          fhEtaPhi05PhotonConv(0),
 fhConvDeltaEta(0),                 fhConvDeltaPhi(0),              fhConvDeltaEtaPhi(0), 
@@ -584,10 +584,8 @@ void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
           //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
           Int_t ancPDG    = 0;
           Int_t ancStatus = 0;
-          TLorentzVector momentum;
-          TVector3 prodVertex;
           Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(cluster1->GetLabel(), cluster2->GetLabel(), 
-                                                                      GetReader(), ancPDG, ancStatus, momentum, prodVertex);
+                                                                      GetReader(), ancPDG, ancStatus, fMomentum, fProdVertex);
           
           // printf("AliAnaPhotonConvInCalo::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
           //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
@@ -603,12 +601,12 @@ void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
               fhConvPtMCConversion         ->Fill( pairM, calo->Pt()+calo2->Pt());
               fhConvDispersionMCConversion ->Fill( cluster1->GetDispersion(), cluster2->GetDispersion());
               fhConvM02MCConversion        ->Fill( cluster1->GetM02(), cluster2->GetM02());
-              fhConvDistMCConversion       ->Fill( convDist , prodVertex.Mag() );
-              fhConvDistMCConversion       ->Fill( convDist2, prodVertex.Mag() );
+              fhConvDistMCConversion       ->Fill( convDist , fProdVertex.Mag() );
+              fhConvDistMCConversion       ->Fill( convDist2, fProdVertex.Mag() );
               
               if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
-                fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
-                fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
+                fhConvDistMCConversionCuts->Fill( convDist , fProdVertex.Mag() );
+                fhConvDistMCConversionCuts->Fill( convDist2, fProdVertex.Mag() );
               }
               
             }              
@@ -661,8 +659,8 @@ void  AliAnaPhotonConvInCalo::MakeAnalysisFillAOD()
       //Add to AOD
       if(fAddConvertedPairsToAOD){
         //Create AOD of pair analysis
-        TLorentzVector mpair = *(calo->Momentum())+*(calo2->Momentum());
-        AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
+        fMomentum = *(calo->Momentum())+*(calo2->Momentum());
+        AliAODPWG4Particle aodpair = AliAODPWG4Particle(fMomentum);
         aodpair.SetLabel(calo->GetLabel());
         
         //printf("Index %d, Id %d\n",iaod, calo->GetID());
index 2f9de34..b95d756 100755 (executable)
@@ -80,6 +80,9 @@ class AliAnaPhotonConvInCalo : public AliAnaCaloTrackCorrBaseClass {
   Float_t  fConvDPhiMinCut;               // Select conversion pairs when dphi of pair lager than cut
   Float_t  fConvDPhiMaxCut;               // Select conversion pairs when dphi of pair smaller than cut
 
+  TLorentzVector fMomentum ;              //! Cluster momentum
+  TVector3       fProdVertex;             //! Production vertex
+  
   // Histograms
   TH1F * fhPtPhotonConv   ;               //! Number of identified photon vs transerse momentum 
   TH2F * fhEtaPhiPhotonConv  ;            //! Pseudorapidity vs Phi of identified  photon for transerse momentum > 0.5, for converted
index cfc3288..39566aa 100755 (executable)
@@ -69,6 +69,9 @@ fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
 fFillBadDistHisto(kFALSE),   fFillSSCombinations(kFALSE),  
 fFillAngleHisto(kFALSE),     fFillAsymmetryHisto(kFALSE),  fFillOriginHisto(0),          fFillArmenterosThetaStar(0),
 fCheckAccInSector(kFALSE),
+fPhotonMom1(),               fPhotonMom1Boost(),           fPhotonMom2(),                fPi0Mom(),
+fProdVertex(),
+
 //Histograms
 fhAverTotECluster(0),        fhAverTotECell(0),            fhAverTotECellvsCluster(0),
 fhEDensityCluster(0),        fhEDensityCell(0),            fhEDensityCellvsCluster(0),
@@ -1262,7 +1265,6 @@ void AliAnaPi0::FillAcceptanceHistograms()
   
   TParticle        * primStack = 0;
   AliAODMCParticle * primAOD   = 0;
-  TLorentzVector lvmeson;
   
   // Get the ESD MC particles container
   AliStack * stack = 0;
@@ -1308,7 +1310,7 @@ void AliAnaPi0::FillAcceptanceHistograms()
       //       prim->GetName(), prim->GetPdgCode());
       
       //Photon kinematics
-      primStack->Momentum(lvmeson);
+      primStack->Momentum(fPi0Mom);
       
       mesonY = 0.5*TMath::Log((primStack->Energy()+primStack->Pz())/(primStack->Energy()-primStack->Pz())) ;
     }
@@ -1332,7 +1334,7 @@ void AliAnaPi0::FillAcceptanceHistograms()
       if(primAOD->E() == TMath::Abs(primAOD->Pz()))  continue ; //Protection against floating point exception
       
       //Photon kinematics
-      lvmeson.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
+      fPi0Mom.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
       
       mesonY = 0.5*TMath::Log((primAOD->E()+primAOD->Pz())/(primAOD->E()-primAOD->Pz())) ;
     }
@@ -1340,10 +1342,10 @@ void AliAnaPi0::FillAcceptanceHistograms()
     // Select only pi0 or eta
     if( pdg != 111 && pdg != 221) continue ;
     
-    mesonPt  = lvmeson.Pt () ;
-    mesonE   = lvmeson.E  () ;
-    mesonYeta= lvmeson.Eta() ;
-    mesonPhi = lvmeson.Phi() ;
+    mesonPt  = fPi0Mom.Pt () ;
+    mesonE   = fPi0Mom.E  () ;
+    mesonYeta= fPi0Mom.Eta() ;
+    mesonPhi = fPi0Mom.Phi() ;
     if( mesonPhi < 0 ) mesonPhi+=TMath::TwoPi();
     mesonPhi *= TMath::RadToDeg();
     
@@ -1449,7 +1451,6 @@ void AliAnaPi0::FillAcceptanceHistograms()
     
     if(iphot1 < 0 || iphot1 >= nprim || iphot2 < 0 || iphot2 >= nprim) continue ;
     
-    TLorentzVector lv1, lv2;
     Int_t pdg1 = 0;
     Int_t pdg2 = 0;
     Bool_t inacceptance1 = kTRUE;
@@ -1465,8 +1466,8 @@ void AliAnaPi0::FillAcceptanceHistograms()
       pdg1 = phot1->GetPdgCode();
       pdg2 = phot2->GetPdgCode();
       
-      phot1->Momentum(lv1);
-      phot2->Momentum(lv2);
+      phot1->Momentum(fPhotonMom1);
+      phot2->Momentum(fPhotonMom2);
       
       // Check if photons hit the Calorimeter acceptance
       if(IsRealCaloAcceptanceOn())
@@ -1486,8 +1487,8 @@ void AliAnaPi0::FillAcceptanceHistograms()
       pdg1 = phot1->GetPdgCode();
       pdg2 = phot2->GetPdgCode();
       
-      lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
-      lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
+      fPhotonMom1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
+      fPhotonMom2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
       
       // Check if photons hit the Calorimeter acceptance
       if(IsRealCaloAcceptanceOn())
@@ -1502,25 +1503,25 @@ void AliAnaPi0::FillAcceptanceHistograms()
     // Check if photons hit desired acceptance in the fidutial borders fixed in the analysis
     if(IsFiducialCutOn())
     {
-      if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(lv1.Eta(), lv1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
-      if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(lv2.Eta(), lv2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
+      if( inacceptance1 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom1.Eta(), fPhotonMom1.Phi(), GetCalorimeter()) ) inacceptance1 = kFALSE ;
+      if( inacceptance2 && !GetFiducialCut()->IsInFiducialCut(fPhotonMom2.Eta(), fPhotonMom2.Phi(), GetCalorimeter()) ) inacceptance2 = kFALSE ;
     }
     
-    if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
+    if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg);
 
     if(GetCalorimeter()==kEMCAL && inacceptance1 && inacceptance2 && fCheckAccInSector)
     {
       Int_t absID1=0;
       Int_t absID2=0;
       
-      Float_t photonPhi1 = lv1.Phi();
-      Float_t photonPhi2 = lv2.Phi();
+      Float_t photonPhi1 = fPhotonMom1.Phi();
+      Float_t photonPhi2 = fPhotonMom2.Phi();
       
       if(photonPhi1 < 0) photonPhi1+=TMath::TwoPi();
       if(photonPhi2 < 0) photonPhi2+=TMath::TwoPi();
       
-      GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv1.Eta(),photonPhi1,absID1);
-      GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(lv2.Eta(),photonPhi2,absID2);
+      GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom1.Eta(),photonPhi1,absID1);
+      GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(fPhotonMom2.Eta(),photonPhi2,absID2);
       
       Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
       Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
@@ -1547,15 +1548,15 @@ void AliAnaPi0::FillAcceptanceHistograms()
       printf("Accepted in %s?: m (%2.2f,%2.2f,%2.2f), p1 (%2.2f,%2.2f,%2.2f), p2 (%2.2f,%2.2f,%2.2f) : in1 %d, in2 %d\n",
              GetCalorimeterString().Data(),
              mesonPt,mesonYeta,mesonPhi,
-             lv1.Pt(),lv1.Eta(),lv1.Phi()*TMath::RadToDeg(),
-             lv2.Pt(),lv2.Eta(),lv2.Phi()*TMath::RadToDeg(),
+             fPhotonMom1.Pt(),fPhotonMom1.Eta(),fPhotonMom1.Phi()*TMath::RadToDeg(),
+             fPhotonMom2.Pt(),fPhotonMom2.Eta(),fPhotonMom2.Phi()*TMath::RadToDeg(),
              inacceptance1, inacceptance2);
 
     
     if(inacceptance1 && inacceptance2)
     {
-      Float_t  asym  = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
-      Double_t angle = lv1.Angle(lv2.Vect());
+      Float_t  asym  = TMath::Abs((fPhotonMom1.E()-fPhotonMom2.E()) / (fPhotonMom1.E()+fPhotonMom2.E()));
+      Double_t angle = fPhotonMom1.Angle(fPhotonMom2.Vect());
       
       if(GetDebug() > 2)
         printf("\t ACCEPTED pdg %d: pt %2.2f, phi %2.2f, eta %2.2f\n",pdg,mesonPt,mesonPhi,mesonYeta);
@@ -1600,24 +1601,23 @@ void AliAnaPi0::FillAcceptanceHistograms()
   
 }
 
-//__________________________________________________________________________________
-void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg,             TLorentzVector meson,
-                                        TLorentzVector daugh1, TLorentzVector daugh2)
+//________________________________________________
+void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg)
 {
   // Fill armenteros plots
   
   // Get pTArm and AlphaArm
-  Float_t momentumSquaredMother = meson.P()*meson.P();
+  Float_t momentumSquaredMother = fPi0Mom.P()*fPi0Mom.P();
   Float_t momentumDaughter1AlongMother = 0.;
   Float_t momentumDaughter2AlongMother = 0.;
   
   if (momentumSquaredMother > 0.)
   {
-    momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
-    momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
+    momentumDaughter1AlongMother = (fPhotonMom1.Px()*fPi0Mom.Px() + fPhotonMom1.Py()*fPi0Mom.Py()+ fPhotonMom1.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
+    momentumDaughter2AlongMother = (fPhotonMom2.Px()*fPi0Mom.Px() + fPhotonMom2.Py()*fPi0Mom.Py()+ fPhotonMom2.Pz()*fPi0Mom.Pz()) / sqrt(momentumSquaredMother);
   }
   
-  Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
+  Float_t momentumSquaredDaughter1 = fPhotonMom1.P()*fPhotonMom1.P();
   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
   
   Float_t pTArm = 0.;
@@ -1625,14 +1625,14 @@ void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg,             TLorentzVector me
     pTArm = sqrt(ptArmSquared);
   
   Float_t alphaArm = 0.;
-  if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
+  if(momentumDaughter1AlongMother + momentumDaughter2AlongMother > 0)
     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
   
-  TLorentzVector daugh1Boost = daugh1;
-  daugh1Boost.Boost(-meson.BoostVector());
-  Float_t  cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
+  fPhotonMom1Boost = fPhotonMom1;
+  fPhotonMom1Boost.Boost(-fPi0Mom.BoostVector());
+  Float_t  cosThStar=TMath::Cos(fPhotonMom1Boost.Vect().Angle(fPi0Mom.Vect()));
   
-  Float_t en   = meson.Energy();
+  Float_t en   = fPi0Mom.Energy();
   Int_t   ebin = -1;
   if(en > 8  && en <= 12) ebin = 0;
   if(en > 12 && en <= 16) ebin = 1;
@@ -1655,7 +1655,7 @@ void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg,             TLorentzVector me
   if(GetDebug() > 2 )
   {
     Float_t asym = 0;
-    if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
+    if(fPhotonMom1.E()+fPhotonMom2.E() > 0) asym = TMath::Abs(fPhotonMom1.E()-fPhotonMom2.E())/(fPhotonMom1.E()+fPhotonMom2.E());
 
     printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
          en,alphaArm,pTArm,cosThStar,asym);
@@ -1676,10 +1676,8 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
   
   Int_t ancPDG    = 0;
   Int_t ancStatus = 0;
-  TLorentzVector ancMomentum;
-  TVector3 prodVertex;
   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
-                                                              GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
+                                                              GetReader(), ancPDG, ancStatus,fPi0Mom, fProdVertex);
   
   Int_t momindex  = -1;
   Int_t mompdg    = -1;
@@ -1721,8 +1719,8 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
               {
                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
-                fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
-                if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
+                fhMCPi0MassPtTrue[index]->Fill(fPi0Mom.Pt(),mass);
+                if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt);
               }//pass the different cuts
             }// pid bit cut loop
           }// icell loop
@@ -1730,11 +1728,11 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
       }//Multi cut ana sim
       else
       {
-        fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
+        fhMCPi0MassPtTrue[0]->Fill(fPi0Mom.Pt(),mass);
         
         if(mass < 0.17 && mass > 0.1)
         {
-          fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
+          fhMCPi0PtTruePtRec[0]->Fill(fPi0Mom.Pt(),pt);
           
           //Int_t uniqueId = -1;
           if(GetReader()->ReadStack())
@@ -1799,8 +1797,8 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
               {
                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
-                fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
-                if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
+                fhMCEtaMassPtTrue[index]->Fill(fPi0Mom.Pt(),mass);
+                if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(fPi0Mom.Pt(),pt);
               }//pass the different cuts
             }// pid bit cut loop
           }// icell loop
@@ -1808,8 +1806,8 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
       } //Multi cut ana sim
       else
       {
-        fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
-        if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
+        fhMCEtaMassPtTrue[0]->Fill(fPi0Mom.Pt(),mass);
+        if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(fPi0Mom.Pt(),pt);
         
         if(GetReader()->ReadStack())
         {
@@ -2004,7 +2002,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
     
     //Get the momentum of this cluster
-    TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
+    fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
     
     //Get (Super)Module number of this cluster
     module1 = GetModuleNumber(p1);
@@ -2076,31 +2074,32 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
       
       //Get the momentum of this cluster
-      TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
+      fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
+      
       //Get module number
       module2       = GetModuleNumber(p2);
       
       //---------------------------------
       // Get pair kinematics
       //---------------------------------
-      Double_t m    = (photon1 + photon2).M() ;
-      Double_t pt   = (photon1 + photon2).Pt();
-      Double_t deta = photon1.Eta() - photon2.Eta();
-      Double_t dphi = photon1.Phi() - photon2.Phi();
+      Double_t m    = (fPhotonMom1 + fPhotonMom2).M() ;
+      Double_t pt   = (fPhotonMom1 + fPhotonMom2).Pt();
+      Double_t deta = fPhotonMom1.Eta() - fPhotonMom2.Eta();
+      Double_t dphi = fPhotonMom1.Phi() - fPhotonMom2.Phi();
       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
       
       if(GetDebug() > 2)
-        printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
+        printf(" E: fPhotonMom1 %f, fPhotonMom2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (fPhotonMom1 + fPhotonMom2).E(),m,a);
       
       //--------------------------------
       // Opening angle selection
       //--------------------------------
       //Check if opening angle is too large or too small compared to what is expected  
-      Double_t angle   = photon1.Angle(photon2.Vect());
-      if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05))
+      Double_t angle   = fPhotonMom1.Angle(fPhotonMom2.Vect());
+      if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
       {
         if(GetDebug() > 2)
-          printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
+          printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (fPhotonMom1+fPhotonMom2).E());
         continue;
       }
       
@@ -2328,7 +2327,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
         
         //Get kinematics of cluster and (super) module of this cluster
-        TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
+        fPhotonMom1.SetPxPyPzE(p1->Px(),p1->Py(),p1->Pz(),p1->E());
         module1 = GetModuleNumber(p1);
         
         //---------------------------------
@@ -2339,17 +2338,17 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
           
           //Get kinematics of second cluster and calculate those of the pair
-          TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
-          m           = (photon1+photon2).M() ; 
-          Double_t pt = (photon1 + photon2).Pt();
+          fPhotonMom2.SetPxPyPzE(p2->Px(),p2->Py(),p2->Pz(),p2->E());
+          m           = (fPhotonMom1+fPhotonMom2).M() ; 
+          Double_t pt = (fPhotonMom1 + fPhotonMom2).Pt();
           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
           
           //Check if opening angle is too large or too small compared to what is expected
-          Double_t angle   = photon1.Angle(photon2.Vect());
-          if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05))
+          Double_t angle   = fPhotonMom1.Angle(fPhotonMom2.Vect());
+          if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((fPhotonMom1+fPhotonMom2).E(),angle+0.05))
           {
             if(GetDebug() > 2)
-              printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
+              printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (fPhotonMom1+fPhotonMom2).E());
             continue;
           }
           
@@ -2361,7 +2360,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
           } 
           
           if(GetDebug() > 2)
-            printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
+            printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: fPhotonMom1 %2.2f, fPhotonMom2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
                    p1->Pt(), p2->Pt(), pt,m,a);        
           
           //In case we want only pairs in same (super) module, check their origin.
index 30360d9..9e17c54 100755 (executable)
@@ -130,8 +130,7 @@ class AliAnaPi0 : public AliAnaCaloTrackCorrBaseClass {
                                              Double_t mass,    Double_t pt,     Double_t asym,
                                              Double_t deta,    Double_t dphi);
   
-  void         FillArmenterosThetaStar(Int_t pdg,             TLorentzVector meson,
-                                       TLorentzVector daugh1, TLorentzVector daugh2);
+  void         FillArmenterosThetaStar(Int_t pdg);
 
   
   private:
@@ -171,6 +170,12 @@ class AliAnaPi0 : public AliAnaCaloTrackCorrBaseClass {
   
   Bool_t   fCheckAccInSector;          // Check that the decay pi0 falls in the same SM or sector
   
+  TLorentzVector fPhotonMom1;          //! photon cluster momentum
+  TLorentzVector fPhotonMom1Boost;     //! photon cluster momentum
+  TLorentzVector fPhotonMom2;          //! photon cluster momentum
+  TLorentzVector fPi0Mom;              //! pi0 cluster momentum
+  TVector3       fProdVertex;          //! production vertex
+  
   //Histograms
   
   //Event characterization
@@ -320,7 +325,7 @@ class AliAnaPi0 : public AliAnaCaloTrackCorrBaseClass {
   AliAnaPi0(              const AliAnaPi0 & api0) ; // cpy ctor
   AliAnaPi0 & operator = (const AliAnaPi0 & api0) ; // cpy assignment
   
-  ClassDef(AliAnaPi0,28)
+  ClassDef(AliAnaPi0,29)
 } ;
 
 
index 4cbe777..4a5401f 100755 (executable)
@@ -60,6 +60,8 @@ fFillEMCALBCHistograms(0),
 fFillAllNLMHistograms(0),
 fInputAODGammaConvName(""),
 fCheckSplitDistToBad(0),
+fMomentum(),  fMomentum1(),  fMomentum2(),
+fMomentum12(),fPrimaryMom(), fGrandMotherMom(),
 // Histograms
 fhPt(0),                            fhE(0),
 fhPtEta(0),                         fhPtPhi(0),                         fhEtaPhi(0),
@@ -87,6 +89,7 @@ fhSplitE(0),                        fhSplitPt(0),
 fhSplitPtEta(0),                    fhSplitPtPhi(0),
 fhNLocMaxSplitPt(0),
 fhPtDecay(0),
+
 // Shower shape histos
 fhPtDispersion(0),                  fhPtLambda0(0),                     fhPtLambda0NoSplitCut(0),
 fhPtLambda1(0),                     fhPtLambda0NoTRD(0),                fhPtLambda0FracMaxCellCut(0),
@@ -482,15 +485,15 @@ void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster *
 
 
 //______________________________________________________________________________________________
-void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
+void AliAnaPi0EbE::FillRejectedClusterHistograms(Int_t mctag, Int_t nMaxima)
 {
   // Fill histograms that do not pass the identification (SS case only)
   
-  Float_t ener  = mom.E();
-  Float_t pt    = mom.Pt();
-  Float_t phi   = mom.Phi();
+  Float_t ener  = fMomentum.E();
+  Float_t pt    = fMomentum.Pt();
+  Float_t phi   = fMomentum.Phi();
   if(phi < 0) phi+=TMath::TwoPi();
-  Float_t eta = mom.Eta();
+  Float_t eta = fMomentum.Eta();
 
   fhPtReject     ->Fill(pt);
   fhEReject      ->Fill(ener);
@@ -2616,22 +2619,16 @@ Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
 }
 
 //__________________________________________________________________
-void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
-                                       AliAODPWG4Particle * photon2,
+void AliAnaPi0EbE::HasPairSameMCMother(Int_t label1 , Int_t label2,
+                                       Int_t tag1   , Int_t tag2,
                                        Int_t & label, Int_t & tag)
 {
   // Check the labels of pair in case mother was same pi0 or eta
   // Set the new AOD accordingly
   
-  Int_t  label1 = photon1->GetLabel();
-  Int_t  label2 = photon2->GetLabel();
   
   if(label1 < 0 || label2 < 0 || label1 == label2) return ;
   
-  //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
-  //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
-  Int_t tag1 = photon1->GetTag();
-  Int_t tag2 = photon2->GetTag();
   
   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
@@ -2687,16 +2684,13 @@ void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
     {
       label = label1;
       
-      TLorentzVector mom1 = *(photon1->Momentum());
-      TLorentzVector mom2 = *(photon2->Momentum());
-      
-      Double_t angle = mom2.Angle(mom1.Vect());
-      Double_t mass  = (mom1+mom2).M();
-      Double_t epair = (mom1+mom2).E();
+      Double_t angle = fMomentum2.Angle(fMomentum1.Vect());
+      Double_t mass  = (fMomentum1+fMomentum2).M();
+      Double_t epair = (fMomentum1+fMomentum2).E();
       
       if(pdg1==111)
       {
-        //printf("Real pi0 pair: pt %f, mass %f\n",(mom1+mom2).Pt(),mass);
+        //printf("Real pi0 pair: pt %f, mass %f\n",(fMomentum1+fMomentum2).Pt(),mass);
         fhMassPairMCPi0 ->Fill(epair,mass);
         fhAnglePairMCPi0->Fill(epair,angle);
         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
@@ -2788,10 +2782,6 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
   //Read photon list from AOD, produced in class AliAnaPhoton
   //Check if 2 photons have the mass of the pi0.
   
-  TLorentzVector mom1;
-  TLorentzVector mom2;
-  TLorentzVector mom ;
-  
   if(!GetInputAODBranch())
   {
     AliFatal(Form("No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data()));
@@ -2816,7 +2806,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
     }
     
-    mom1 = *(photon1->Momentum());
+    fMomentum1 = *(photon1->Momentum());
     
     //Get original cluster, to recover some information
     Int_t iclus = -1;
@@ -2853,7 +2843,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
       }
       
-      mom2 = *(photon2->Momentum());
+      fMomentum2 = *(photon2->Momentum());
       
       //Get original cluster, to recover some information
       Int_t iclus2 = -1;
@@ -2882,7 +2872,9 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       Int_t label   =-1;
       if(IsDataMC())
       {
-        HasPairSameMCMother(photon1, photon2, label, tag) ;
+        HasPairSameMCMother(photon1->GetLabel(), photon2->GetLabel(),
+                            photon1->GetTag()  , photon2->GetTag(),
+                            label, tag) ;
         mcIndex = GetMCIndex(tag);
       }
       
@@ -2890,11 +2882,11 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       Int_t nMaxima1 = photon1->GetNLM();
       Int_t nMaxima2 = photon2->GetNLM();
       
-      mom = mom1+mom2;
+      fMomentum = fMomentum1+fMomentum2;
       
-      Double_t mass  = mom.M();
-      Double_t epair = mom.E();
-      Float_t ptpair = mom.Pt();
+      Double_t mass  = fMomentum.M();
+      Double_t epair = fMomentum.E();
+      Float_t ptpair = fMomentum.Pt();
       
       if(fFillAllNLMHistograms)
       {
@@ -2934,11 +2926,11 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       //
       // Select good pair (good phi, pt cuts, aperture and invariant mass)
       //
-      if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,GetCalorimeter())) continue;
+      if(!GetNeutralMesonSelection()->SelectPair(fMomentum1, fMomentum2,GetCalorimeter())) continue;
       
       if(GetDebug()>1)
         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",
-               mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
+               fMomentum.Pt(), fMomentum.Phi()*TMath::RadToDeg(), fMomentum.Eta());
       
       //
       // Tag both photons as decay if not done before
@@ -2950,7 +2942,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       {
         if( GetDebug() > 1 )
           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT1 %2.2f; bit requested %d; decay bit1: In %d, ",
-                 mom1.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit1);
+                 fMomentum1.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit1);
         
         GetNeutralMesonSelection()->SetDecayBit(bit1);
         photon1->SetDecayTag(bit1);
@@ -2962,7 +2954,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
 
         //Fill some histograms about shower shape
         if(fFillSelectClHisto && cluster1 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
-          FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
+          FillSelectedClusterHistograms(cluster1, fMomentum1.Pt(), nMaxima1, photon1->GetTag());
         
         if(IsDataMC())
         {
@@ -2982,7 +2974,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       {
         if( GetDebug() > 1 )
           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter - pT2 %2.2f; bit requested %d; decay bit2: In %d, ",
-                 mom2.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit2);
+                 fMomentum2.Pt(), GetNeutralMesonSelection()->GetDecayBit(), bit2);
         
         GetNeutralMesonSelection()->SetDecayBit(bit2);
         photon2->SetDecayTag(bit2);
@@ -2994,7 +2986,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         
         //Fill some histograms about shower shape
         if(fFillSelectClHisto && cluster2 && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
-          FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
+          FillSelectedClusterHistograms(cluster2, fMomentum2.Pt(), nMaxima2, photon2->GetTag());
         
         if(IsDataMC())
         {
@@ -3018,7 +3010,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       FillPileUpHistograms(ptpair,((cluster1->GetTOF()+cluster2->GetTOF())*1e9)/2,cluster1);
       
       //Create AOD for analysis
-      AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
+      AliAODPWG4Particle pi0 = AliAODPWG4Particle(fMomentum);
       
       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
@@ -3055,10 +3047,6 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
   //Check if 2 photons have the mass of the pi0.
   
-  TLorentzVector mom1;
-  TLorentzVector mom2;
-  TLorentzVector mom ;
-  
   // Check calorimeter input
   if(!GetInputAODBranch())
   {
@@ -3100,7 +3088,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++)
   {
     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
-    mom1 = *(photon1->Momentum());
+    fMomentum1 = *(photon1->Momentum());
     
     // Do analysis only when one of the decays is isolated
     // Run AliAnaParticleIsolation before
@@ -3125,13 +3113,13 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
       }
       
-      mom2 = *(photon2->Momentum());
+      fMomentum2 = *(photon2->Momentum());
       
-      mom = mom1+mom2;
+      fMomentum = fMomentum1+fMomentum2;
       
-      Double_t mass  = mom.M();
-      Double_t epair = mom.E();
-      Float_t ptpair = mom.Pt();
+      Double_t mass  = fMomentum.M();
+      Double_t epair = fMomentum.E();
+      Float_t ptpair = fMomentum.Pt();
       
       Int_t nMaxima = photon1->GetNLM();
       if(fFillAllNLMHistograms)
@@ -3153,7 +3141,9 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         Int_t  label2 = photon2->GetLabel();
         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(),"CTS"));
         
-        HasPairSameMCMother(photon1, photon2, label, tag) ;
+        HasPairSameMCMother(photon1->GetLabel(), photon2->GetLabel(),
+                            photon1->GetTag()  , photon2->GetTag(),
+                            label, tag) ;
         mcIndex = GetMCIndex(tag);
       }
       
@@ -3165,10 +3155,10 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
       //
       // Select good pair (good phi, pt cuts, aperture and invariant mass)
       //
-      if(!GetNeutralMesonSelection()->SelectPair(mom1, mom2,GetCalorimeter())) continue ;
+      if(!GetNeutralMesonSelection()->SelectPair(fMomentum1, fMomentum2,GetCalorimeter())) continue ;
       
       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",
-                                mom.Pt(), mom.Phi()*TMath::RadToDeg(), mom.Eta());
+                                fMomentum.Pt(), fMomentum.Phi()*TMath::RadToDeg(), fMomentum.Eta());
       
       //
       // Tag both photons as decay if not done before
@@ -3185,7 +3175,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         
         //Fill some histograms about shower shape
         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
-          FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
+          FillSelectedClusterHistograms(cluster, fMomentum1.Pt(), nMaxima, photon1->GetTag());
         
         if(IsDataMC())
         {
@@ -3214,11 +3204,11 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
       
       // Fill histograms to undertand pile-up before other cuts applied
       // Remember to relax time cuts in the reader
-      if(cluster) FillPileUpHistograms(mom.Pt(),cluster->GetTOF()*1e9,cluster);
+      if(cluster) FillPileUpHistograms(fMomentum.Pt(),cluster->GetTOF()*1e9,cluster);
       
       //Create AOD for analysis
       
-      AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
+      AliAODPWG4Particle pi0 = AliAODPWG4Particle(fMomentum);
       
       if     ( (GetNeutralMesonSelection()->GetParticle()).Contains("Pi0") ) pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
       else if( (GetNeutralMesonSelection()->GetParticle()).Contains("Eta") ) pi0.SetIdentifiedParticleType(AliCaloPID::kEta);
@@ -3276,7 +3266,6 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     return;
   }
 
-  TLorentzVector mom ;
   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
   {
     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
@@ -3292,25 +3281,25 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     Double_t vertex[]={0,0,0};
     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
     {
-      calo->GetMomentum(mom,GetVertex(evtIndex)) ;
+      calo->GetMomentum(fMomentum,GetVertex(evtIndex)) ;
     }//Assume that come from vertex in straight line
     else
     {
-      calo->GetMomentum(mom,vertex) ;
+      calo->GetMomentum(fMomentum,vertex) ;
     }
     
     //If too small or big pt, skip it
-    if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
+    if(fMomentum.E() < GetMinEnergy() || fMomentum.E() > GetMaxEnergy() ) continue ;
     
     //Check acceptance selection
     if(IsFiducialCutOn())
     {
-      Bool_t in = GetFiducialCut()->IsInFiducialCut(mom.Eta(),mom.Phi(),GetCalorimeter()) ;
+      Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),GetCalorimeter()) ;
       if(! in ) continue ;
     }
     
     if(GetDebug() > 1)
-      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
+      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",fMomentum.Pt(),fMomentum.Phi(),fMomentum.Eta());
     
     //Play with the MC stack if available
     //Check origin of the candidates
@@ -3328,7 +3317,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
     if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
-      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      //FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
  
@@ -3337,7 +3326,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     //If too low number of cells, skip it
     if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
     {
-      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      //FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
@@ -3350,7 +3339,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     Double_t tof = calo->GetTOF()*1e9;
     if(tof < fTimeCutMin || tof > fTimeCutMax)
     {
-      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      //FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
 
@@ -3361,12 +3350,14 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     Int_t    absId1   =-1, absId2   =-1;
     Float_t  distbad1 =-1, distbad2 =-1;
     Bool_t   fidcut1  = 0, fidcut2  = 0;
-    TLorentzVector    l1, l2;
 
     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
                                                                                    GetVertex(evtIndex),nMaxima,
-                                                                                   mass,angle,l1,l2,absId1,absId2,
-                                                                                   distbad1,distbad2,fidcut1,fidcut2) ;
+                                                                                   mass,angle,
+                                                                                   fMomentum1,fMomentum2,
+                                                                                   absId1,absId2,
+                                                                                   distbad1,distbad2,
+                                                                                   fidcut1,fidcut2) ;
     
     
     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
@@ -3379,14 +3370,14 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
         Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
                calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
       
-      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      //FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
     //Skip events with too few or too many  NLM
     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
     {
-      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      //FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
@@ -3396,35 +3387,35 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     //Skip matched clusters with tracks
     if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
     {
-      FillRejectedClusterHistograms(mom,tag,nMaxima);
+      FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
     Float_t l0 = calo->GetM02();
-    Float_t e1 = l1.Energy();
-    Float_t e2 = l2.Energy();
-    TLorentzVector l12 = l1+l2;
-    Float_t ptSplit = l12.Pt();
+    Float_t e1 = fMomentum1.Energy();
+    Float_t e2 = fMomentum2.Energy();
+    fMomentum12 = fMomentum1+fMomentum2;
+    Float_t ptSplit = fMomentum12.Pt();
     Float_t  eSplit = e1+e2;
 
     //mass of all clusters
-    fhMass       ->Fill(mom.E() ,mass);
-    fhMassPt     ->Fill(mom.Pt(),mass);
+    fhMass       ->Fill(fMomentum.E() ,mass);
+    fhMassPt     ->Fill(fMomentum.Pt(),mass);
     fhMassSplitPt->Fill(ptSplit ,mass);
-    fhPtLambda0NoSplitCut->Fill(mom.Pt(),l0);
+    fhPtLambda0NoSplitCut->Fill(fMomentum.Pt(),l0);
     
     // Asymmetry of all clusters
     Float_t asy =-10;
     
     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
-    fhAsymmetry->Fill(mom.E(),asy);
+    fhAsymmetry->Fill(fMomentum.E(),asy);
     
     // Divide NLM in 3 cases, 1 local maxima, 2 local maxima, more than 2 local maxima
     Int_t indexMax = -1;
     if     (nMaxima==1) indexMax = 0 ;
     else if(nMaxima==2) indexMax = 1 ;
     else                indexMax = 2 ;
-    fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
+    fhMassPtLocMax[indexMax]->Fill(fMomentum.Pt(),mass);
     
     Int_t   mcIndex   =-1;
     Int_t   noverlaps = 0;
@@ -3436,7 +3427,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       Bool_t ok      = kFALSE;
       Int_t  mcLabel = calo->GetLabel();
       
-      TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
+      fPrimaryMom = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
       
       Int_t mesonLabel = -1;
       
@@ -3444,13 +3435,13 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       {
         if(mcIndex == kmcPi0)
         {
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
-          if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
+          if(fGrandMotherMom.E() > 0 && ok) ptprim =  fGrandMotherMom.Pt();
         }
         else
         {
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
-          if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
+          if(fGrandMotherMom.E() > 0 && ok) ptprim = fGrandMotherMom.Pt();
         }
       }
             
@@ -3458,21 +3449,21 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       Int_t overpdg[nlabels];
       noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
  
-      fhMCMassPt     [mcIndex]->Fill(mom.Pt(),mass);
+      fhMCMassPt     [mcIndex]->Fill(fMomentum.Pt(),mass);
       fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
       if(mcIndex==kmcPi0)
       {
-        fhMCPi0PtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCPi0PtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
         fhMCPi0SplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
-        fhMCPi0PtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCPi0PtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
         fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
 
       }
       else if(mcIndex==kmcEta)
       {
-        fhMCEtaPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCEtaPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
         fhMCEtaSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
-        fhMCEtaPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCEtaPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
         fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
       }
 
@@ -3480,90 +3471,90 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       {
         if(mcIndex==kmcPi0)
         {
-          fhMCPi0PtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCPi0PtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
           fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
         }
         else if(mcIndex==kmcEta)
         {
-          fhMCEtaPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCEtaPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
           fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
         }
         
-        fhMassNoOverlap       ->Fill(mom.E() ,mass);
-        fhMassPtNoOverlap     ->Fill(mom.Pt(),mass);
+        fhMassNoOverlap       ->Fill(fMomentum.E() ,mass);
+        fhMassPtNoOverlap     ->Fill(fMomentum.Pt(),mass);
         fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
         
-        fhMCMassPtNoOverlap     [mcIndex]->Fill(mom.Pt(),mass);
+        fhMCMassPtNoOverlap     [mcIndex]->Fill(fMomentum.Pt(),mass);
         fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
       }
 
-      fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
+      fhMCPtAsymmetry[mcIndex]->Fill(fMomentum.Pt(),asy);
     }
     
     // If cluster does not pass pid, not pi0/eta, skip it.
     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
     {
       if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
-      FillRejectedClusterHistograms(mom,tag,nMaxima);
+      FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
     {
       if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
-      FillRejectedClusterHistograms(mom,tag,nMaxima);
+      FillRejectedClusterHistograms(tag,nMaxima);
       continue ;
     }
     
     if(GetDebug() > 1)
       Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
-             mom.Pt(), idPartType);
+             fMomentum.Pt(), idPartType);
         
     //Mass and asymmetry of selected pairs
-    fhSelectedAsymmetry  ->Fill(mom.E() ,asy );
-    fhSelectedMass       ->Fill(mom.E() ,mass);
-    fhSelectedMassPt     ->Fill(mom.Pt(),mass);
+    fhSelectedAsymmetry  ->Fill(fMomentum.E() ,asy );
+    fhSelectedMass       ->Fill(fMomentum.E() ,mass);
+    fhSelectedMassPt     ->Fill(fMomentum.Pt(),mass);
     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
-    if(fFillAllNLMHistograms) fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
+    if(fFillAllNLMHistograms) fhSelectedMassPtLocMax[indexMax]->Fill(fMomentum.Pt(),mass);
 
     Int_t   nSM  = GetModuleNumber(calo);
     if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0 && fFillAllNLMHistograms)
     {
-      fhSelectedMassPtLocMaxSM   [indexMax][nSM]->Fill(mom.Pt(),mass);
-      fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),l0  );
+      fhSelectedMassPtLocMaxSM   [indexMax][nSM]->Fill(fMomentum.Pt(),mass);
+      fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(fMomentum.Pt(),l0  );
     }
 
     if(IsDataMC())
     {
       if(mcIndex==kmcPi0)
       {
-        fhMCPi0SelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCPi0SelectedPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
         fhMCPi0SelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
-        fhMCPi0SelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCPi0SelectedPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
         fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
       }
       else if(mcIndex==kmcEta)
       {
-        fhMCEtaSelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCEtaSelectedPtRecoPtPrim                     ->Fill(fMomentum.Pt(),ptprim);
         fhMCEtaSelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
-        fhMCEtaSelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCEtaSelectedPtRecoPtPrimLocMax     [indexMax]->Fill(fMomentum.Pt(),ptprim);
         fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
       }
       
       if(noverlaps==0)
       {
-        fhSelectedMassNoOverlap       ->Fill(mom.E() ,mass);
-        fhSelectedMassPtNoOverlap     ->Fill(mom.Pt(),mass);
+        fhSelectedMassNoOverlap       ->Fill(fMomentum.E() ,mass);
+        fhSelectedMassPtNoOverlap     ->Fill(fMomentum.Pt(),mass);
         fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
         
         if(mcIndex==kmcPi0)
         {
-          fhMCPi0SelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCPi0SelectedPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
           fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
         }
         else if(mcIndex==kmcEta)
         {
-          fhMCEtaSelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCEtaSelectedPtRecoPtPrimNoOverlap     ->Fill(fMomentum.Pt(),ptprim);
           fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
         }
       }
@@ -3571,10 +3562,10 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     
     fhSplitE        ->Fill( eSplit);
     fhSplitPt       ->Fill(ptSplit);
-    Float_t phi = mom.Phi();
+    Float_t phi = fMomentum.Phi();
     if(phi<0) phi+=TMath::TwoPi();
     fhSplitPtPhi    ->Fill(ptSplit,phi);
-    fhSplitPtEta    ->Fill(ptSplit,mom.Eta());
+    fhSplitPtEta    ->Fill(ptSplit,fMomentum.Eta());
     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
     
     //Check split-clusters with good time window difference
@@ -3594,45 +3585,45 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       fhMCSplitE        [mcIndex]->Fill( eSplit);
       fhMCSplitPt       [mcIndex]->Fill(ptSplit);
       fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
-      fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,mom.Eta());
+      fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,fMomentum.Eta());
       fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
-      fhMCNLocMaxPt     [mcIndex]->Fill(mom.Pt(),nMaxima);
+      fhMCNLocMaxPt     [mcIndex]->Fill(fMomentum.Pt(),nMaxima);
       
-      fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
+      fhMCSelectedMassPt     [mcIndex]->Fill(fMomentum.Pt(),mass);
       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
-      fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
+      fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(fMomentum.Pt(),mass);
 
       if(noverlaps==0)
       {
-        fhMCSelectedMassPtNoOverlap     [mcIndex]->Fill(mom.Pt(),mass);
+        fhMCSelectedMassPtNoOverlap     [mcIndex]->Fill(fMomentum.Pt(),mass);
         fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
       }
     }
     
     // Remove clusters with NLM=x depeding on a minimim energy cut
-    if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
-    if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
-    if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
+    if(nMaxima == 1 && fNLMECutMin[0] > fMomentum.E()) continue;
+    if(nMaxima == 2 && fNLMECutMin[1] > fMomentum.E()) continue;
+    if(nMaxima >  2 && fNLMECutMin[2] > fMomentum.E()) continue;
 
     //Fill some histograms about shower shape
     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
     {
-      FillSelectedClusterHistograms(calo, mom.Pt(), nMaxima, tag, asy);
+      FillSelectedClusterHistograms(calo, fMomentum.Pt(), nMaxima, tag, asy);
     }
     
     // Fill histograms to undertand pile-up before other cuts applied
     // Remember to relax time cuts in the reader
     Double_t tofcluster   = calo->GetTOF()*1e9;
     
-    FillPileUpHistograms(mom.Pt(),tofcluster,calo);
+    FillPileUpHistograms(fMomentum.Pt(),tofcluster,calo);
     
     if(fFillEMCALBCHistograms && GetCalorimeter()==kEMCAL)
-      FillEMCALBCHistograms(mom.E(), mom.Eta(), mom.Phi(), tofcluster);
+      FillEMCALBCHistograms(fMomentum.E(), fMomentum.Eta(), fMomentum.Phi(), tofcluster);
     
     //-----------------------
     //Create AOD for analysis
     
-    AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
+    AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(fMomentum);
     aodpi0.SetLabel(calo->GetLabel());
     
     //Set the indeces of the original caloclusters
@@ -3728,44 +3719,44 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
         Int_t   momlabel  = -1;
         Bool_t  ok        = kFALSE;
         
-        TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
+        fPrimaryMom = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
         if(!ok) continue;
         
         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
         {
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
-          if(grandmom.E() > 0 && ok)
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
+          if(fGrandMotherMom.E() > 0 && ok)
           {
-            efracMC =  grandmom.E()/ener;
+            efracMC =  fGrandMotherMom.E()/ener;
             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
           }
         }
         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
         {
           fhMCPi0DecayPt->Fill(pt);
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
-          if(grandmom.E() > 0 && ok)
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok,momlabel);
+          if(fGrandMotherMom.E() > 0 && ok)
           {
-            efracMC =  mom.E()/grandmom.E();
+            efracMC =  fPrimaryMom.E()/fGrandMotherMom.E();
             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
           }
         }
         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
         {
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
-          if(grandmom.E() > 0 && ok)
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
+          if(fGrandMotherMom.E() > 0 && ok)
           {
-            efracMC =  grandmom.E()/ener;
+            efracMC =  fGrandMotherMom.E()/ener;
             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
           }
         }
         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
         {
           fhMCEtaDecayPt->Fill(pt);
-          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
-          if(grandmom.E() > 0 && ok)
+          fGrandMotherMom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok,momlabel);
+          if(fGrandMotherMom.E() > 0 && ok)
           {
-            efracMC =  mom.E()/grandmom.E();
+            efracMC =  fPrimaryMom.E()/fGrandMotherMom.E();
             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
           }
         }
index 1790a84..1fad6d2 100755 (executable)
@@ -50,7 +50,7 @@ class AliAnaPi0EbE : public AliAnaCaloTrackCorrBaseClass {
   
   void           FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * c) ;
   
-  void           FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima);
+  void           FillRejectedClusterHistograms(Int_t mctag, Int_t nMaxima);
   
   void           FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt,
                                                Int_t nLocMax,        Int_t tag,
@@ -58,8 +58,8 @@ class AliAnaPi0EbE : public AliAnaCaloTrackCorrBaseClass {
     
   void           FillWeightHistograms(AliVCluster *clus);
     
-  void           HasPairSameMCMother(AliAODPWG4Particle * photon1, 
-                                     AliAODPWG4Particle * photon2, 
+  void           HasPairSameMCMother(Int_t label1 , Int_t label2,
+                                     Int_t tag1   , Int_t tag2,
                                      Int_t & label, Int_t & tag);
   
   void           MakeInvMassInCalorimeter() ;
@@ -159,6 +159,13 @@ class AliAnaPi0EbE : public AliAnaCaloTrackCorrBaseClass {
 
   Bool_t         fCheckSplitDistToBad;     // Check the distance to bad channel and to EMCal borders of split clusters
   
+  TLorentzVector fMomentum;                //! cluster/pi0 momentum
+  TLorentzVector fMomentum1;               //! cluster/photon momentum
+  TLorentzVector fMomentum2;               //! cluster/photon momentum
+  TLorentzVector fMomentum12;              //! cluster/pi0 momentum, sum 1+2
+  TLorentzVector fPrimaryMom;              //! primary momentum
+  TLorentzVector fGrandMotherMom;          //! primary momentum
+
   //Histograms
   
   TH1F         * fhPt  ;                   //! Number of identified  pi0/eta vs pT