]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Reader: Add option to remove or not event with primary vertex not reconstructed
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 16 Apr 2011 18:24:38 +0000 (18:24 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 16 Apr 2011 18:24:38 +0000 (18:24 +0000)
Change centrality option from 5 to 20 when centrality is given up to 20
Photon and Pi0: Reduce number of filled histograms by default. Removed histograms can still be recovered via setters (conversion, pi0 invariant mass for different module combinations, event multiplitcity histograms for pp ...)
Macro: Correct some settings (phi range, z vertex event binning), reflect the above changes.

PWG4/PartCorrBase/AliCaloTrackReader.cxx
PWG4/PartCorrBase/AliCaloTrackReader.h
PWG4/PartCorrDep/AliAnaPhoton.cxx
PWG4/PartCorrDep/AliAnaPhoton.h
PWG4/PartCorrDep/AliAnaPi0.cxx
PWG4/PartCorrDep/AliAnaPi0.h
PWG4/macros/AddTaskPi0.C

index 139e641f3a5dd24ad4855182f884de1c2e85039c..56f7b644e52e197c2b5f76a9cd6ae6c735e0c411 100755 (executable)
@@ -77,7 +77,7 @@ ClassImp(AliCaloTrackReader)
     fMixedEvent(NULL), fNMixedEvent(1), fVertex(NULL), 
     fWriteOutputDeltaAOD(kFALSE),fOldAOD(kFALSE),fCaloFilterPatch(kFALSE),
     fEMCALClustersListName(""),fZvtxCut(0.), 
-    fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE),
+    fDoEventSelection(kFALSE),   fDoV0ANDEventSelection(kFALSE), fUseEventsWithPrimaryVertex(kFALSE),
     fTriggerAnalysis (new AliTriggerAnalysis), fCentralityClass("V0M"),fCentralityOpt(10)
    
 {
@@ -479,10 +479,9 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * /*cur
           bV0AND = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0AND);
         //else bV0AND = //FIXME FOR AODs
         if(!bV0AND) return kFALSE;
-        
       }
       
-      if(!CheckForPrimaryVertex()) return kFALSE;
+      if(fUseEventsWithPrimaryVertex && !CheckForPrimaryVertex()) return kFALSE;
       
     }//CaloFilter patch
     else{ 
@@ -494,7 +493,7 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * /*cur
           if(bPileup) return kFALSE;
           
           Bool_t bGoodV = selection[1]; 
-          if(!bGoodV) return kFALSE;
+          if(fUseEventsWithPrimaryVertex && !bGoodV) return kFALSE;
           
           if(fDoV0ANDEventSelection){
             Bool_t bV0AND = selection[2]; 
@@ -512,7 +511,7 @@ Bool_t AliCaloTrackReader::FillInputEvent(const Int_t iEntry, const char * /*cur
       else {
         //printf("AliCaloTrackReader::FillInputEvent() - No clusters in event\n");
         //Remove events with  vertex (0,0,0), bad vertex reconstruction
-        if(TMath::Abs(fVertex[0][0]) < 1.e-6 && TMath::Abs(fVertex[0][1]) < 1.e-6 && TMath::Abs(fVertex[0][2]) < 1.e-6) return kFALSE;
+        if(fUseEventsWithPrimaryVertex && TMath::Abs(fVertex[0][0]) < 1.e-6 && TMath::Abs(fVertex[0][1]) < 1.e-6 && TMath::Abs(fVertex[0][2]) < 1.e-6) return kFALSE;
         
         //First filtered AODs, track multiplicity stored there.  
         fTrackMult = (Int_t) ((AliAODHeader*)fInputEvent->GetHeader())->GetCentrality();
@@ -597,11 +596,11 @@ Int_t AliCaloTrackReader::GetEventCentrality() const {
   //Return current event centrality
   
   if(GetCentrality()){
-    if(fCentralityOpt==100)     return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass);
-    else if(fCentralityOpt==10) return GetCentrality()->GetCentralityClass10(fCentralityClass); 
-    else if(fCentralityOpt==5)  return GetCentrality()->GetCentralityClass5(fCentralityClass);
+    if(fCentralityOpt==100)     return (Int_t) GetCentrality()->GetCentralityPercentile(fCentralityClass); // 100 bins max
+    else if(fCentralityOpt==10) return GetCentrality()->GetCentralityClass10(fCentralityClass);// 10 bins max
+    else if(fCentralityOpt==20) return GetCentrality()->GetCentralityClass5(fCentralityClass); // 20 bins max
     else {
-      printf("AliAnaPartCorrBaseClass::Unknown centrality option %d, use 5, 10 or 100\n",fCentralityOpt);
+      printf("AliAnaPartCorrBaseClass::Unknown centrality option %d, use 10, 20 or 100\n",fCentralityOpt);
       return 0;
     } 
   }
@@ -1108,7 +1107,7 @@ Bool_t AliCaloTrackReader::CheckForPrimaryVertex(){
   //Only for ESDs ...
   
   AliESDEvent * event = dynamic_cast<AliESDEvent*> (fInputEvent);
-  if(!event) return kFALSE;
+  if(!event) return kTRUE;
   
   if(event->GetPrimaryVertexTracks()->GetNContributors() > 0) {
     return kTRUE;
index 63b19fe8ab3b719eb9e54d78650fe152b458ac3f..e4df4ef469d697f79407220d446ec2a6be3bbf60 100755 (executable)
@@ -201,6 +201,12 @@ public:
   void    SwitchOnV0ANDSelection()                { fDoV0ANDEventSelection = kTRUE  ; }
   void    SwitchOffV0ANDSelection()               { fDoV0ANDEventSelection = kFALSE ; }
   Bool_t  IsV0ANDEventSelectionDone()       const { return fDoV0ANDEventSelection   ; } 
+
+  void    SwitchOnPrimaryVertexSelection()        { fUseEventsWithPrimaryVertex = kTRUE  ; }
+  void    SwitchOffPrimaryVertexSelection()       { fUseEventsWithPrimaryVertex = kFALSE ; }
+  Bool_t  IsPrimaryVertexSelectionDone()    const { return fUseEventsWithPrimaryVertex   ; } 
+  
+  
   
   // Track selection
   ULong_t GetTrackStatus()                  const {return fTrackStatus      ; }
@@ -391,19 +397,20 @@ public:
   Int_t            fV0ADC[2]    ;        // Integrated V0 signal
   Int_t            fV0Mul[2]    ;        // Integrated V0 Multiplicity
 
-  Bool_t           fCaloFilterPatch;     // CaloFilter patch
-  TString          fEMCALClustersListName; //Alternative list of clusters produced elsewhere and not from InputEvent
-  Float_t          fZvtxCut ;             // Cut on vertex position  
-  Bool_t           fDoEventSelection;    // Select events depending on V0, pileup, vertex well reconstructed, at least 1 track ...
-  Bool_t           fDoV0ANDEventSelection; // Select events depending on V0, fDoEventSelection should be on
-  AliTriggerAnalysis* fTriggerAnalysis;  // Access to trigger selection algorithm for V0AND calculation
+  Bool_t           fCaloFilterPatch;             // CaloFilter patch
+  TString          fEMCALClustersListName;       // Alternative list of clusters produced elsewhere and not from InputEvent
+  Float_t          fZvtxCut ;                     // Cut on vertex position  
+  Bool_t           fDoEventSelection;            // Select events depending on V0, pileup, vertex well reconstructed, at least 1 track ...
+  Bool_t           fDoV0ANDEventSelection;       // Select events depending on V0, fDoEventSelection should be on
+  Bool_t           fUseEventsWithPrimaryVertex ; // Select events with primary vertex
+  AliTriggerAnalysis* fTriggerAnalysis;          // Access to trigger selection algorithm for V0AND calculation
   
   //Centrality
   TString          fCentralityClass;     // Name of selected centrality class     
   Int_t            fCentralityOpt;       // Option for the returned value of the centrality, possible options 5, 10, 100
   Int_t            fCentralityBin[2];    // Minimum and maximum value of the centrality for the analysis
   
-  ClassDef(AliCaloTrackReader,28)
+  ClassDef(AliCaloTrackReader,29)
 } ;
 
 
index 39a7982c757384149cbcac1cd1e2cd102aab45b8..9aa529d63dd82255dae6cf477d4f62d23150202a 100755 (executable)
@@ -54,11 +54,11 @@ ClassImp(AliAnaPhoton)
     AliAnaPartCorrBaseClass(), fCalorimeter(""), 
     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
     fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
-    fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
+    fCheckConversion(kFALSE), fRemoveConvertedPair(kFALSE), fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
     fConvAsymCut(1.), fConvDEtaCut(2.),fConvDPhiMinCut(-1.), fConvDPhiMaxCut(7.), 
     //fhVertex(0), 
     fhNtraNclu(0), fhNCellsPt(0),
-    fhPtPhoton(0),     fhPhiPhoton(0),       fhEtaPhoton(0),       fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
+    fhEPhoton(0),      fhPtPhoton(0),  fhPhiPhoton(0),  fhEtaPhoton(0),  fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
     fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
     fhConvDeltaEta(0), fhConvDeltaPhi(0),    fhConvDeltaEtaPhi(0), fhConvAsym(0),     fhConvPt(0),
     //MC
@@ -163,19 +163,24 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
   fhNCellsPt->SetYTitle("# of cells in cluster");
   outputContainer->Add(fhNCellsPt);  
   
-  fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+  fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax); 
+  fhEPhoton->SetYTitle("N");
+  fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
+  outputContainer->Add(fhEPhoton) ;   
+  
+  fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax); 
   fhPtPhoton->SetYTitle("N");
   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
   outputContainer->Add(fhPtPhoton) ; 
   
   fhPhiPhoton  = new TH2F
-    ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
   fhPhiPhoton->SetYTitle("#phi (rad)");
   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhPhiPhoton) ; 
   
   fhEtaPhoton  = new TH2F
-    ("hEtaPhoton","#eta_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
   fhEtaPhoton->SetYTitle("#eta");
   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhEtaPhoton) ;
@@ -185,62 +190,64 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
   fhEtaPhiPhoton->SetXTitle("#eta");
   outputContainer->Add(fhEtaPhiPhoton) ;
-  
-  fhEtaPhi05Photon  = new TH2F
-  ("hEtaPhi05Photon","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
-  fhEtaPhi05Photon->SetYTitle("#phi (rad)");
-  fhEtaPhi05Photon->SetXTitle("#eta");
-  outputContainer->Add(fhEtaPhi05Photon) ;
-  
+  if(GetMinPt() < 0.5){
+    fhEtaPhi05Photon  = new TH2F
+    ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax); 
+    fhEtaPhi05Photon->SetYTitle("#phi (rad)");
+    fhEtaPhi05Photon->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhi05Photon) ;
+  }
   
   //Conversion
-
-  fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax); 
-  fhPtPhotonConv->SetYTitle("N");
-  fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
-  outputContainer->Add(fhPtPhotonConv) ; 
-  
-  fhEtaPhiPhotonConv  = new TH2F
-  ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
-  fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
-  fhEtaPhiPhotonConv->SetXTitle("#eta");
-  outputContainer->Add(fhEtaPhiPhotonConv) ;
-  
-  fhEtaPhi05PhotonConv  = new TH2F
-  ("hEtaPhi05PhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
-  fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
-  fhEtaPhi05PhotonConv->SetXTitle("#eta");
-  outputContainer->Add(fhEtaPhi05PhotonConv) ;
-  
-  fhConvDeltaEta  = new TH2F
-  ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5); 
-  fhConvDeltaEta->SetYTitle("#Delta #eta");
-  fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
-  outputContainer->Add(fhConvDeltaEta) ;
-  
-  fhConvDeltaPhi  = new TH2F
-  ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5); 
-  fhConvDeltaPhi->SetYTitle("#Delta #phi");
-  fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
-  outputContainer->Add(fhConvDeltaPhi) ;
-  
-  fhConvDeltaEtaPhi  = new TH2F
-  ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
-  fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
-  fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
-  outputContainer->Add(fhConvDeltaEtaPhi) ;
-  
-  fhConvAsym  = new TH2F
-  ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1); 
-  fhConvAsym->SetYTitle("Asymmetry");
-  fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
-  outputContainer->Add(fhConvAsym) ;  
-  
-  fhConvPt  = new TH2F
-  ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.); 
-  fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
-  fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
-  outputContainer->Add(fhConvPt) ;
+  if(fCheckConversion){
+    fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax); 
+    fhPtPhotonConv->SetYTitle("N");
+    fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtPhotonConv) ; 
+    
+    fhEtaPhiPhotonConv  = new TH2F
+    ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+    fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
+    fhEtaPhiPhotonConv->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhiPhotonConv) ;
+    if(GetMinPt() < 0.5){
+      fhEtaPhi05PhotonConv  = new TH2F
+      ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax); 
+      fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
+      fhEtaPhi05PhotonConv->SetXTitle("#eta");
+      outputContainer->Add(fhEtaPhi05PhotonConv) ;
+    }
+    
+    fhConvDeltaEta  = new TH2F
+    ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5); 
+    fhConvDeltaEta->SetYTitle("#Delta #eta");
+    fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaEta) ;
+    
+    fhConvDeltaPhi  = new TH2F
+    ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5); 
+    fhConvDeltaPhi->SetYTitle("#Delta #phi");
+    fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaPhi) ;
+    
+    fhConvDeltaEtaPhi  = new TH2F
+    ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+    fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
+    fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
+    outputContainer->Add(fhConvDeltaEtaPhi) ;
+    
+    fhConvAsym  = new TH2F
+    ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1); 
+    fhConvAsym->SetYTitle("Asymmetry");
+    fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvAsym) ;  
+    
+    fhConvPt  = new TH2F
+    ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.); 
+    fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
+    fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvPt) ;
+  }
   
   if(IsDataMC()){
     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50); 
@@ -451,195 +458,196 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
     outputContainer->Add(fhEtaUnknown) ;
        
-    
-    fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
-    fhPtConversionTagged->SetYTitle("N");
-    fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
-    outputContainer->Add(fhPtConversionTagged) ; 
-    
-    fhPtAntiNeutronTagged  = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
-    fhPtAntiNeutronTagged->SetYTitle("N");
-    fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
-    outputContainer->Add(fhPtAntiNeutronTagged) ; 
-    
-    fhPtAntiProtonTagged  = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
-    fhPtAntiProtonTagged->SetYTitle("N");
-    fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
-    outputContainer->Add(fhPtAntiProtonTagged) ; 
-
-    fhPtUnknownTagged  = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
-    fhPtUnknownTagged->SetYTitle("N");
-    fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
-    outputContainer->Add(fhPtUnknownTagged) ;     
-    
-    fhConvDeltaEtaMCConversion  = new TH2F
-    ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5); 
-    fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
-    fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaEtaMCConversion) ;
-    
-    fhConvDeltaPhiMCConversion  = new TH2F
-    ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5); 
-    fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
-    fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaPhiMCConversion) ;
-    
-    fhConvDeltaEtaPhiMCConversion  = new TH2F
-    ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
-    fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
-    fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
-    outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
-    
-    fhConvAsymMCConversion  = new TH2F
-    ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1); 
-    fhConvAsymMCConversion->SetYTitle("Asymmetry");
-    fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvAsymMCConversion) ;
-    
-    fhConvPtMCConversion  = new TH2F
-    ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.); 
-    fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
-    fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvPtMCConversion) ;    
-    
-    fhConvDispersionMCConversion  = new TH2F
-    ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.); 
-    fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
-    fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
-    outputContainer->Add(fhConvDispersionMCConversion) ;   
-    
-    fhConvM02MCConversion  = new TH2F
-    ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
-    fhConvM02MCConversion->SetYTitle("M02 cluster 1");
-    fhConvM02MCConversion->SetXTitle("M02 cluster 2");
-    outputContainer->Add(fhConvM02MCConversion) ;           
-    
-    fhConvDeltaEtaMCAntiNeutron  = new TH2F
-    ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5); 
-    fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
-    fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
-    
-    fhConvDeltaPhiMCAntiNeutron  = new TH2F
-    ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5); 
-    fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
-    fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
-    
-    fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
-    ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
-    fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
-    fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
-    outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;    
-    
-    fhConvAsymMCAntiNeutron  = new TH2F
-    ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1); 
-    fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
-    fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvAsymMCAntiNeutron) ;
-    
-    fhConvPtMCAntiNeutron  = new TH2F
-    ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.); 
-    fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
-    fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvPtMCAntiNeutron) ;    
-    
-    fhConvDispersionMCAntiNeutron  = new TH2F
-    ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.); 
-    fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
-    fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
-    outputContainer->Add(fhConvDispersionMCAntiNeutron) ;       
-    
-    fhConvM02MCAntiNeutron  = new TH2F
-    ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
-    fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
-    fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
-    outputContainer->Add(fhConvM02MCAntiNeutron) ;  
-    
-    fhConvDeltaEtaMCAntiProton  = new TH2F
-    ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5); 
-    fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
-    fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
-    
-    fhConvDeltaPhiMCAntiProton  = new TH2F
-    ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5); 
-    fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
-    fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
-    
-    fhConvDeltaEtaPhiMCAntiProton  = new TH2F
-    ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
-    fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
-    fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
-    outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;    
-    
-    fhConvAsymMCAntiProton  = new TH2F
-    ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1); 
-    fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
-    fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvAsymMCAntiProton) ;
-    
-    fhConvPtMCAntiProton  = new TH2F
-    ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.); 
-    fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
-    fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvPtMCAntiProton) ;
-    
-    fhConvDispersionMCAntiProton  = new TH2F
-    ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.); 
-    fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
-    fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
-    outputContainer->Add(fhConvDispersionMCAntiProton) ;       
-    
-    fhConvM02MCAntiProton  = new TH2F
-    ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
-    fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
-    fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
-    outputContainer->Add(fhConvM02MCAntiProton) ;       
-    
-    fhConvDeltaEtaMCString  = new TH2F
-    ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5); 
-    fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
-    fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaEtaMCString) ;
-    
-    fhConvDeltaPhiMCString  = new TH2F
-    ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5); 
-    fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
-    fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvDeltaPhiMCString) ;
-    
-    fhConvDeltaEtaPhiMCString  = new TH2F
-    ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
-    fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
-    fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
-    outputContainer->Add(fhConvDeltaEtaPhiMCString) ;    
-    
-    fhConvAsymMCString  = new TH2F
-    ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1); 
-    fhConvAsymMCString->SetYTitle("Asymmetry");
-    fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvAsymMCString) ;
-    
-    fhConvPtMCString  = new TH2F
-    ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.); 
-    fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
-    fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
-    outputContainer->Add(fhConvPtMCString) ;
-    
-    fhConvDispersionMCString  = new TH2F
-    ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
-    fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
-    fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
-    outputContainer->Add(fhConvDispersionMCString) ;       
-    
-    fhConvM02MCString  = new TH2F
-    ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
-    fhConvM02MCString->SetYTitle("M02 cluster 1");
-    fhConvM02MCString->SetXTitle("M02 cluster 2");
-    outputContainer->Add(fhConvM02MCString) ;       
-    
+    if(fCheckConversion){  
+      fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
+      fhPtConversionTagged->SetYTitle("N");
+      fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtConversionTagged) ; 
+      
+      
+      fhPtAntiNeutronTagged  = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
+      fhPtAntiNeutronTagged->SetYTitle("N");
+      fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtAntiNeutronTagged) ; 
+      
+      fhPtAntiProtonTagged  = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
+      fhPtAntiProtonTagged->SetYTitle("N");
+      fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtAntiProtonTagged) ; 
+      
+      fhPtUnknownTagged  = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
+      fhPtUnknownTagged->SetYTitle("N");
+      fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+      outputContainer->Add(fhPtUnknownTagged) ;     
+      
+      fhConvDeltaEtaMCConversion  = new TH2F
+      ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5); 
+      fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
+      fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaEtaMCConversion) ;
+      
+      fhConvDeltaPhiMCConversion  = new TH2F
+      ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5); 
+      fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
+      fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaPhiMCConversion) ;
+      
+      fhConvDeltaEtaPhiMCConversion  = new TH2F
+      ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+      fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
+      fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
+      outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
+      
+      fhConvAsymMCConversion  = new TH2F
+      ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1); 
+      fhConvAsymMCConversion->SetYTitle("Asymmetry");
+      fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvAsymMCConversion) ;
+      
+      fhConvPtMCConversion  = new TH2F
+      ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.); 
+      fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
+      fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvPtMCConversion) ;    
+      
+      fhConvDispersionMCConversion  = new TH2F
+      ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.); 
+      fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
+      fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
+      outputContainer->Add(fhConvDispersionMCConversion) ;   
+      
+      fhConvM02MCConversion  = new TH2F
+      ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+      fhConvM02MCConversion->SetYTitle("M02 cluster 1");
+      fhConvM02MCConversion->SetXTitle("M02 cluster 2");
+      outputContainer->Add(fhConvM02MCConversion) ;           
+      
+      fhConvDeltaEtaMCAntiNeutron  = new TH2F
+      ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5); 
+      fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
+      fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
+      
+      fhConvDeltaPhiMCAntiNeutron  = new TH2F
+      ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5); 
+      fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+      fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
+      
+      fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
+      ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+      fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+      fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
+      outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;    
+      
+      fhConvAsymMCAntiNeutron  = new TH2F
+      ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1); 
+      fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
+      fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvAsymMCAntiNeutron) ;
+      
+      fhConvPtMCAntiNeutron  = new TH2F
+      ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.); 
+      fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
+      fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvPtMCAntiNeutron) ;    
+      
+      fhConvDispersionMCAntiNeutron  = new TH2F
+      ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.); 
+      fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
+      fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
+      outputContainer->Add(fhConvDispersionMCAntiNeutron) ;       
+      
+      fhConvM02MCAntiNeutron  = new TH2F
+      ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+      fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
+      fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
+      outputContainer->Add(fhConvM02MCAntiNeutron) ;  
+      
+      fhConvDeltaEtaMCAntiProton  = new TH2F
+      ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5); 
+      fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
+      fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
+      
+      fhConvDeltaPhiMCAntiProton  = new TH2F
+      ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5); 
+      fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
+      fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
+      
+      fhConvDeltaEtaPhiMCAntiProton  = new TH2F
+      ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+      fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
+      fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
+      outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;    
+      
+      fhConvAsymMCAntiProton  = new TH2F
+      ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1); 
+      fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
+      fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvAsymMCAntiProton) ;
+      
+      fhConvPtMCAntiProton  = new TH2F
+      ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.); 
+      fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
+      fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvPtMCAntiProton) ;
+      
+      fhConvDispersionMCAntiProton  = new TH2F
+      ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.); 
+      fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
+      fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
+      outputContainer->Add(fhConvDispersionMCAntiProton) ;       
+      
+      fhConvM02MCAntiProton  = new TH2F
+      ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+      fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
+      fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
+      outputContainer->Add(fhConvM02MCAntiProton) ;       
+      
+      fhConvDeltaEtaMCString  = new TH2F
+      ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5); 
+      fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
+      fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaEtaMCString) ;
+      
+      fhConvDeltaPhiMCString  = new TH2F
+      ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5); 
+      fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
+      fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvDeltaPhiMCString) ;
+      
+      fhConvDeltaEtaPhiMCString  = new TH2F
+      ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+      fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
+      fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
+      outputContainer->Add(fhConvDeltaEtaPhiMCString) ;    
+      
+      fhConvAsymMCString  = new TH2F
+      ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1); 
+      fhConvAsymMCString->SetYTitle("Asymmetry");
+      fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvAsymMCString) ;
+      
+      fhConvPtMCString  = new TH2F
+      ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.); 
+      fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
+      fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
+      outputContainer->Add(fhConvPtMCString) ;
+      
+      fhConvDispersionMCString  = new TH2F
+      ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+      fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
+      fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
+      outputContainer->Add(fhConvDispersionMCString) ;       
+      
+      fhConvM02MCString  = new TH2F
+      ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+      fhConvM02MCString->SetYTitle("M02 cluster 1");
+      fhConvM02MCString->SetXTitle("M02 cluster 2");
+      outputContainer->Add(fhConvM02MCString) ;       
+    }
     
   }//Histos with MC
     
@@ -684,6 +692,7 @@ void AliAnaPhoton::InitParameters()
        
   fRejectTrackMatch       = kTRUE ;
   fCheckConversion        = kFALSE;
+  fRemoveConvertedPair    = kFALSE;
   fAddConvertedPairsToAOD = kFALSE;
        
 }
@@ -713,10 +722,13 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
   TLorentzVector mom, mom2 ;
   Int_t nCaloClusters = pl->GetEntriesFast();
   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
-  Bool_t * indexConverted = new Bool_t[nCaloClusters];
-  for (Int_t i = 0; i < nCaloClusters; i++) 
-    indexConverted[i] = kFALSE;
-       
+  Bool_t * indexConverted = 0x0;
+  if(fCheckConversion){
+    indexConverted = new Bool_t[nCaloClusters];
+    for (Int_t i = 0; i < nCaloClusters; i++) 
+      indexConverted[i] = kFALSE;
+       }
+  
   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
 
   //----------------------------------------------------
@@ -898,12 +910,12 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
     //--------------------------------------------------------------------------------------
 
     // Do analysis only if there are more than one cluster
-    if( nCaloClusters > 1){
+    if( nCaloClusters > 1 && fCheckConversion){
       Bool_t bConverted = kFALSE;
       Int_t id2 = -1;
                  
       //Check if set previously as converted couple, if so skip its use.
-      if (fCheckConversion && indexConverted[icalo]) continue;
+      if (indexConverted[icalo]) continue;
                  
       // Second cluster loop
       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
@@ -1064,7 +1076,7 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
         }
         
         //Do not add the current calocluster
-        if(fCheckConversion) continue;
+        if(fRemoveConvertedPair) continue;
         else {
           //printf("TAGGED\n");
           //Tag this cluster as likely conversion
@@ -1161,19 +1173,19 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
          Float_t etacluster = ph->Eta();
          Float_t ecluster   = ph->E();
          
+    fhEPhoton   ->Fill(ecluster);
          fhPtPhoton  ->Fill(ptcluster);
-    if(ph->IsTagged())fhPtPhotonConv->Fill(ptcluster);
          fhPhiPhoton ->Fill(ptcluster,phicluster);
-         fhEtaPhoton ->Fill(ptcluster,etacluster);
-    if(ptcluster > 0.5){
-      fhEtaPhiPhoton ->Fill(etacluster, phicluster);
-      if(ph->IsTagged())fhEtaPhiPhotonConv->Fill(etacluster, phicluster);
-    }
-    else {
-      fhEtaPhi05Photon ->Fill(etacluster, phicluster);
-      if(ph->IsTagged())fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
+         fhEtaPhoton ->Fill(ptcluster,etacluster);    
+    if(ecluster > 0.5)         fhEtaPhiPhoton ->Fill(etacluster, phicluster);
+    else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
+    
+    if(fCheckConversion &&ph->IsTagged()){
+      fhPtPhotonConv->Fill(ptcluster);
+      if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
+      else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
     }
-
+    
     //.......................................
          //Play with the MC data if available
          if(IsDataMC()){
@@ -1231,7 +1243,7 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
         fhPtAntiNeutron  ->Fill(ptcluster);
         fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
         fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
-        if(ph->IsTagged()) fhPtAntiNeutronTagged ->Fill(ptcluster);
+        if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
 
       }
       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
@@ -1239,14 +1251,14 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
         fhPtAntiProton  ->Fill(ptcluster);
         fhPhiAntiProton ->Fill(ptcluster,phicluster);
         fhEtaAntiProton ->Fill(ptcluster,etacluster);
-        if(ph->IsTagged()) fhPtAntiProtonTagged ->Fill(ptcluster);
+        if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
 
       }      
            else{
              fhPtUnknown  ->Fill(ptcluster);
              fhPhiUnknown ->Fill(ptcluster,phicluster);
              fhEtaUnknown ->Fill(ptcluster,etacluster);
-        if(ph->IsTagged()) fhPtUnknownTagged ->Fill(ptcluster);
+        if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
 
              
         //              printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
index ecece04ce742692f0cfebe7a9ec56f07b63ae236..65c941dcadceed75bd38a286d0c01be83a632e96 100755 (executable)
@@ -86,17 +86,21 @@ class AliAnaPhoton : public AliAnaPartCorrBaseClass {
 
   // ** Conversion pair analysis **
   
-  Float_t  GetMassCut()                     const { return fMassCut ; }
-  void     SetMassCut(Float_t m)                  { fMassCut    = m ; }
+  Float_t  GetMassCut()                     const { return fMassCut           ; }
+  void     SetMassCut(Float_t m)                  { fMassCut    = m           ; }
   
   Bool_t   IsCheckConversionOn()            const { return fCheckConversion   ; }
   void     SwitchOnConversionChecker()            { fCheckConversion = kTRUE  ; }
   void     SwitchOffConversionChecker()           { fCheckConversion = kFALSE ; }  
        
   Bool_t   AreConvertedPairsInAOD()         const { return fAddConvertedPairsToAOD   ; }
-  void     SwitchOnAdditionConvertedPairsToAOD()  { fAddConvertedPairsToAOD = kTRUE  ; }
+  void     SwitchOnAdditionConvertedPairsToAOD()  { fAddConvertedPairsToAOD = kTRUE  ; fCheckConversion = kTRUE ; }
   void     SwitchOffAdditionConvertedPairsToAOD() { fAddConvertedPairsToAOD = kFALSE ; }  
        
+  Bool_t   AreConvertedPairsRemoved()       const { return fRemoveConvertedPair      ; }
+  void     SwitchOnConvertedPairsRemoval()        { fRemoveConvertedPair  = kTRUE    ; fCheckConversion = kTRUE ; }
+  void     SwitchOffConvertedPairsRemoval()       { fRemoveConvertedPair  = kFALSE   ; }    
+  
   void     SetConvAsymCut(Float_t c)              { fConvAsymCut = c    ; }
   Float_t  GetConvAsymCut()                 const { return fConvAsymCut ; }
   
@@ -121,6 +125,7 @@ class AliAnaPhoton : public AliAnaPartCorrBaseClass {
   
   //Conversion pairs selection cuts
   Bool_t   fCheckConversion;             // Combine pairs of clusters with mass close to 0
+  Bool_t   fRemoveConvertedPair;         // Combine pairs of clusters with mass close to 0
   Bool_t   fAddConvertedPairsToAOD;      // Put Converted pairs in AOD
   Float_t  fMassCut;                     // Mass cut for the conversion pairs selection  
   Float_t  fConvAsymCut;                 // Select conversion pairs when asymmetry is smaller than cut
@@ -131,6 +136,7 @@ class AliAnaPhoton : public AliAnaPartCorrBaseClass {
   //Histograms 
   TH2F * fhNtraNclu;                     //! track multiplicity distribution vs cluster multiplicity
   TH2F * fhNCellsPt;                     //! number of cells in cluster vs pt 
+  TH1F * fhEPhoton    ;                  //! Number of identified photon vs energy
   TH1F * fhPtPhoton   ;                  //! Number of identified photon vs transerse momentum 
   TH2F * fhPhiPhoton  ;                  //! Azimuthal angle of identified  photon vs transerse momentum 
   TH2F * fhEtaPhoton  ;                  //! Pseudorapidity of identified  photon vs transerse momentum 
@@ -237,7 +243,7 @@ class AliAnaPhoton : public AliAnaPartCorrBaseClass {
   TH2F * fhConvDispersionMCString;       //! Small mass cluster pairs, dispersion of cluster 1 vs cluster 2, origin of both clusters is string
   TH2F * fhConvM02MCString;              //! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is string
 
-   ClassDef(AliAnaPhoton,10)
+   ClassDef(AliAnaPhoton,11)
 
 } ;
  
index 359b5cc1d397bdfe3d5e0e5cf279f8fbf80b8419..43b90775dc972117018c6a203f8c3e53ec86a94d 100755 (executable)
@@ -58,8 +58,10 @@ ClassImp(AliAnaPi0)
 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
 fNmaxMixEv(0), fCalorimeter(""),
-fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
-fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),  fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
+fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), 
+fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
+fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),  
+fMakeInvPtPlots(kFALSE), fSameSM(kFALSE), fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
 fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverCellEBins(kFALSE), fUseAverClusterEBins(kFALSE),
 fUseAverClusterEDenBins(0), //fUseAverClusterPairRBins(0), fUseAverClusterPairRWeightBins(0), fUseEMaxBins(0),
 fFillBadDistHisto(kFALSE),
@@ -75,7 +77,7 @@ fhRe1(0x0),      fhMi1(0x0),       fhRe2(0x0),      fhMi2(0x0),      fhRe3(0x0),
 fhReInvPt1(0x0), fhMiInvPt1(0x0),  fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
 fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
 fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),  
-fhEvents(0x0),   fhCentrality(0x0),
+fhEvents(0x0),   fhCentrality(0x0),fhCentralityNoPair(0x0),
 fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
 fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
 fhPrimPi0OpeningAngle(0x0), fhPrimPi0CosOpeningAngle(0x0),
@@ -131,19 +133,19 @@ void AliAnaPi0::InitParameters()
 
   fMultiCutAna = kFALSE;
   
-  fNPtCuts = 3;
+  fNPtCuts = 1;
   fPtCuts[0] = 0.; fPtCuts[1] = 0.3;   fPtCuts[2] = 0.5;
   for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
   
-  fNAsymCuts = 4;
-  fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.8; fAsymCuts[2] = 0.6;   fAsymCuts[3] = 0.1;    
+  fNAsymCuts = 2;
+  fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; //  fAsymCuts[3] = 0.1;    
   for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
 
-  fNCellNCuts = 3;
+  fNCellNCuts = 1;
   fCellNCuts[0] = 0; fCellNCuts[1] = 1;   fCellNCuts[2] = 2;   
   for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i]  = 0;
 
-  fNPIDBits = 2;
+  fNPIDBits = 1;
   fPIDBits[0] = 0;   fPIDBits[1] = 2; //  fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut,  dispersion, neutral, dispersion&&neutral
   for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
 
@@ -167,8 +169,8 @@ TObjString * AliAnaPi0::GetAnalysisCuts()
   parList+=onePar ;
   snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
   parList+=onePar ;
-  snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
-           fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
+  snprintf(onePar,buffersize,"Pair in same Module: %d ; Fill Different SM histos %d; CheckConversions %d; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
+           fSameSM, fFillSMCombinations, fCheckConversion, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
   parList+=onePar ;
   snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
   parList+=onePar ;
@@ -277,98 +279,103 @@ TList * AliAnaPi0::GetCreateOutputObjects()
   Int_t ntrmmax   = GetHistoTrackMultiplicityMax();
   Int_t ntrmmin   = GetHistoTrackMultiplicityMin(); 
   
-  fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
-  fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
-  outputContainer->Add(fhAverTotECluster) ;
-  
-  fhAverTotECell    = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
-  fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
-  outputContainer->Add(fhAverTotECell) ;
-  
-  fhAverTotECellvsCluster    = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
-  fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
-  fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
-  outputContainer->Add(fhAverTotECellvsCluster) ;
-  
-  fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
-  fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-  outputContainer->Add(fhEDensityCluster) ;
-  
-  fhEDensityCell    = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
-  fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
-  outputContainer->Add(fhEDensityCell) ;
-  
-  fhEDensityCellvsCluster    = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
-  fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
-  fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-  outputContainer->Add(fhEDensityCellvsCluster) ;
-  
-//  fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
-//  fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
-//  outputContainer->Add(fhClusterPairDist) ;
-//  
-//  fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
-//  fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
-//  outputContainer->Add(fhClusterPairDistWeight) ;
-//   
-//  fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
-//  fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-//  outputContainer->Add(fhAverClusterPairDist) ;
-//  
-//  fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
-//  fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
-//  outputContainer->Add(fhAverClusterPairDistWeight) ;
-//  
-//  fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
-//  fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-//  fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-//  outputContainer->Add(fhAverClusterPairDistvsAverE) ;
-//  
-//  fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
-//  fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
-//  fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-//  outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
-  
-//  fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
-//  fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
-//  fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
-//  outputContainer->Add(fhAverClusterPairDistvsN) ;
-//  
-//  fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
-//  fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
-//  fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
-//  outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
-  
-//  fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
-//  fhMaxEvsClustMult->SetXTitle("E_{max}");
-//  fhMaxEvsClustMult->SetYTitle("N_{cluster}");
-//  outputContainer->Add(fhMaxEvsClustMult) ;
-//  
-//  fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
-//  fhMaxEvsClustEDen->SetXTitle("E_{max}");
-//  fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
-//  outputContainer->Add(fhMaxEvsClustEDen) ;
-  
-  fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-  fhReConv->SetXTitle("p_{T} (GeV/c)");
-  fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-  outputContainer->Add(fhReConv) ;
-  
-  fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-  fhReConv2->SetXTitle("p_{T} (GeV/c)");
-  fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-  outputContainer->Add(fhReConv2) ;
-  
-  if(fDoOwnMix){
-    fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-    fhMiConv->SetXTitle("p_{T} (GeV/c)");
-    fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-    outputContainer->Add(fhMiConv) ;
-  
-    fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-    fhMiConv2->SetXTitle("p_{T} (GeV/c)");
-    fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-    outputContainer->Add(fhMiConv2) ;
+  if(fNCentrBin > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins)){
+    
+    fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
+    fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
+    outputContainer->Add(fhAverTotECluster) ;
+    
+    fhAverTotECell    = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
+    fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
+    outputContainer->Add(fhAverTotECell) ;
+    
+    fhAverTotECellvsCluster    = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
+    fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
+    fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
+    outputContainer->Add(fhAverTotECellvsCluster) ;
+    
+    fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
+    fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+    outputContainer->Add(fhEDensityCluster) ;
+    
+    fhEDensityCell    = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
+    fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
+    outputContainer->Add(fhEDensityCell) ;
+    
+    fhEDensityCellvsCluster    = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
+    fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
+    fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+    outputContainer->Add(fhEDensityCellvsCluster) ;
+    
+    //  fhClusterPairDist = new TH1F("hClusterPairDist","Distance between clusters",250,0,750) ;
+    //  fhClusterPairDist->SetXTitle("#sqrt{(x_{1}-x_{2})^2+(z_{1}-z_{2})^2} (cm)");
+    //  outputContainer->Add(fhClusterPairDist) ;
+    //  
+    //  fhClusterPairDistWeight = new TH1F("hClusterPairDistWeighted","Distance between clusters, weighted by pair energy",200,0,400) ;
+    //  fhClusterPairDistWeight->SetXTitle("#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2}) (cm)");
+    //  outputContainer->Add(fhClusterPairDistWeight) ;
+    //   
+    //  fhAverClusterPairDist = new TH1F("hAverClusterPairDist","Average distance between clusters",250,0,750) ;
+    //  fhAverClusterPairDist->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+    //  outputContainer->Add(fhAverClusterPairDist) ;
+    //  
+    //  fhAverClusterPairDistWeight = new TH1F("hAverClusterPairDistWeighted","Average distance between clusters, weighted by pair energy",200,0,400) ;
+    //  fhAverClusterPairDistWeight->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
+    //  outputContainer->Add(fhAverClusterPairDistWeight) ;
+    //  
+    //  fhAverClusterPairDistvsAverE = new TH2F("hAverClusterPairDistvsAverE","Average distance between clusters",250,0,750,200,0,50) ;
+    //  fhAverClusterPairDistvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+    //  fhAverClusterPairDistvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+    //  outputContainer->Add(fhAverClusterPairDistvsAverE) ;
+    //  
+    //  fhAverClusterPairDistWeightvsAverE = new TH2F("hAverClusterPairDistWeightedvsAverE","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
+    //  fhAverClusterPairDistWeightvsAverE->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^2+(z_{1}E_{1}-z_{2}E_{2})^2}/ (E_{1}+E_{2})) / N_{pairs} (cm/GeV)");
+    //  fhAverClusterPairDistWeightvsAverE->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+    //  outputContainer->Add(fhAverClusterPairDistWeightvsAverE) ;
+    
+    //  fhAverClusterPairDistvsN = new TH2F("hAverClusterPairDistvsN","Average distance between clusters",250,0,750,200,0,50) ;
+    //  fhAverClusterPairDistvsN->SetXTitle("#Sigma (#sqrt{(x_{1}-x_{2})^{2}+(z_{1}-z_{2})^{2}}) / N_{pairs} (cm)");
+    //  fhAverClusterPairDistvsN->SetYTitle("N_{cluster}");
+    //  outputContainer->Add(fhAverClusterPairDistvsN) ;
+    //  
+    //  fhAverClusterPairDistWeightvsN = new TH2F("hAverClusterPairDistWeightedvsN","Average distance between clusters, weighted by pair energy",200,0,400,200,0,50) ;
+    //  fhAverClusterPairDistWeightvsN->SetXTitle("#Sigma (#sqrt{(x_{1}E_{1}-x_{2}E_{2})^{2}+(z_{1}E_{1}-z_{2}E_{2})^{2}}/ (E_{1}+E_{2})) / N_{pairs} (cm)");
+    //  fhAverClusterPairDistWeightvsN->SetYTitle("N_{cluster}");
+    //  outputContainer->Add(fhAverClusterPairDistWeightvsN) ;
+    
+    //  fhMaxEvsClustMult = new TH2F("hMaxEvsClustMult","",nptbins,ptmin,ptmax,50,0,50) ;
+    //  fhMaxEvsClustMult->SetXTitle("E_{max}");
+    //  fhMaxEvsClustMult->SetYTitle("N_{cluster}");
+    //  outputContainer->Add(fhMaxEvsClustMult) ;
+    //  
+    //  fhMaxEvsClustEDen = new TH2F("hMaxEvsClustEDen","",nptbins,ptmin,ptmax,200,0,50) ;
+    //  fhMaxEvsClustEDen->SetXTitle("E_{max}");
+    //  fhMaxEvsClustEDen->SetYTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
+    //  outputContainer->Add(fhMaxEvsClustEDen) ;
+  }//counting and average histograms
+  
+  if(fCheckConversion){
+    fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+    fhReConv->SetXTitle("p_{T} (GeV/c)");
+    fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+    outputContainer->Add(fhReConv) ;
+    
+    fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+    fhReConv2->SetXTitle("p_{T} (GeV/c)");
+    fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+    outputContainer->Add(fhReConv2) ;
+    
+    if(fDoOwnMix){
+      fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMiConv->SetXTitle("p_{T} (GeV/c)");
+      fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMiConv) ;
+      
+      fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMiConv2->SetXTitle("p_{T} (GeV/c)");
+      fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMiConv2) ;
+    }
   }
   
   for(Int_t ic=0; ic<fNCentrBin; ic++){
@@ -604,10 +611,16 @@ TList * AliAnaPi0::GetCreateOutputObjects()
   fhEvents->SetZTitle("RP bin");
   outputContainer->Add(fhEvents) ;
        
-  fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",fNCentrBin*10,0.,1.*fNCentrBin) ;
-  fhCentrality->SetXTitle("Centrality bin");
-  outputContainer->Add(fhCentrality) ;
-
+  if(fNCentrBin>1){
+    fhCentrality=new TH1D("hCentralityBin","Number of events in centrality bin",fNCentrBin,0.,1.*fNCentrBin) ;
+    fhCentrality->SetXTitle("Centrality bin");
+    outputContainer->Add(fhCentrality) ;
+    
+    fhCentralityNoPair=new TH1D("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",fNCentrBin,0.,1.*fNCentrBin) ;
+    fhCentralityNoPair->SetXTitle("Centrality bin");
+    outputContainer->Add(fhCentralityNoPair) ;
+  }
+  
   fhRealOpeningAngle  = new TH2D
   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi()); 
   fhRealOpeningAngle->SetYTitle("#theta(rad)");
@@ -873,78 +886,80 @@ TList * AliAnaPi0::GetCreateOutputObjects()
     }
   }
   
-  TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
-  for(Int_t imod=0; imod<fNModules; imod++){
-    //Module dependent invariant mass
-    snprintf(key, buffersize,"hReMod_%d",imod) ;
-    snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
-    fhReMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-    fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
-    fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-    outputContainer->Add(fhReMod[imod]) ;
-    if(fCalorimeter=="PHOS"){
-      snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
-      snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
-      fhReDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-      fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
-      fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-      outputContainer->Add(fhReDiffPHOSMod[imod]) ;
-    }
-    else{//EMCAL
-      if(imod<fNModules/2){
-        snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
-        snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
-        fhReSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
-        fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
-      }
-      if(imod<fNModules-2){
-        snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
-        snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
-        fhReSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
-        fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
-      }
-    }//EMCAL
-    
-    if(fDoOwnMix){ 
-      snprintf(key, buffersize,"hMiMod_%d",imod) ;
-      snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
-      fhMiMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-      fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
-      fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-      outputContainer->Add(fhMiMod[imod]) ;
-      
+  if(fFillSMCombinations){
+    TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
+    for(Int_t imod=0; imod<fNModules; imod++){
+      //Module dependent invariant mass
+      snprintf(key, buffersize,"hReMod_%d",imod) ;
+      snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
+      fhReMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
+      fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhReMod[imod]) ;
       if(fCalorimeter=="PHOS"){
-        snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
-        snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
-        fhMiDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
-        fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
-      }//PHOS
+        snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
+        snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
+        fhReDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+        fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
+        fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+        outputContainer->Add(fhReDiffPHOSMod[imod]) ;
+      }
       else{//EMCAL
         if(imod<fNModules/2){
-          snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
-          snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
-          fhMiSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
-          fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
+          snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
+          snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
+          fhReSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+          fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
         }
         if(imod<fNModules-2){
-          snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
-          snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
-          fhMiSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
-          fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
+          snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
+          snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
+          fhReSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+          fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
         }
-      }//EMCAL      
-     }
-  }
+      }//EMCAL
+      
+      if(fDoOwnMix){ 
+        snprintf(key, buffersize,"hMiMod_%d",imod) ;
+        snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
+        fhMiMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+        fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
+        fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+        outputContainer->Add(fhMiMod[imod]) ;
+        
+        if(fCalorimeter=="PHOS"){
+          snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
+          snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
+          fhMiDiffPHOSMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
+          fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
+        }//PHOS
+        else{//EMCAL
+          if(imod<fNModules/2){
+            snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
+            snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
+            fhMiSameSectorEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+            fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
+          }
+          if(imod<fNModules-2){
+            snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
+            snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
+            fhMiSameSideEMCALMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
+            fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
+          }
+        }//EMCAL      
+      }// own mix
+    }//loop combinations
+  } // SM combinations
       
 //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
 //  
@@ -1519,6 +1534,92 @@ void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1,  const Int_t i
     }
 }  
 
+//____________________________________________________________________________________________________________________________________________________
+void AliAnaPi0::CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) {
+// Count the number of clusters and cells, deposited energy, and do some averages in case multiplicity bins dependent on such numbers
+// are requested
+  if(fCalorimeter=="EMCAL"){ 
+    nClus = GetEMCALClusters()  ->GetEntriesFast();
+    nCell = GetEMCALCells()->GetNumberOfCells();
+    for(Int_t icl=0; icl < nClus; icl++) {
+      Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
+      eClusTot +=  e1;
+      //      if(e1 > emax) emax = e1;
+      //      ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
+      //      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
+      //        Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
+      //        ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
+      //        rtmp  =  TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
+      //        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
+      //        rxz  += rtmp;  
+      //        rxzw += rtmpw;
+      //        ncomb++;
+      //        fhClusterPairDist      ->Fill(rtmp);
+      //        fhClusterPairDistWeight->Fill(rtmpw);
+      //        //printf("Distance: %f; weighted  %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
+      //
+      //      }// second cluster loop
+    }// first cluster
+    
+    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetEMCALCells()->GetAmplitude(jce);
+      }
+  else {                     
+    nClus = GetPHOSClusters()->GetEntriesFast();
+    nCell = GetPHOSCells()   ->GetNumberOfCells();
+    for(Int_t icl=0; icl < nClus; icl++) {
+      Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
+      eClusTot +=  e1;
+      //      ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
+      //      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
+      //        Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
+      //        ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
+      //        rtmp  = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
+      //        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
+      //        rxz  += rtmp;  
+      //        rxzw += rtmpw;
+      //        ncomb++;
+      //        fhClusterPairDist      ->Fill(rtmp);
+      //        fhClusterPairDistWeight->Fill(rtmpw);
+      //      }// second cluster loop
+    }// first cluster
+    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetPHOSCells()->GetAmplitude(jce);
+      }
+  if(GetDebug() > 1) 
+    printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
+    
+    //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
+    eDenClus = eClusTot/nClus;
+    eDenCell = eCellTot/nCell;
+    fhEDensityCluster      ->Fill(eDenClus);
+    fhEDensityCell         ->Fill(eDenCell);
+    fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
+    //Fill the average number of cells or clusters per SM
+    eClusTot /=fNModules;
+    eCellTot /=fNModules;
+    fhAverTotECluster      ->Fill(eClusTot);
+    fhAverTotECell         ->Fill(eCellTot);
+    fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
+    //printf("Average Cluster: E %f, density %f;  Average Cell E %f, density  %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
+    
+    //  //Average weighted pair distance
+    //  rxz  /= ncomb;                             
+    //  rxzw /= ncomb;                             
+    //
+    //  fhAverClusterPairDist             ->Fill(rxz );
+    //  fhAverClusterPairDistWeight       ->Fill(rxzw);
+    //  fhAverClusterPairDistvsAverE      ->Fill(rxz ,eDenClus);
+    //  fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
+    //  fhAverClusterPairDistvsN          ->Fill(rxz ,nClus);
+    //  fhAverClusterPairDistWeightvsN    ->Fill(rxzw,nClus);
+    //  
+    //  //emax
+    //  fhMaxEvsClustEDen->Fill(emax,eDenClus);
+    //  fhMaxEvsClustMult->Fill(emax,nPhot);
+    
+    //printf("Average Distance: %f; weighted  %f\n ",rxz,rxzw);
+    
+}
+
 //____________________________________________________________________________________________________________________________________________________
 void AliAnaPi0::MakeAnalysisFillHistograms() 
 {
@@ -1549,94 +1650,23 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
 //  Float_t pos2[3];
 //  Float_t emax     = 0;
   
+  if(fNCentrBin > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins))
+    CountAndGetAverages(nClus,nCell,eClusTot,eCellTot,eDenClus,eDenCell);
+
+  
   if(GetDebug() > 1) 
     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
 
   //If less than photon 2 entries in the list, skip this event
-  if(nPhot < 2 ) return ; 
-
-  // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
-  // are requested
-  if(fCalorimeter=="EMCAL"){ 
-    nClus = GetEMCALClusters()  ->GetEntriesFast();
-    nCell = GetEMCALCells()->GetNumberOfCells();
-    for(Int_t icl=0; icl < nClus; icl++) {
-      Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
-      eClusTot +=  e1;
-//      if(e1 > emax) emax = e1;
-//      ((AliVCluster*)GetEMCALClusters()->At(icl))->GetPosition(pos1);
-//      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
-//        Float_t e2 = ((AliVCluster*)GetEMCALClusters()->At(icl2))->E();
-//        ((AliVCluster*)GetEMCALClusters()->At(icl2))->GetPosition(pos2);
-//        rtmp  =  TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
-//        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
-//        rxz  += rtmp;  
-//        rxzw += rtmpw;
-//        ncomb++;
-//        fhClusterPairDist      ->Fill(rtmp);
-//        fhClusterPairDistWeight->Fill(rtmpw);
-//        //printf("Distance: %f; weighted  %f\n ",rtmp,rtmp/(e1+((AliVCluster*)GetEMCALClusters()->At(icl2))->E()));
-//
-//      }// second cluster loop
-    }// first cluster
+  if(nPhot < 2 ) {
     
-    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetEMCALCells()->GetAmplitude(jce);
-  }
-  else {                     
-    nClus = GetPHOSClusters()->GetEntriesFast();
-    nCell = GetPHOSCells()   ->GetNumberOfCells();
-    for(Int_t icl=0; icl < nClus; icl++) {
-      Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
-      eClusTot +=  e1;
-//      ((AliVCluster*)GetPHOSClusters()->At(icl))->GetPosition(pos1);
-//      for(Int_t icl2=icl+1; icl2 < nClus; icl2++) {
-//        Float_t e2 = ((AliVCluster*)GetPHOSClusters()->At(icl2))->E();
-//        ((AliVCluster*)GetPHOSClusters()->At(icl2))->GetPosition(pos2);
-//        rtmp  = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0]) + (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
-//        rtmpw =  TMath::Sqrt((pos1[0]*e1-pos2[0]*e2)*(pos1[0]*e1-pos2[0]*e2) + (pos1[2]*e1-pos2[2]*e2)*(pos1[2]*e1-pos2[2]*e2))/(e1+e2);
-//        rxz  += rtmp;  
-//        rxzw += rtmpw;
-//        ncomb++;
-//        fhClusterPairDist      ->Fill(rtmp);
-//        fhClusterPairDistWeight->Fill(rtmpw);
-//      }// second cluster loop
-    }// first cluster
-    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetPHOSCells()->GetAmplitude(jce);
+    if(GetDebug() > 2)
+      printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
+    
+    if(fNCentrBin > 1) fhCentralityNoPair->Fill(GetEventCentrality() * fNCentrBin / GetReader()->GetCentralityOpt());
+    
+    return ;
   }
-  if(GetDebug() > 1) 
-    printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
-
-  //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
-  eDenClus = eClusTot/nClus;
-  eDenCell = eCellTot/nCell;
-  fhEDensityCluster      ->Fill(eDenClus);
-  fhEDensityCell         ->Fill(eDenCell);
-  fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
-  //Fill the average number of cells or clusters per SM
-  eClusTot /=fNModules;
-  eCellTot /=fNModules;
-  fhAverTotECluster      ->Fill(eClusTot);
-  fhAverTotECell         ->Fill(eCellTot);
-  fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
-  //printf("Average Cluster: E %f, density %f;  Average Cell E %f, density  %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
-
-//  //Average weighted pair distance
-//  rxz  /= ncomb;                             
-//  rxzw /= ncomb;                             
-//
-//  fhAverClusterPairDist             ->Fill(rxz );
-//  fhAverClusterPairDistWeight       ->Fill(rxzw);
-//  fhAverClusterPairDistvsAverE      ->Fill(rxz ,eDenClus);
-//  fhAverClusterPairDistWeightvsAverE->Fill(rxzw,eDenClus);
-//  fhAverClusterPairDistvsN          ->Fill(rxz ,nClus);
-//  fhAverClusterPairDistWeightvsN    ->Fill(rxzw,nClus);
-//  
-//  //emax
-//  fhMaxEvsClustEDen->Fill(emax,eDenClus);
-//  fhMaxEvsClustMult->Fill(emax,nPhot);
-  
-  //printf("Average Distance: %f; weighted  %f\n ",rxz,rxzw);
-
   
   //Init variables
   Int_t module1         = -1;
@@ -1723,8 +1753,11 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
 //        //printf("cluster pair rW average %f, bin %d \n",rxzw,curCentrBin);
 //      }     
       else { //Event centrality
-          curCentrBin = GetEventCentrality();
-        //printf("curCentrBin %d\n",curCentrBin);
+          // Centrality task returns at maximum 10, 20 or 100, depending on option chosen and 
+          // number of bins, the bin has to be corrected
+          curCentrBin = GetEventCentrality() * fNCentrBin / GetReader()->GetCentralityOpt(); 
+          if(GetDebug() > 0 )printf("AliAnaPi0::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
+                                    curCentrBin, GetEventCentrality(), fNCentrBin, GetReader()->GetCentralityOpt());
       }
       
       if (curCentrBin < 0 || curCentrBin >= fNCentrBin){ 
@@ -1741,7 +1774,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
       
       //Fill event bin info
       fhEvents    ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
-      fhCentrality->Fill(curCentrBin);
+      if(fNCentrBin > 1) fhCentrality->Fill(curCentrBin);
       currentEvtIndex = evtIndex1 ; 
       if(GetDebug() > 1) 
         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
@@ -1809,7 +1842,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
       //-------------------------------------------------------------------------------------------------
       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
       //-------------------------------------------------------------------------------------------------
-      if(a < fAsymCuts[0]){
+      if(a < fAsymCuts[0] && fFillSMCombinations){
         if(module1==module2 && module1 >=0 && module1<fNModules)
           fhReMod[module1]->Fill(pt,m) ;
 
@@ -1938,7 +1971,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
                     else if(module1==1)  fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
                     else if(module1==2)  fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
                     else if(module1==3)  fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
-                    else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
+                    //else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
                   }
                 }
               }// pid bit cut loop
@@ -2014,7 +2047,7 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
           //-------------------------------------------------------------------------------------------------
           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
           //-------------------------------------------------------------------------------------------------          
-          if(a < fAsymCuts[0]){
+          if(a < fAsymCuts[0] && fFillSMCombinations){
             if(module1==module2 && module1 >=0 && module1<fNModules)
               fhMiMod[module1]->Fill(pt,m) ;
 
@@ -2145,20 +2178,25 @@ void AliAnaPi0::ReadHistograms(TList* outputList)
   if(!fhMiInvPt1) fhMiInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
   if(!fhMiInvPt2) fhMiInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
   if(!fhMiInvPt3) fhMiInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;   
-  if(!fhReMod)    fhReMod     = new TH2D*[fNModules]   ;       
-  if(!fhReDiffPHOSMod)       fhReDiffPHOSMod        = new TH2D*[fNModules] ;   
-  if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ; 
-  if(!fhReSameSideEMCALMod)  fhReSameSideEMCALMod   = new TH2D*[fNModules-2] ; 
-  if(!fhMiMod)    fhMiMod     = new TH2D*[fNModules]   ;       
-  if(!fhMiDiffPHOSMod)       fhMiDiffPHOSMod        = new TH2D*[fNModules] ;   
-  if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ; 
-  if(!fhMiSameSideEMCALMod)  fhMiSameSideEMCALMod   = new TH2D*[fNModules-2] ; 
-    
-  fhReConv  = (TH2D*) outputList->At(index++);
-  fhMiConv  = (TH2D*) outputList->At(index++);
-  fhReConv2 = (TH2D*) outputList->At(index++);
-  fhMiConv2 = (TH2D*) outputList->At(index++);
-
+  
+  if(fFillSMCombinations){
+    if(!fhReMod)    fhReMod     = new TH2D*[fNModules]   ;     
+    if(!fhReDiffPHOSMod)       fhReDiffPHOSMod        = new TH2D*[fNModules] ; 
+    if(!fhReSameSectorEMCALMod)fhReSameSectorEMCALMod = new TH2D*[fNModules/2] ;       
+    if(!fhReSameSideEMCALMod)  fhReSameSideEMCALMod   = new TH2D*[fNModules-2] ;       
+    if(!fhMiMod)    fhMiMod     = new TH2D*[fNModules]   ;     
+    if(!fhMiDiffPHOSMod)       fhMiDiffPHOSMod        = new TH2D*[fNModules] ; 
+    if(!fhMiSameSectorEMCALMod)fhMiSameSectorEMCALMod = new TH2D*[fNModules/2] ;       
+    if(!fhMiSameSideEMCALMod)  fhMiSameSideEMCALMod   = new TH2D*[fNModules-2] ;       
+  }
+  
+  if(fCheckConversion){
+    fhReConv  = (TH2D*) outputList->At(index++);
+    fhMiConv  = (TH2D*) outputList->At(index++);
+    fhReConv2 = (TH2D*) outputList->At(index++);
+    fhMiConv2 = (TH2D*) outputList->At(index++);
+  }
+  
   for(Int_t ic=0; ic<fNCentrBin; ic++){
     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
@@ -2212,7 +2250,8 @@ void AliAnaPi0::ReadHistograms(TList* outputList)
   }// multi cut analysis 
   
   fhEvents     = (TH3D *) outputList->At(index++); 
-  fhCentrality = (TH1D *) outputList->At(index++); 
+  if(fNCentrBin)fhCentrality       = (TH1D *) outputList->At(index++); 
+  if(fNCentrBin)fhCentralityNoPair = (TH1D *) outputList->At(index++); 
 
   fhRealOpeningAngle     = (TH2D*)  outputList->At(index++);
   fhRealCosOpeningAngle  = (TH2D*)  outputList->At(index++);
index 76680a85e0d4ba9312fd59b988563cea4286665f..ed6e72e05b8289eba5ab554e57f040a05f61a624 100755 (executable)
@@ -62,6 +62,8 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
 
   virtual Int_t GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)  ;
 
+  void CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) ;
+  
   //Setters for parameters of event buffers
   void SetNCentrBin(Int_t n=5)    {fNCentrBin=n ;} //number of bins in centrality 
 //void SetNZvertBin(Int_t n=5)    {fNZvertBin=n ;} //number of bins for vertex position
@@ -81,8 +83,8 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   void SwitchOnCellEBins()        {fUseAverCellEBins = kTRUE  ; }
   void SwitchOffCellEBins()       {fUseAverCellEBins = kFALSE ; }
 
-  void SwitchOnClusterEDenBins()     {fUseAverClusterEDenBins = kTRUE  ; }
-  void SwitchOffClusterEDenBins()    {fUseAverClusterEDenBins = kFALSE ; }
+  void SwitchOnClusterEDenBins()  {fUseAverClusterEDenBins = kTRUE  ; }
+  void SwitchOffClusterEDenBins() {fUseAverClusterEDenBins = kFALSE ; }
   
 //  void SwitchOnClusterPairRBins()     {fUseAverClusterPairRBins = kTRUE  ; }
 //  void SwitchOffClusterPairRBins()    {fUseAverClusterPairRBins = kFALSE ; }
@@ -110,11 +112,14 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   void SwitchOffOwnMix()             {fDoOwnMix = kFALSE ; }
 
   //------------------------------------------
-  //Do analysis only with clusters in same SM
+  //Do analysis only with clusters in same SM or different combinations of SM
   //------------------------------------------
   void SwitchOnSameSM()              {fSameSM = kTRUE    ; }
   void SwitchOffSameSM()             {fSameSM = kFALSE   ; }
   
+  void SwitchOnSMCombinations()   {fFillSMCombinations = kTRUE  ; }
+  void SwitchOffSMCombinations()  {fFillSMCombinations = kFALSE ; }
+  
   //-------------------------------
   //Histogram filling options off by default
   //-------------------------------
@@ -151,6 +156,9 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   void SwitchOnMultipleCutAnalysisInSimulation()   {fMultiCutAnaSim = kTRUE;}
   void SwitchOffMultipleCutAnalysisInSimulation()  {fMultiCutAnaSim = kFALSE;}
 
+  void     SwitchOnConversionChecker()             { fCheckConversion = kTRUE  ; }
+  void     SwitchOffConversionChecker()            { fCheckConversion = kFALSE ; }  
+  
   
   private:
   Bool_t   fDoOwnMix;            // Do combinatorial background not the one provided by the frame
@@ -181,6 +189,8 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   //Switchs of different analysis options
   Bool_t   fMakeInvPtPlots;      // D plots with inverse pt weight
   Bool_t   fSameSM;              // Select only pairs in same SM;
+  Bool_t   fFillSMCombinations;  // Fill histograms with different cluster pairs in SM combinations
+  Bool_t   fCheckConversion;     // Fill histograms with tagged photons as conversion
   Bool_t   fUseTrackMultBins;    // Use track multiplicity and not centrality bins
   Bool_t   fUsePhotonMultBins;   // Use photon multiplicity and not centrality bins
   Bool_t   fUseAverCellEBins;    // Use cell average energy and not centrality bins
@@ -258,8 +268,9 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   TH2D *  fhRePtAsymEta ;           //! REAL two-photon pt vs asymmetry, close to eta mass
 
   TH3D * fhEvents;                  //! Number of events per centrality, RP, zbin
-  TH1D * fhCentrality;              //! Simple TH1D histogram with finer binning to check centrality
-  
+  TH1D * fhCentrality;              //! Histogram with centrality bins with at least one pare
+  TH1D * fhCentralityNoPair;        //! Histogram with centrality bins with no pair
+
   // Pair opening angle
   TH2D * fhRealOpeningAngle ;       //! Opening angle of pair versus pair energy
   TH2D * fhRealCosOpeningAngle ;    //! Cosinus of opening angle of pair version pair energy
@@ -307,7 +318,7 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   TH2D *  fhMCPi0PtOrigin ;         //! Mass of reoconstructed pi0 pairs  in calorimeter vs mother
   TH2D *  fhMCEtaPtOrigin ;         //! Mass of reoconstructed pi0 pairs  in calorimeter vs mother
 
-  ClassDef(AliAnaPi0,17)
+  ClassDef(AliAnaPi0,18)
 } ;
 
 
index fe67e60f271828166720ca3da982db87d7f7165e..8fbed4849a9c13a1727404359848e4cb4d0d1f05 100644 (file)
@@ -76,6 +76,17 @@ AliAnalysisTaskParticleCorrelation *AddTaskPi0(TString data, TString calorimeter
   //In case of AODs created only for calorimeters, and track information filtered
   //CTS is off when calling this method
   //reader->SwitchOnCaloFilterPatch();
+  //Event selection
+  if     (col=="pp"  ) {
+    reader->SwitchOnEventSelection(); //remove pileup by default
+    reader->SwitchOffV0ANDSelection() ; // and besides v0 AND
+    reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
+  }
+  else if(col=="PbPb") {
+    reader->SwitchOffEventSelection(); //remove pileup by default
+    reader->SwitchOffV0ANDSelection() ; // and besides v0 AND
+    reader->SwitchOffPrimaryVertexSelection(); // and besides primary vertex
+  }
   
   reader->SetZvertexCut(10.);
   
@@ -166,8 +177,8 @@ AliAnalysisTaskParticleCorrelation *AddTaskPi0(TString data, TString calorimeter
   anaphoton->AddToHistogramsName("AnaPhotonCorr_");
   //Set Histograms bins and ranges
   anaphoton->SetHistoPtRangeAndNBins(0, 30, 150) ;
-  if(year==2010)anaphoton->SetHistoPhiRangeAndNBins(78, 122*TMath::DegToRad(), 44*TMath::DegToRad()) ;
-  else          anaphoton->SetHistoPhiRangeAndNBins(78, 182*TMath::DegToRad(), 104*TMath::DegToRad()) ;
+  if(year==2010)anaphoton->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 122*TMath::DegToRad(), 78) ;
+  else          anaphoton->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 182*TMath::DegToRad(), 108) ;
   anaphoton->SetHistoEtaRangeAndNBins(-0.72, 0.72, 144) ;
   if(kPrintSettings) anaphoton->Print("");
   
@@ -179,17 +190,16 @@ AliAnalysisTaskParticleCorrelation *AddTaskPi0(TString data, TString calorimeter
   anapi0->SetDebug(-1);//10 for lots of messages
   anapi0->SetInputAODName(Form("Photons%s",calorimeter.Data()));
   anapi0->SetCalorimeter(calorimeter);
-  
   anapi0->SwitchOffMultipleCutAnalysis(); 
   //anapi0->SetNPtCuts(2);
-  anapi0->SetNAsymCuts(3);
+  //anapi0->SetNAsymCuts(2);
   //anapi0->SetNNCellCuts(2);
-  anapi0->SetNPIDBits(1);
+  //anapi0->SetNPIDBits(1);
   
   //anapi0->SetPtCutsAt(0,0.3); anapi0->SetPtCutsAt(1,0.5);
   //anapi0->SetAsymCutsAt(0,0.1);anapi0->SetAsymCutsAt(1,0.5);
   //anapi0->SetNCellCutsAt(0,1); anapi0->SetNCellCutsAt(1,2);
-  anapi0->SetPIDBitsAt(0,0); //No Cut
+  //anapi0->SetPIDBitsAt(0,0); //No Cut
   //anapi0->SetPIDBitsAt(1,2); //Dispersion Cut
 
   
@@ -205,9 +215,15 @@ AliAnalysisTaskParticleCorrelation *AddTaskPi0(TString data, TString calorimeter
        
   //settings for pp collision mixing
   anapi0->SwitchOnOwnMix(); //Off when mixing done with general mixing frame
-  if     (col=="pp"  ) anapi0->SetNCentrBin(1);
-  else if(col=="PbPb") anapi0->SetNCentrBin(10);
-  anapi0->SetNZvertBin(1);
+  if     (col=="pp"  ) {
+    anapi0->SetNCentrBin(1);
+    anapi0->SwitchOnSMCombinations();
+  }
+  else if(col=="PbPb") {
+    anapi0->SetNCentrBin(10);
+    anapi0->SwitchOffSMCombinations();
+  }
+  anapi0->SetNZvertBin(10);
   anapi0->SetNRPBin(1);
   anapi0->SetNMaxEvMix(50);
   
@@ -220,8 +236,8 @@ AliAnalysisTaskParticleCorrelation *AddTaskPi0(TString data, TString calorimeter
   }
   anapi0->SetHistoPtRangeAndNBins(0, 30, 150) ;    
   anapi0->SetHistoEtaRangeAndNBins(-0.72, 0.72, 144) ;
-  if(year==2010)anapi0->SetHistoPhiRangeAndNBins(78, 122*TMath::DegToRad(), 44*TMath::DegToRad()) ;
-  else          anapi0->SetHistoPhiRangeAndNBins(78, 182*TMath::DegToRad(), 104*TMath::DegToRad()) ;
+  if(year==2010)anaphoton->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 122*TMath::DegToRad(), 78) ;
+  else          anaphoton->SetHistoPhiRangeAndNBins(78*TMath::DegToRad(), 182*TMath::DegToRad(), 108) ;
 
   anapi0->SetHistoMassRangeAndNBins(0., 1., 200) ;
   anapi0->SetHistoAsymmetryRangeAndNBins(0., 1. , 100) ;