]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/GammaConv/AliAnalysisTaskGammaConversion.cxx
Bug fixes for rotation method
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / AliAnalysisTaskGammaConversion.cxx
index 5562a87aa2961399762c15e18fb0f8322b4f397b..d7348d3545ca97bac0d6069ec830b4669d228521 100644 (file)
 //#include "AliCFContainer.h"   // for CF
 #include "AliESDInputHandler.h"
 #include "AliAnalysisManager.h"
+#include "AliAODPWG4Particle.h"
+#include "AliAODPWG4ParticleCorrelation.h"
 #include "AliGammaConversionAODObject.h"
+#include "AliAODConversionParticle.h"
 #include "AliGammaConversionBGHandler.h"
 #include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
 #include "AliKFVertex.h"
 #include "AliGenDPMjetEventHeader.h"
 #include "AliGenEventHeader.h"
 #include <AliMCEventHandler.h>
+#include "TRandom3.h"
+#include "AliTriggerAnalysis.h"
+#include "AliCentrality.h"
+
 class AliESDTrackCuts;
 class AliCFContainer;
 class AliCFManager;
@@ -121,17 +128,27 @@ AliAnalysisTaskSE(),
   fAODPi0(NULL),
   fAODOmega(NULL),
   fAODBranchName("GammaConv"),
+  fOutputAODClassName("AliAODConversionParticle"),
+  fKFCreateAOD(kTRUE),
   fKFForceAOD(kFALSE),
   fKFDeltaAODFileName(""),
   fDoNeutralMesonV0MCCheck(kFALSE),
   fUseTrackMultiplicityForBG(kTRUE),
   fMoveParticleAccordingToVertex(kFALSE),
   fApplyChi2Cut(kFALSE),
-  nRandomEventsForBG(15),
-  nDegreesPMBackground(15),
+  fNRandomEventsForBG(15),
+  fNDegreesPMBackground(15),
   fDoRotation(kTRUE),
   fCheckBGProbability(kTRUE),
-  fKFReconstructedGammasV0Index()
+  fKFReconstructedGammasV0Index(),
+  fRemovePileUp(kFALSE),
+  fSelectV0AND(kFALSE),
+  fTriggerAnalysis(NULL),
+  fMultiplicity(0),
+  fUseMultiplicity(0), 
+  fUseMultiplicityBin(0),
+  fUseCentrality(0), 
+  fUseCentralityBin(0)
 {
   // Default constructor
 
@@ -209,17 +226,27 @@ AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name)
   fAODPi0(NULL),
   fAODOmega(NULL),
   fAODBranchName("GammaConv"),
+  fOutputAODClassName("AliAODConversionParticle"),
+  fKFCreateAOD(kTRUE),
   fKFForceAOD(kFALSE),
   fKFDeltaAODFileName(""),
   fDoNeutralMesonV0MCCheck(kFALSE),
   fUseTrackMultiplicityForBG(kTRUE),
   fMoveParticleAccordingToVertex(kFALSE),
   fApplyChi2Cut(kFALSE),
-  nRandomEventsForBG(15),
-  nDegreesPMBackground(15),
+  fNRandomEventsForBG(15),
+  fNDegreesPMBackground(15),
   fDoRotation(kTRUE),
   fCheckBGProbability(kTRUE),
-  fKFReconstructedGammasV0Index()
+  fKFReconstructedGammasV0Index(),
+  fRemovePileUp(kFALSE),
+  fSelectV0AND(kFALSE),
+  fTriggerAnalysis(NULL),
+  fMultiplicity(0),
+  fUseMultiplicity(0), 
+  fUseMultiplicityBin(0),
+  fUseCentrality(0), 
+  fUseCentralityBin(0)
 {
   // Common I/O in slot 0
   DefineInput (0, TChain::Class());
@@ -278,6 +305,11 @@ AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
   }
   fAODOmega = NULL;
 
+  if(fTriggerAnalysis) {
+    delete fTriggerAnalysis;
+  }
+
+
 }
 
 
@@ -329,7 +361,8 @@ void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
 
 // Using standard function  for setting Cuts
   Bool_t selectPrimaries=kTRUE;
-  fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(selectPrimaries);
+  fEsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
+  fEsdTrackCuts->SetMaxDCAToVertexZ(2);
   fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
   fEsdTrackCuts->SetPtRange(0.15);
   
@@ -452,11 +485,14 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
     fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
   }
  
- if(fKFRecalculatedGammasTClone == NULL){
 if(fKFRecalculatedGammasTClone == NULL){
     fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
   }
 
-       
+  if(fTriggerAnalysis== NULL){
+    fTriggerAnalysis = new AliTriggerAnalysis;
+  }
+
   //clear TClones
   fKFReconstructedGammasTClone->Delete();
   fCurrentEventPosElectronTClone->Delete();
@@ -493,6 +529,19 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
     if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
   }
   */
+  Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0A);
+  Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(fV0Reader->GetESDEvent(), AliTriggerAnalysis::kV0C);
+  Bool_t v0AND = v0A && v0C;
+
+  if(fSelectV0AND && !v0AND){
+    eventQuality=5;
+    fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
+    if(fDoMCTruth){
+      CheckMesonProcessTypeEventQuality(eventQuality);
+    }
+
+    return;
+  }
 
   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
     //    cout<< "Event not taken"<< endl;
@@ -504,6 +553,7 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
     return; // aborts if the primary vertex does not have contributors.
   }
 
+   
 
   if(!fV0Reader->CheckForPrimaryVertexZ() ){
     eventQuality=2;
@@ -513,9 +563,65 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
     }
     return;
   }
+
+  if(fV0Reader->GetESDEvent()->GetPrimaryVertexTracks()->GetNContributors()>0) {
+    fHistograms->FillHistogram("ESD_GlobalPrimaryVtxZ",fV0Reader->GetESDEvent()->GetPrimaryVertex()->GetZ());
+  }else{
+    if(fV0Reader->GetESDEvent()->GetPrimaryVertexSPD()->GetNContributors()>0) {
+      fHistograms->FillHistogram("ESD_SPDPrimaryVtxZ",fV0Reader->GetESDEvent()->GetPrimaryVertex()->GetZ());
+    }
+  }
+
+  if(fRemovePileUp && fV0Reader->GetESDEvent()->IsPileupFromSPD()) {
+    eventQuality=4;
+    fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
+    return;
+  }
+
+  fMultiplicity = fEsdTrackCuts->CountAcceptedTracks(fV0Reader->GetESDEvent());
+
+  if( CalculateMultiplicityBin() != fUseMultiplicityBin){
+    eventQuality=6;
+    fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
+    return;
+  }
+
+  if(fV0Reader->GetIsHeavyIon()){
+    if(fUseCentrality>0){
+      AliCentrality *esdCentrality = fV0Reader->GetESDEvent()->GetCentrality();
+      Int_t centralityC = -1;
+
+      if(fUseCentrality==1){
+       centralityC = esdCentrality->GetCentralityClass10("V0M");
+       if( centralityC != fUseCentralityBin ){
+         eventQuality=7;
+         fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
+         return;
+       }
+      }
+
+      if(fUseCentrality==2){
+       centralityC = esdCentrality->GetCentralityClass10("CL1");
+       if( centralityC != fUseCentralityBin ){
+         eventQuality=7;
+         fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
+         return;
+       }
+      }
+    }
+  }
   eventQuality=3;
   fHistograms->FillHistogram("ESD_EventQuality",eventQuality);
 
+
+
+  fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",fMultiplicity);
+  if (fV0Reader->GetNumberOfContributorsVtx()>=1){
+    fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",fMultiplicity);
+  } 
+
+
+
   // Process the MC information
   if(fDoMCTruth){
     ProcessMCData();
@@ -529,8 +635,11 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
   
 
   //Fill Gamma AOD
-  FillAODWithConversionGammas() ; 
-       
+  if(fKFCreateAOD) {
+    FillAODWithConversionGammas() ; 
+  }
+
+
   // Process reconstructed gammas
   if(fDoNeutralMeson == kTRUE){
     ProcessGammasForNeutralMesonAnalysis();
@@ -585,6 +694,7 @@ void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
 // }
 
 void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evtQ){
+  // Check meson process type event quality
   fStack= MCEvent()->Stack();
   fGCMCEvent=MCEvent();
 
@@ -597,7 +707,17 @@ void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evt
     if(particle->GetPdgCode()!=111){     //Pi0
       continue;
     }
-    if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )   continue;
+
+    Double_t rapidity;
+    if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
+      rapidity=0;
+    }
+    else{
+      rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
+    }  
+
+    if(TMath::Abs(rapidity) > fV0Reader->GetRapidityMesonCut() ) continue; 
+
     if(evtQ==1){
       switch(GetProcessType(fGCMCEvent)){
       case  kProcSD:
@@ -629,6 +749,37 @@ void AliAnalysisTaskGammaConversion::CheckMesonProcessTypeEventQuality(Int_t evt
       }
     }
 
+   if(evtQ==4){
+      switch(GetProcessType(fGCMCEvent)){
+      case  kProcSD:
+       fHistograms->FillHistogram("MC_SD_EvtQ4_Pi0_Pt", particle->Pt());
+       break;
+      case  kProcDD:
+       fHistograms->FillHistogram("MC_DD_EvtQ4_Pi0_Pt", particle->Pt());
+       break;
+      case  kProcND:
+       fHistograms->FillHistogram("MC_ND_EvtQ4_Pi0_Pt", particle->Pt());
+       break;
+      default:
+       AliError("Unknown Process");
+      }
+    }
+
+   if(evtQ==5){
+      switch(GetProcessType(fGCMCEvent)){
+      case  kProcSD:
+       fHistograms->FillHistogram("MC_SD_EvtQ5_Pi0_Pt", particle->Pt());
+       break;
+      case  kProcDD:
+       fHistograms->FillHistogram("MC_DD_EvtQ5_Pi0_Pt", particle->Pt());
+       break;
+      case  kProcND:
+       fHistograms->FillHistogram("MC_ND_EvtQ5_Pi0_Pt", particle->Pt());
+       break;
+      default:
+       AliError("Unknown Process");
+      }
+    }
 
   }
 
@@ -657,8 +808,8 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
     return; // aborts if the primary vertex does not have contributors.
   }
-
-  for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
+  for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
+    //  for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
     TParticle* particle = (TParticle *)fStack->Particle(iTracks);
 
 
@@ -669,70 +820,70 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
     }
                
     ///////////////////////Begin Chic Analysis/////////////////////////////
-               
-    if(particle->GetPdgCode() == 443){//Is JPsi        
-      if(particle->GetNDaughters()==2){
-       if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
-          TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
-
-         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
-         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
-         if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
-           fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
+    if(fDoChic) {
+      if(particle->GetPdgCode() == 443){//Is JPsi      
+       if(particle->GetNDaughters()==2){
+         if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
+            TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
+
+           TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
+           TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
+           if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
+             fHistograms->FillTable("Table_Electrons",3);//e+ e-  from J/Psi inside acceptance
                                        
-         if( TMath::Abs(daug0->Eta()) < 0.9){
-           if(daug0->GetPdgCode() == -11)
-             fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
-           else
-             fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
+           if( TMath::Abs(daug0->Eta()) < 0.9){
+             if(daug0->GetPdgCode() == -11)
+               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
+             else
+               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
                                                
-         }
-         if(TMath::Abs(daug1->Eta()) < 0.9){
-           if(daug1->GetPdgCode() == -11)
-             fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
-           else
-             fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
+           }
+           if(TMath::Abs(daug1->Eta()) < 0.9){
+             if(daug1->GetPdgCode() == -11)
+               fHistograms->FillTable("Table_Electrons",1);//e+  from J/Psi inside acceptance
+             else
+               fHistograms->FillTable("Table_Electrons",2);//e-   from J/Psi inside acceptance
+           }
          }
        }
       }
-    }
-    //              const int CHI_C0   = 10441;
-    //              const int CHI_C1   = 20443;
-    //              const int CHI_C2   = 445
-    if(particle->GetPdgCode() == 22){//gamma from JPsi
-      if(particle->GetMother(0) > -1){
-       if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
-          fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
-          fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
-         if(TMath::Abs(particle->Eta()) < 1.2)
-           fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
+      //              const int CHI_C0   = 10441;
+      //              const int CHI_C1   = 20443;
+      //              const int CHI_C2   = 445
+      if(particle->GetPdgCode() == 22){//gamma from JPsi
+       if(particle->GetMother(0) > -1){
+         if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
+            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
+            fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
+           if(TMath::Abs(particle->Eta()) < 1.2)
+             fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
+         }
        }
       }
-    }
-    if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
-      if( particle->GetNDaughters() == 2){
-       TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
-       TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
+      if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
+       if( particle->GetNDaughters() == 2){
+         TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
+         TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
                                
-       if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
-         if( daug0->GetPdgCode() == 443){
-           TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
-           TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
-           if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
-             fHistograms->FillTable("Table_Electrons",18);
+         if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
+           if( daug0->GetPdgCode() == 443){
+             TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
+             TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
+             if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
+               fHistograms->FillTable("Table_Electrons",18);
                                                
-         }//if
-         else if (daug1->GetPdgCode() == 443){
-           TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
-           TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
-           if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
-             fHistograms->FillTable("Table_Electrons",18);
-         }//else if
-       }//gamma o Jpsi
-      }//GetNDaughters
+           }//if
+           else if (daug1->GetPdgCode() == 443){
+             TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
+             TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
+             if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
+               fHistograms->FillTable("Table_Electrons",18);
+           }//else if
+         }//gamma o Jpsi
+       }//GetNDaughters
+      }
     }
                
-               
     /////////////////////End Chic Analysis////////////////////////////
                
                
@@ -800,6 +951,9 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
         case 331: // Eta'
            fHistograms->FillHistogram("MC_DecayEtapGamma_Pt", particle->Pt());
            break;
+        case 333: // Phi
+           fHistograms->FillHistogram("MC_DecayPhiGamma_Pt", particle->Pt());
+           break;
         }
       }
        
@@ -822,7 +976,10 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
       
        fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);                                        // generated gamma
       }                
-      if(particle->GetMother(0) < 0){   // direct gamma
+      if(particle->GetMother(0) < 0 || //Phojet p+p -> Direct Photons have no mother
+        ((particle->GetMother(0) > -1) && 
+         (TMath::Abs(fStack->Particle(particle->GetMother(0))->GetPdgCode()) < 10)) //Pythia p+p -> Direct Photons have quarks as mother
+        ){   // direct gamma
        fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
        fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
        fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
@@ -910,6 +1067,8 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
       Int_t zBin    = fHistograms->GetZBin(ePos->Vz());
       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());
       Double_t rFMD=30;
+      Double_t rITSTPCMin=50;
+      Double_t rITSTPCMax=80;
 
       TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());       
       
@@ -949,6 +1108,12 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
        fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
       }
 
+      if(ePos->R()>rITSTPCMin  && ePos->R()<rITSTPCMax){
+       TString nameMCMappingITSTPCPhiInZ="";
+       nameMCMappingITSTPCPhiInZ.Form("MC_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
+       fHistograms->FillHistogram(nameMCMappingITSTPCPhiInZ, vtxPos.Phi());
+      }
+
       TString nameMCMappingRInZ="";
       nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
       fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
@@ -1015,6 +1180,14 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
     // process motherparticles (2 gammas as daughters)
     // the motherparticle had already to pass the R and the eta cut, but no line cut.
     // the line cut is just valid for the conversions!
+
+    // OWN primary Pi0 debug ////////////////////////////////////////////////////////////////////////////////////////////
+    if (particle->GetPdgCode()==111){
+      if( TMath::Abs(rapidity) < fV0Reader->GetRapidityMesonCut() ){
+       fHistograms->FillHistogram("MC_Pi0_Pt_vs_Rapid_allDaughters", particle->Pt(),rapidity);
+      }
+    }
+    // end OWN primary Pi0 debug ////////////////////////////////////////////////////////////////////////////////////////
                
     if(particle->GetNDaughters() == 2){
                        
@@ -1161,7 +1334,7 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
              fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
 
              Double_t alfa=0.;
-             if((daughter0->Energy()+daughter1->Energy())!= 0.){
+             if((daughter0->Energy()+daughter1->Energy()) > 0.){
                alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
              }
              fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
@@ -1196,6 +1369,12 @@ void AliAnalysisTaskGammaConversion::ProcessMCData(){
            fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
            fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
 
+           Double_t alfa=0.;
+           if((daughter0->Energy()+daughter1->Energy()) > 0.){
+             alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
+             }
+           fHistograms->FillHistogram("MC_Eta_alpha",alfa);
+
          }
                                        
        }
@@ -1425,9 +1604,18 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
     fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
     fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
-    if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 && fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
-      Double_t EclsToF= (Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
-      fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),EclsToF );
+    Double_t eClsToF= 0;
+    if(!fV0Reader->GetUseCorrectedTPCClsInfo()){
+      if(fV0Reader->GetNegativeTracknTPCFClusters()!=0 ){
+       eClsToF=(Double_t)fV0Reader->GetNegativeTracknTPCClusters()/(Double_t)fV0Reader->GetNegativeTracknTPCFClusters();
+      }
+    }else{
+      eClsToF= fV0Reader->GetNegativeESDTrack()->GetTPCClusterInfo(2,0,fV0Reader->GetFirstTPCRow(fV0Reader->GetXYRadius()));
+    }
+    fHistograms->FillHistogram("ESD_E_nTPCClustersToFP", fV0Reader->GetNegativeTrackP(),eClsToF );
+    fHistograms->FillHistogram("ESD_E_nTPCClustersToFR", fV0Reader->GetXYRadius(),eClsToF );
+
+    if(fV0Reader->GetNegativeTracknTPCClusters()!=0 ){
       fHistograms->FillHistogram("ESD_E_TPCchi2", fV0Reader->GetNegativeTrackTPCchi2()/(Double_t)fV0Reader->GetNegativeTracknTPCClusters());
     }
 
@@ -1439,13 +1627,24 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
     fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
     fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
-    if(fV0Reader->GetPositiveTracknTPCFClusters()!=0 && (Double_t)fV0Reader->GetPositiveTracknTPCClusters()!=0 ){
-      Double_t PclsToF= (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
-      fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), PclsToF);
+    Double_t pClsToF= 0;
+    if(!fV0Reader->GetUseCorrectedTPCClsInfo()){
+      if(fV0Reader->GetPositiveTracknTPCFClusters()!=0){
+       pClsToF = (Double_t)fV0Reader->GetPositiveTracknTPCClusters()/(Double_t)fV0Reader->GetPositiveTracknTPCFClusters();
+      }
+    }else{
+      pClsToF= fV0Reader->GetPositiveESDTrack()->GetTPCClusterInfo(2,0,fV0Reader->GetFirstTPCRow(fV0Reader->GetXYRadius()));
+    }
+
+    fHistograms->FillHistogram("ESD_P_nTPCClustersToFP",fV0Reader->GetPositiveTrackP(), pClsToF);
+    fHistograms->FillHistogram("ESD_P_nTPCClustersToFR",fV0Reader->GetXYRadius(), pClsToF);
+
+    if(fV0Reader->GetPositiveTracknTPCClusters()!=0){
       fHistograms->FillHistogram("ESD_P_TPCchi2", fV0Reader->GetPositiveTrackTPCchi2()/(Double_t)fV0Reader->GetPositiveTracknTPCClusters());
     }
 
 
+
     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
@@ -1496,7 +1695,8 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
     Int_t zBin    = fHistograms->GetZBin(fV0Reader->GetZ());
     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
     Double_t rFMD=30;
-
+    Double_t rITSTPCMin=50;
+    Double_t rITSTPCMax=80;
 
 
     //    Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
@@ -1533,6 +1733,11 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
       fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
     }
 
+    if(fV0Reader->GetXYRadius()>rITSTPCMin && fV0Reader->GetXYRadius()<rITSTPCMax){
+      TString nameESDMappingITSTPCPhiInZ="";
+      nameESDMappingITSTPCPhiInZ.Form("ESD_Conversion_Mapping_ITSTPC_Phi_in_Z_%02d",zBin);
+      fHistograms->FillHistogram(nameESDMappingITSTPCPhiInZ, vtxConv.Phi());
+    }
 
     TString nameESDMappingRInZ="";
     nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
@@ -1574,26 +1779,65 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
                
     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
     if(fDoMCTruth){
+      TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
+      TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
                        
       if(fV0Reader->HasSameMCMother() == kFALSE){
+       fHistograms->FillHistogram("ESD_TrueConvCombinatorial_R", fV0Reader->GetXYRadius());
+       fHistograms->FillHistogram("ESD_TrueConvCombinatorial_Pt", fV0Reader->GetMotherCandidatePt());
+       if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
+         fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_R", fV0Reader->GetXYRadius());
+         fHistograms->FillHistogram("ESD_TrueConvCombinatorialElec_Pt", fV0Reader->GetMotherCandidatePt());
+       }
        continue;
       }
-      TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
-      TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
+      // Moved up to check true electron background
+      //      TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
+      //      TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
 
       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
+       fHistograms->FillHistogram("ESD_TrueConvHadronicBck_R", fV0Reader->GetXYRadius());
+       fHistograms->FillHistogram("ESD_TrueConvHadronicBck_Pt", fV0Reader->GetMotherCandidatePt());
        continue;
       }
 
       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
        continue;
       }
-      if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
-         (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays 
-       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
-         fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
+
+      UInt_t statusPos = fV0Reader->GetPositiveESDTrack()->GetStatus(); 
+      UInt_t statusNeg = fV0Reader->GetNegativeESDTrack()->GetStatus(); 
+      UChar_t itsPixelPos = fV0Reader->GetPositiveESDTrack()->GetITSClusterMap();
+      UChar_t itsPixelNeg = fV0Reader->GetNegativeESDTrack()->GetITSClusterMap();
+
+      // Using the UniqueID Phojet does not get the Dalitz right
+      //      if( (negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4) ||
+      //  (negativeMC->GetUniqueID() == 0 && positiveMC->GetUniqueID() ==0) ){// fill r distribution for Dalitz decays 
+      if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
+       fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
+       fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_R", fV0Reader->GetXYRadius());
+       //--------Histos for HFE 
+
+       if(statusPos & AliESDtrack::kTOFpid){
+         fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_R", fV0Reader->GetXYRadius());
+         if( TESTBIT(itsPixelPos, 0) ){ 
+           fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SinglePos_kFirst_R", fV0Reader->GetXYRadius());
+         }
+       }
+       if(statusNeg & AliESDtrack::kTOFpid){
+         fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_R", fV0Reader->GetXYRadius());
+         if( TESTBIT(itsPixelNeg, 0) ){ 
+           fHistograms->FillHistogram("ESD_TrueConvDalitzPi0_SingleNeg_kFirst_R", fV0Reader->GetXYRadius());
+         }
        }
+       //--------------------------------------------------------
+
       }
+      if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 221){ //eta
+       fHistograms->FillHistogram("ESD_TrueConvDalitzEta_R", fV0Reader->GetXYRadius());
+      }
+
+      //}
 
       if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
        continue;
@@ -1627,12 +1871,30 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
        fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
        fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
                                
-                               
+       fHistograms->FillHistogram("ESD_TrueConversion_E_nTPCClustersToFR", fV0Reader->GetXYRadius(),eClsToF );
+       fHistograms->FillHistogram("ESD_TrueConversion_P_nTPCClustersToFR",fV0Reader->GetXYRadius(), pClsToF);
+
        fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
        fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
        fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
        fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
 
+       //----Histos for HFE-------------------------------------- 
+
+       if(statusPos & AliESDtrack::kTOFpid){
+         fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
+         if( TESTBIT(itsPixelPos, 0) ){ 
+           fHistograms->FillHistogram("ESD_TrueConversion_SinglePos_kFirst_R", positiveMC->R(),fV0Reader->GetPositiveMCParticle()->Pt());
+         }
+       }
+       if(statusNeg & AliESDtrack::kTOFpid){
+         fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
+         if( TESTBIT(itsPixelNeg, 0) ){ 
+            fHistograms->FillHistogram("ESD_TrueConversion_SingleNeg_kFirst_R", negativeMC->R(),fV0Reader->GetNegativeMCParticle()->Pt());
+          }
+       }
+       //--------------------------------------------------------
+
        fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
        fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
        fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
@@ -1691,9 +1953,9 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
  //    AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
        UInt_t kTRDoutN =  (statusN & AliESDtrack::kTRDout);
 
-       Int_t ITSclsE= fV0Reader->GetNegativeTracknITSClusters();
+       Int_t nITSclsE= fV0Reader->GetNegativeTracknITSClusters();
        // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
-        switch(ITSclsE){
+        switch(nITSclsE){
          case 0: // 0 ITS clusters
                fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
            break;
@@ -1737,9 +1999,9 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
 //     AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
        UInt_t kTRDoutP =  (statusP & AliESDtrack::kTRDout);
        
-       Int_t ITSclsP = fV0Reader->GetPositiveTracknITSClusters();
+       Int_t nITSclsP = fV0Reader->GetPositiveTracknITSClusters();
        // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
-        switch(ITSclsP){
+        switch(nITSclsP){
          case 0: // 0 ITS clusters
                fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
            break;
@@ -1802,43 +2064,79 @@ void AliAnalysisTaskGammaConversion::ProcessV0s(){
   fV0Reader->ResetV0IndexNumber();
 }
 
+void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODPWG4Particle & particle) {
+  //See header file for documentation
+
+  Int_t i = branch->GetEntriesFast();
+  if(! (fOutputAODClassName.Contains("Correlation")) ) {
+    new((*branch)[i])  AliAODPWG4Particle(particle);
+  } else {
+    new((*branch)[i])  AliAODPWG4ParticleCorrelation(particle);
+  }
+}
+
+void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliGammaConversionAODObject & particle) {
+  //See header file for documentation
+
+  Int_t i = branch->GetEntriesFast();
+  new((*branch)[i])  AliGammaConversionAODObject(particle);
+}
+
+void AliAnalysisTaskGammaConversion::AddToAODBranch(TClonesArray * branch, AliAODConversionParticle & particle) {
+  //See header file for documentation
+
+  Int_t i = branch->GetEntriesFast();
+  new((*branch)[i])  AliAODConversionParticle(particle);
+}
+
+
 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
   // Fill AOD with reconstructed Gamma
        
   for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
-    //  for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
-    //Create AOD particle object from AliKFParticle
-    //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
-    //but this means that I have to work a little bit more in my side.
-    //AODPWG4Particle objects are simpler and lighter, I think
-    /*
     AliKFParticle * gammakf = dynamic_cast<AliKFParticle*>(fKFReconstructedGammasTClone->At(gammaIndex));
-    AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
-    //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
-    gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
-    gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
-    gamma.SetPdg(AliPID::kEleCon); //photon id
-    gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
-    gamma.SetChi2(gammakf->Chi2());
-    Int_t i = fAODBranch->GetEntriesFast();
-    new((*fAODBranch)[i])  AliAODPWG4Particle(gamma);
-    */
-
-    AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
-    AliGammaConversionAODObject aodObject;
-    aodObject.SetPx(gammakf->GetPx());
-    aodObject.SetPy(gammakf->GetPy());
-    aodObject.SetPz(gammakf->GetPz());
-    aodObject.SetLabel1(fElectronv1[gammaIndex]);
-    aodObject.SetLabel2(fElectronv2[gammaIndex]);
-    aodObject.SetChi2(gammakf->Chi2());
-    aodObject.SetE(gammakf->E());
-    Int_t i = fAODGamma->GetEntriesFast();
-    new((*fAODGamma)[i])  AliGammaConversionAODObject(aodObject);
-  }
+    
+    if(fOutputAODClassName.Contains("AliAODPWG4Particle")) {
+      AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
+      //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
+      gamma.SetTrackLabel( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
+      gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
+      gamma.SetPdg(AliPID::kEleCon); //photon id
+      gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
+      //PH    gamma.SetChi2(gammakf->Chi2());
+      
+      AddToAODBranch(fAODGamma, gamma);
+      
+    } else if(fOutputAODClassName.Contains("ConversionParticle")) {
+      TLorentzVector momentum(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
+      AliAODConversionParticle gamma = AliAODConversionParticle(momentum);
+      //gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
+      gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] ); //How to get the MC label of the 2 electrons that form the gamma?
+      //gamma.SetPdg(AliPID::kEleCon); //photon id
+      //gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
+      gamma.SetChi2(gammakf->Chi2());
+      gamma.SetTrackLabels( fElectronv1[gammaIndex], fElectronv2[gammaIndex] );
+      gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));    
+      AddToAODBranch(fAODGamma, gamma);
+      
+      
+
+    } else {
+      AliGammaConversionAODObject gamma;
+      gamma.SetPx(gammakf->GetPx());
+      gamma.SetPy(gammakf->GetPy());
+      gamma.SetPz(gammakf->GetPz());
+      gamma.SetE(gammakf->GetE());
+      gamma.SetLabel1(fElectronv1[gammaIndex]);
+      gamma.SetLabel2(fElectronv2[gammaIndex]);
+      gamma.SetChi2(gammakf->Chi2());
+      gamma.SetE(gammakf->E());
+      gamma.SetESDEvent(dynamic_cast<AliESDEvent*>(InputEvent()));
+      AddToAODBranch(fAODGamma, gamma);
+    }
        
+  }
 }
-
 void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
   // omega meson analysis pi0+gamma decay
   for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
@@ -1857,7 +2155,7 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
       omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
 
       if ( massOmegaCandidate > 733 && massOmegaCandidate < 833 ) {
-       AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
+       //AddOmegaToAOD(&omegaCandidate, massOmegaCandidate, firstPi0Index, firstGammaIndex);
       }
       
       fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
@@ -1892,7 +2190,7 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
          omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
 
          if ( massOmegaCandidatePipPinPi0 > 733 && massOmegaCandidatePipPinPi0 < 833 ) {
-           AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
+           // AddOmegaToAOD(&omegaCandidatePipPinPi0, massOmegaCandidatePipPinPi0, -1, -1);
          }
          
          fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
@@ -1901,9 +2199,10 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
          //  delete omegaCandidatePipPinPi0;
        }
       }
-    } // checking ig gammajet because in that case the chargedparticle list is created
 
+      if (posPiKF) delete posPiKF; posPiKF=NULL;     if (negPiKF) delete negPiKF; negPiKF=NULL;
 
+    } // checking ig gammajet because in that case the chargedparticle list is created
 
   }
 
@@ -1956,19 +2255,19 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
 }
 
 
-void AliAnalysisTaskGammaConversion::AddOmegaToAOD(AliKFParticle * omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
+void AliAnalysisTaskGammaConversion::AddOmegaToAOD(const AliKFParticle * const omegakf, Double_t mass, Int_t omegaDaughter, Int_t gammaDaughter) {
   //See header file for documentation
   AliGammaConversionAODObject omega;
-  omega.SetPx(omegakf->Px());
-  omega.SetPy(omegakf->Py());
-  omega.SetPz(omegakf->Pz());
+  omega.SetPx(omegakf->GetPx());
+  omega.SetPy(omegakf->GetPy());
+  omega.SetPz(omegakf->GetPz());
   omega.SetChi2(omegakf->GetChi2());
-  omega.SetE(omegakf->E());
+  omega.SetE(omegakf->GetE());
   omega.SetIMass(mass);
   omega.SetLabel1(omegaDaughter);
-  //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
+  // //dynamic_cast<AliAODPWG4Particle*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
   omega.SetLabel2(gammaDaughter);
-  new((*fAODOmega)[fAODOmega->GetEntriesFast()])  AliGammaConversionAODObject(omega);
+  AddToAODBranch(fAODOmega, omega);
 }
 
 
@@ -2053,6 +2352,13 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
            fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
            fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
            fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
+           if(massTwoGammaCandidate>0.1 && massTwoGammaCandidate<0.15){
+             fHistograms->FillHistogram("ESD_Mother_alfa_Pi0", alfa);
+           }
+           if(massTwoGammaCandidate>0.5 && massTwoGammaCandidate<0.57){
+             fHistograms->FillHistogram("ESD_Mother_alfa_Eta", alfa);
+           }
+
            fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!
            fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
            fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
@@ -2114,6 +2420,12 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
                      gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
                    }
                  }
+                 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==111){
+                     gamma1MotherLabel=-111;
+                  }
+                 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==221){
+                     gamma1MotherLabel=-221;
+                  }
                }
              }
              Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
@@ -2133,7 +2445,14 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
                        gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
                      }
                    }
-                 }
+                   if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==111){
+                      gamma2MotherLabel=-111;
+                    }
+                   if(fV0Reader->GetMotherMCParticle()->GetPdgCode() ==221){
+                      gamma2MotherLabel=-221;
+                    }
+                   
+                 }
                }
                if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
                  if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
@@ -2159,9 +2478,12 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
 
                    if(!isRealPi0 && !isRealEta){
                      if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
-                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }else{
-                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
+                     }
+                     if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221){
+                       fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }
                    }
 
@@ -2180,9 +2502,12 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
                    }
                    if(!isRealPi0 && !isRealEta){
                      if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
-                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }else{
-                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
+                     }
+                     if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221){
+                       fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }
                    }
                  }
@@ -2197,14 +2522,17 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
                      fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
                      fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
                      if(gamma1MotherLabel > fV0Reader->GetMCStack()->GetNprimary()){
-                       fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TruePi0Sec_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }
                    }
                    if(!isRealPi0 && !isRealEta){
                      if(gamma1MotherLabel>-1 && gamma2MotherLabel>-1){
-                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckGG_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }else{
-                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass",massTwoGammaCandidate);
+                       fHistograms->FillHistogram("ESD_TrueBckCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
+                     }
+                     if(gamma1MotherLabel==-111 || gamma2MotherLabel==-111 || gamma1MotherLabel==-221 || gamma2MotherLabel==-221 ){
+                       fHistograms->FillHistogram("ESD_TruePi0DalitzCont_InvMass_vs_Pt",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
                      }
                    }
                  }
@@ -2233,16 +2561,17 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
 
            Double_t lowMassPi0=0.1;
            Double_t highMassPi0=0.15;
-           if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
+           if ( ( massTwoGammaCandidate > lowMassPi0) && (massTwoGammaCandidate < highMassPi0) ){
              new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()])  AliKFParticle(*twoGammaCandidate);
              fGammav1.push_back(firstGammaIndex);
              fGammav2.push_back(secondGammaIndex);
-             AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
+             if( fKFCreateAOD ) {
+               AddPionToAOD(twoGammaCandidate, massTwoGammaCandidate, firstGammaIndex, secondGammaIndex);
+             }
            }
          }
 
        }
-         //}
        delete twoGammaCandidate;
     }
   }
@@ -2250,24 +2579,31 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
 
 void AliAnalysisTaskGammaConversion::AddPionToAOD(AliKFParticle * pionkf, Double_t mass, Int_t daughter1, Int_t daughter2) {
   //See header file for documentation
-  AliGammaConversionAODObject pion;
-  pion.SetPx(pionkf->Px());
-  pion.SetPy(pionkf->Py());
-  pion.SetPz(pionkf->Pz());
-  pion.SetChi2(pionkf->GetChi2());
-  pion.SetE(pionkf->E());
-  pion.SetIMass(mass);
-  pion.SetLabel1(daughter1);
-  //dynamic_cast<AliGammaConversionAODObject*>(fAODBranch->At(daughter1))->SetTagged(kTRUE);
-  pion.SetLabel2(daughter2);
-  new((*fAODPi0)[fAODPi0->GetEntriesFast()])  AliGammaConversionAODObject(pion);
-
+  if(fOutputAODClassName.Contains("AODObject")) {
+    AliGammaConversionAODObject pion;
+    pion.SetPx(pionkf->GetPx());
+    pion.SetPy(pionkf->GetPy());
+    pion.SetPz(pionkf->GetPz());
+    pion.SetChi2(pionkf->GetChi2());
+    pion.SetE(pionkf->GetE());
+    pion.SetIMass(mass);
+    pion.SetLabel1(daughter1);
+    pion.SetLabel2(daughter2);
+    AddToAODBranch(fAODPi0, pion);
+  } else {
+    TLorentzVector momentum(pionkf->Px(),pionkf->Py(),pionkf->Pz(), pionkf->E());
+    AliAODConversionParticle pion = AliAODConversionParticle(momentum);
+    pion.SetTrackLabels( daughter1, daughter2 );
+    pion.SetChi2(pionkf->GetChi2());
+    AddToAODBranch(fAODPi0, pion);
+    
+  }
 }
 
 
-
-void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
 /*
+  void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
+
   // see header file for documentation
   // Analyse Pi0 with one photon from Phos and 1 photon from conversions
        
@@ -2290,29 +2626,29 @@ void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysi
 
 
   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
-    clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
-    if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
-    clu ->GetMomentum(pPHOS ,vtx);
-    for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
-      AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
-      gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
-      gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
-      pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
-      fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
-      fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
-
-      TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
-      TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
-      Double_t opanConvPHOS= v3D0.Angle(v3D1);
-      if ( opanConvPHOS < 0.35){
-       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
-      }else{
-       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
-      }
+  clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
+  if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
+  clu ->GetMomentum(pPHOS ,vtx);
+  for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
+  AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
+  gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
+  gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
+  pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
+  fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
+
+  TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
+  TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
+  Double_t opanConvPHOS= v3D0.Angle(v3D1);
+  if ( opanConvPHOS < 0.35){
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
+  }else{
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
+  }
 
-    }
+  }
 
-    //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
+  //  Now the LorentVector pPHOS is obtained and can be paired with the converted proton
   }
   //==== End of the PHOS cluster selection ============
   TLorentzVector pEMCAL;
@@ -2320,122 +2656,142 @@ void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysi
   TLorentzVector pi0GammaConvEMCAL;
   TLorentzVector pi0GammaConvEMCALBck;
 
-   for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
-    clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
-    if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
-    if (clu->GetNCells() <= 1) continue;
-    if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
-
-    clu ->GetMomentum(pEMCAL ,vtx);
-    for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
-      AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
-      gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
-                            twoGammaDecayCandidateDaughter0->Py(),
-                            twoGammaDecayCandidateDaughter0->Pz(),0.);
-      gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
-      pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
-      fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
-      fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
-      TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
-                   twoGammaDecayCandidateDaughter0->Py(),
-                   twoGammaDecayCandidateDaughter0->Pz());
-      TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
-
+  for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
+  clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
+  if ( !clu->IsEMCAL()  || clu->E()<0.1 ) continue;
+  if (clu->GetNCells() <= 1) continue;
+  if ( clu->GetTOF()*1e9 < 550  || clu->GetTOF()*1e9 > 750) continue;
 
-      Double_t opanConvEMCAL= v3D0.Angle(v3D1);
-      if ( opanConvEMCAL < 0.35){
-       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
-      }else{
-       fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
-      }
+  clu ->GetMomentum(pEMCAL ,vtx);
+  for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
+  AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
+  gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
+  twoGammaDecayCandidateDaughter0->Py(),
+  twoGammaDecayCandidateDaughter0->Pz(),0.);
+  gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
+  pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
+  fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
+  TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
+  twoGammaDecayCandidateDaughter0->Py(),
+  twoGammaDecayCandidateDaughter0->Pz());
+  TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
+
+
+  Double_t opanConvEMCAL= v3D0.Angle(v3D1);
+  if ( opanConvEMCAL < 0.35){
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
+  }else{
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
+  }
 
-    }
-    if(fCalculateBackground){
-      for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
-       AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
-       for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
-         AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
-         gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
-                                   previousGoodV0.Py(),
-                                   previousGoodV0.Pz(),0.);
-         pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
-         fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
-         fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
-                                    pi0GammaConvEMCALBck.Pt());
-       }
-      }
+  }
+  if(fCalculateBackground){
+  for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
+  AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
+  for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
+  AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
+  gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
+  previousGoodV0.Py(),
+  previousGoodV0.Pz(),0.);
+  pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
+  fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
+  fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
+  pi0GammaConvEMCALBck.Pt());
+  }
+  }
       
-      //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
-    } // end of checking if background photons are available
-   }
+  //  Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
+  } // end of checking if background photons are available
+  }
   //==== End of the PHOS cluster selection ============
+
+  }
 */
-}
 
-void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle *particle,AliGammaConversionBGHandler::GammaConversionVertex *vertex){
+void AliAnalysisTaskGammaConversion::MoveParticleAccordingToVertex(AliKFParticle * particle,const AliGammaConversionBGHandler::GammaConversionVertex *vertex){
   //see header file for documentation
 
-  Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
-  Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
-  Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
+   Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
+   Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
+   Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
   
   //  cout<<"dx, dy, dz: ["<<dx<<","<<dy<<","<<dz<<"]"<<endl;
-  particle->X() = particle->GetX() - dx;
-  particle->Y() = particle->GetY() - dy;
-  particle->Z() = particle->GetZ() - dz;
+   particle->X() = particle->GetX() - dx;
+   particle->Y() = particle->GetY() - dy;
+   particle->Z() = particle->GetZ() - dz;
 }
 
 void AliAnalysisTaskGammaConversion::RotateKFParticle(AliKFParticle * kfParticle,Double_t angle){
+  // Before rotate needs to be moved to position 0,0,0, ; move back after rotation
+  Double_t dx = fESDEvent->GetPrimaryVertex()->GetX()-0.;
+  Double_t dy = fESDEvent->GetPrimaryVertex()->GetY()-0.;
+  Double_t dz = fESDEvent->GetPrimaryVertex()->GetZ()-0.;
   
+  kfParticle->X() = kfParticle->GetX() - dx;
+  kfParticle->Y() = kfParticle->GetY() - dy;
+  kfParticle->Z() = kfParticle->GetZ() - dz;
+
+
+  // Rotate the kf particle
   Double_t c = cos(angle);
   Double_t s = sin(angle);
   
-  Double_t A[7][ 7];
-  for( Int_t i=0; i<7; i++ ){
-    for( Int_t j=0; j<7; j++){
-      A[i][j] = 0;
+  Double_t mA[8][ 8];
+  for( Int_t i=0; i<8; i++ ){
+    for( Int_t j=0; j<8; j++){
+      mA[i][j] = 0;
     }
   }
-  for( int i=0; i<7; i++ ){
-    A[i][i] = 1;
+  for( int i=0; i<8; i++ ){
+    mA[i][i] = 1;
   }
-  A[0][0] =  c;   A[0][1] = s;
-  A[1][0] = -s;  A[1][1] = c;
-  A[3][3] =  c;  A[3][4] = s;
-  A[4][3] = -s;  A[4][4] = c;
+  mA[0][0] =  c;  mA[0][1] = s;
+  mA[1][0] = -s;  mA[1][1] = c;
+  mA[3][3] =  c;  mA[3][4] = s;
+  mA[4][3] = -s;  mA[4][4] = c;
   
-  Double_t AC[7][7];
-  Double_t Ap[7];
+  Double_t mAC[8][8];
+  Double_t mAp[8];
   
-  for( Int_t i=0; i<7; i++ ){
-    Ap[i] = 0;
-    for( Int_t k=0; k<7; k++){
-      Ap[i]+=A[i][k] * kfParticle->GetParameter(k);
+  for( Int_t i=0; i<8; i++ ){
+    mAp[i] = 0;
+    for( Int_t k=0; k<8; k++){
+      mAp[i]+=mA[i][k] * kfParticle->GetParameter(k);
     }
   }
   
-  for( Int_t i=0; i<7; i++){
-    kfParticle->Parameter(i) = Ap[i];
+  for( Int_t i=0; i<8; i++){
+    kfParticle->Parameter(i) = mAp[i];
   }
 
-  for( Int_t i=0; i<7; i++ ){
-    for( Int_t j=0; j<7; j++ ){
-      AC[i][j] = 0;
-      for( Int_t k=0; k<7; k++ ){
-       AC[i][j]+= A[i][k] * kfParticle->GetCovariance(k,j);
+  for( Int_t i=0; i<8; i++ ){
+    for( Int_t j=0; j<8; j++ ){
+      mAC[i][j] = 0;
+      for( Int_t k=0; k<8; k++ ){
+       mAC[i][j]+= mA[i][k] * kfParticle->GetCovariance(k,j);
       }
     }
   }
 
-  for( Int_t i=0; i<7; i++ ){
+  for( Int_t i=0; i<8; i++ ){
     for( Int_t j=0; j<=i; j++ ){
       Double_t xx = 0;
-      for( Int_t k=0; k<7; k++){
-       xx+= AC[i][k]*A[j][k];
+      for( Int_t k=0; k<8; k++){
+       xx+= mAC[i][k]*mA[j][k];
       }
       kfParticle->Covariance(i,j) = xx;
     }
   }
+
+  Double_t dx1 = 0.-fESDEvent->GetPrimaryVertex()->GetX();
+  Double_t dy1 = 0.-fESDEvent->GetPrimaryVertex()->GetY();
+  Double_t dz1 = 0.-fESDEvent->GetPrimaryVertex()->GetZ();
+  
+  kfParticle->X() = kfParticle->GetX() - dx1;
+  kfParticle->Y() = kfParticle->GetY() - dy1;
+  kfParticle->Z() = kfParticle->GetZ() - dz1;
+
 }
 
 
@@ -2457,11 +2813,12 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
   }
 
   if(fDoRotation == kTRUE){
-    TRandom3 *random = new TRandom3();
+    TRandom3 *random = new TRandom3(0);
+
     for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
       AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent)); 
       for(Int_t iCurrent2=iCurrent+1;iCurrent2<currentEventV0s->GetEntriesFast();iCurrent2++){
-       for(Int_t nRandom=0;nRandom<nRandomEventsForBG;nRandom++){
+       for(Int_t nRandom=0;nRandom<fNRandomEventsForBG;nRandom++){
        
          AliKFParticle currentEventGoodV02 = *(AliKFParticle *)(currentEventV0s->At(iCurrent2));
 
@@ -2479,10 +2836,10 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
            delete backgroundCandidateProb;
          }
        
-         Double_t nRadiansPM = nDegreesPMBackground*TMath::Pi()/180;
+         Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
          
          Double_t rotationValue = random->Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
-
+         
          RotateKFParticle(&currentEventGoodV02,rotationValue);
 
          AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,currentEventGoodV02);
@@ -2491,7 +2848,6 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
          Double_t widthBG = 0.;
          Double_t chi2BG =10000.;      
          backgroundCandidate->GetMass(massBG,widthBG);
-
          //      if(backgroundCandidate->GetNDF()>0){
          chi2BG = backgroundCandidate->GetChi2();
          if((chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson())  || fApplyChi2Cut == kFALSE){
@@ -2539,6 +2895,12 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
              fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
              fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
 
+             if(massBG>0.1 && massBG<0.15){
+               fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
+             }
+             if(massBG>0.5 && massBG<0.57){
+               fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
+             }
 
              if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(currentEventGoodV02.GetEta())<0.9 ){
                fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
@@ -2609,6 +2971,7 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
            Double_t widthBG = 0.;
            Double_t chi2BG =10000.;    
            backgroundCandidate->GetMass(massBG,widthBG);
+
            //  if(backgroundCandidate->GetNDF()>0){
            //    chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
            chi2BG = backgroundCandidate->GetChi2();
@@ -2661,6 +3024,12 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
                fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
                fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
 
+               if(massBG>0.1 && massBG<0.15){
+                 fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
+               }
+               if(massBG>0.5 && massBG<0.57){
+                 fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
+               }
 
                if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
                  fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
@@ -2726,6 +3095,7 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
              Double_t widthBG = 0.;
              Double_t chi2BG =10000.;  
              backgroundCandidate->GetMass(massBG,widthBG);
+
              /*            if(backgroundCandidate->GetNDF()>0){
                            chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
                            {//remember to remove
@@ -2782,7 +3152,13 @@ void AliAnalysisTaskGammaConversion::CalculateBackground(){
 
                  fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
 
-                 
+                 if(massBG>0.1 && massBG<0.15){
+                   fHistograms->FillHistogram("ESD_Background_alfa_Pi0", alfa);
+                 }
+                 if(massBG>0.5 && massBG<0.57){
+                   fHistograms->FillHistogram("ESD_Background_alfa_Eta", alfa);
+                 }
+
                  if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 &&  TMath::Abs(previousGoodV0.GetEta())<0.9 ){
                    fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
                    fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
@@ -2854,7 +3230,7 @@ void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
 }
 
 //____________________________________________________________________
-Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
+Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(const AliESDtrack *const track)
 {
 //
 // check whether particle has good DCAr(Pt) impact
@@ -2895,14 +3271,16 @@ void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
       numberOfESDTracks++;
     }
   }
-  fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
-
-  if (fV0Reader->GetNumberOfContributorsVtx()>=1){
-    fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
-  } 
+// Moved to UserExec using CountAcceptedTracks function. runjet is not needed by default
+//   fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
+//   cout<<"esdtracks::"<< numberOfESDTracks<<endl;
+//   if (fV0Reader->GetNumberOfContributorsVtx()>=1){
+//     fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
+//   } 
 }
 void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
-  
+  //recalculates v0 for gamma
+
  Double_t massE=0.00051099892;
  TLorentzVector curElecPos;
  TLorentzVector curElecNeg;
@@ -3284,29 +3662,33 @@ void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
 
 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
 {
-  //AOD
-  if(!fAODGamma) fAODGamma = new TClonesArray("AliGammaConversionAODObject", 0);
-  else fAODGamma->Delete();
-  fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
   
-  if(!fAODPi0) fAODPi0 = new TClonesArray("AliGammaConversionAODObject", 0);
-  else fAODPi0->Delete();
-  fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
-
-  if(!fAODOmega) fAODOmega = new TClonesArray("AliGammaConversionAODObject", 0);
-  else fAODOmega->Delete();
-  fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
-
-  //If delta AOD file name set, add in separate file. Else add in standard aod file. 
-  if(fKFDeltaAODFileName.Length() > 0) {
-    AddAODBranch("TClonesArray", &fAODGamma, fKFDeltaAODFileName.Data());
-    AddAODBranch("TClonesArray", &fAODPi0, fKFDeltaAODFileName.Data());
-    AddAODBranch("TClonesArray", &fAODOmega, fKFDeltaAODFileName.Data());
-    AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fKFDeltaAODFileName.Data());
-  } else  {
-    AddAODBranch("TClonesArray", &fAODGamma);
-    AddAODBranch("TClonesArray", &fAODPi0);
-    AddAODBranch("TClonesArray", &fAODOmega);
+  if(fKFCreateAOD) {
+
+    //AOD
+    if(!fAODGamma) fAODGamma = new TClonesArray(fOutputAODClassName.Data(), 0);
+    else fAODGamma->Delete();
+    fAODGamma->SetName(Form("%s_gamma", fAODBranchName.Data()));
+    
+    if(!fAODPi0) fAODPi0 = new TClonesArray(fOutputAODClassName.Data(), 0);
+    else fAODPi0->Delete();
+    fAODPi0->SetName(Form("%s_Pi0", fAODBranchName.Data()));
+    
+    if(!fAODOmega) fAODOmega = new TClonesArray(fOutputAODClassName.Data(), 0);
+    else fAODOmega->Delete();
+    fAODOmega->SetName(Form("%s_Omega", fAODBranchName.Data()));
+    
+    //If delta AOD file name set, add in separate file. Else add in standard aod file. 
+    if(GetDeltaAODFileName().Length() > 0) {
+      AddAODBranch("TClonesArray", &fAODGamma, GetDeltaAODFileName().Data());
+      AddAODBranch("TClonesArray", &fAODPi0, GetDeltaAODFileName().Data());
+      AddAODBranch("TClonesArray", &fAODOmega, GetDeltaAODFileName().Data());
+      AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(GetDeltaAODFileName().Data());
+    } else  {
+      AddAODBranch("TClonesArray", &fAODGamma);
+      AddAODBranch("TClonesArray", &fAODPi0);
+      AddAODBranch("TClonesArray", &fAODOmega);
+    }
   }
 
   // Create the output container
@@ -3340,7 +3722,7 @@ void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
   fOutputContainer->SetName(GetName());
 }
 
-Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
+Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(const TParticle* const daughter0, const TParticle* const daughter1) const{
   //helper function
   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
@@ -3895,7 +4277,7 @@ void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negat
 }
 
 
-void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
+void  AliAnalysisTaskGammaConversion::GetPID(const AliESDtrack *track, Stat_t &pid, Stat_t &weight)
 {
   // see header file for documentation
   pid = -1;
@@ -3943,7 +4325,7 @@ void  AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, St
   pid = ipid;
   weight = max;
 }
-double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
+double AliAnalysisTaskGammaConversion::GetSigmaToVertex(const AliESDtrack* t)
 {
   // Calculates the number of sigma to the vertex.
        
@@ -4004,8 +4386,7 @@ TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *con
        
   //  return tlVtrack;
 }
-
-Int_t AliAnalysisTaskGammaConversion::GetProcessType(AliMCEvent * mcEvt) {
+Int_t AliAnalysisTaskGammaConversion::GetProcessType(const AliMCEvent * mcEvt) {
 
   // Determine if the event was generated with pythia or phojet and return the process type
 
@@ -4060,3 +4441,21 @@ Int_t AliAnalysisTaskGammaConversion::GetProcessType(AliMCEvent * mcEvt) {
   AliError(Form("Unknown header: %s", htmp->IsA()->GetName()));
   return kProcUnknown;
 }
+
+
+Int_t AliAnalysisTaskGammaConversion::CalculateMultiplicityBin(){
+  // Get Centrality bin
+
+  Int_t multiplicity = 0;
+
+  if ( fUseMultiplicity == 1 ) {
+
+    if (fMultiplicity>= 0 && fMultiplicity<= 5) multiplicity=1;
+    if (fMultiplicity>= 6 && fMultiplicity<= 9) multiplicity=2;
+    if (fMultiplicity>=10 && fMultiplicity<=14) multiplicity=3;
+    if (fMultiplicity>=15 && fMultiplicity<=22) multiplicity=4;
+    if (fMultiplicity>=23 )                     multiplicity=5;
+
+  }
+  return multiplicity;
+}