- Removing alot of histograms, replaced with trees instead.
authorodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Jun 2012 00:55:32 +0000 (00:55 +0000)
committerodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Jun 2012 00:55:32 +0000 (00:55 +0000)
- Moving more of the cut selection into the AliAnalysisEtSelector
  class
- Fixing some coverity issues.
- Some PHOS specific code removed from non-PHOS specific classes.

The move to using trees is reasonable due to the low statistics
needed for the E_T analysis.

At the moment EMCAL code will not run. Sorry! Needs to be fixed
ASAP.

18 files changed:
PWGLF/totEt/AliAnalysisEt.cxx
PWGLF/totEt/AliAnalysisEt.h
PWGLF/totEt/AliAnalysisEtCommon.h
PWGLF/totEt/AliAnalysisEtMonteCarlo.cxx
PWGLF/totEt/AliAnalysisEtMonteCarlo.h
PWGLF/totEt/AliAnalysisEtMonteCarloEmcal.cxx
PWGLF/totEt/AliAnalysisEtMonteCarloPhos.cxx
PWGLF/totEt/AliAnalysisEtMonteCarloPhos.h
PWGLF/totEt/AliAnalysisEtReconstructed.cxx
PWGLF/totEt/AliAnalysisEtReconstructed.h
PWGLF/totEt/AliAnalysisEtReconstructedEmcal.cxx
PWGLF/totEt/AliAnalysisEtReconstructedPhos.cxx
PWGLF/totEt/AliAnalysisEtReconstructedPhos.h
PWGLF/totEt/AliAnalysisEtSelector.cxx
PWGLF/totEt/AliAnalysisEtSelector.h
PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx
PWGLF/totEt/AliAnalysisEtSelectorPhos.h
PWGLF/totEt/AliAnalysisTaskTotEt.cxx

index 1936d2b..3cf6428 100644 (file)
@@ -56,204 +56,138 @@ Float_t AliAnalysisEt::fgRAxis[48]={-2.,-1.,0.,0.0005,0.001,0.0015,0.002,0.0025,
 
 
 AliAnalysisEt::AliAnalysisEt() : AliAnalysisEtCommon()
+                              ,fEventSummaryTree(0)
+                              ,fAcceptedTree(0)
+                              ,fDepositTree(0)
                               ,fTotEt(0)
-                              ,fTotEtAcc(0)
                               ,fTotNeutralEt(0)
-                              ,fTotNeutralEtAcc(0)
                               ,fTotChargedEt(0)
-                              ,fTotChargedEtAcc(0)
                               ,fMultiplicity(0)
                               ,fChargedMultiplicity(0)
                               ,fNeutralMultiplicity(0)
-                              ,fBaryonEt(0)
-                              ,fAntiBaryonEt(0)
-                              ,fMesonEt(0)
                               ,fProtonEt(0)
-                              ,fPionEt(0)
-                              ,fChargedKaonEt(0)
-                              ,fMuonEt(0)
-                              ,fElectronEt(0)
+                              ,fAntiProtonEt(0)
                               ,fNeutronEt(0)
                               ,fAntiNeutronEt(0)
+                              ,fPi0Et(0)
+                              ,fPiPlusEt(0)
+                              ,fPiMinusEt(0)
+                              ,fKPlusEt(0)
+                              ,fKMinusEt(0)
+                              ,fK0sEt(0)
+                              ,fK0lEt(0)
+                              ,fMuMinusEt(0)
+                              ,fMuPlusEt(0)
+                              ,fEMinusEt(0)
+                              ,fEPlusEt(0)
                               ,fGammaEt(0)
-                              ,fProtonEtAcc(0)
-                              ,fPionEtAcc(0)
-                              ,fChargedKaonEtAcc(0)
-                              ,fMuonEtAcc(0)
-                              ,fElectronEtAcc(0)
+                              ,fProtonRemovedEt(0)
+                              ,fAntiProtonRemovedEt(0)
+                              ,fNeutronRemovedEt(0)
+                              ,fAntiNeutronRemovedEt(0)
+                              ,fPi0RemovedEt(0)
+                              ,fPiPlusRemovedEt(0)
+                              ,fPiMinusRemovedEt(0)
+                              ,fKPlusRemovedEt(0)
+                              ,fKMinusRemovedEt(0)
+                              ,fK0sRemovedEt(0)
+                              ,fK0lRemovedEt(0)
+                              ,fMuMinusRemovedEt(0)
+                              ,fMuPlusRemovedEt(0)
+                              ,fEMinusRemovedEt(0)
+                              ,fEPlusRemovedEt(0)
+                              ,fGammaRemovedEt(0)
+                              ,fProtonMult(0)
+                              ,fAntiProtonMult(0)
+                              ,fNeutronMult(0)
+                              ,fAntiNeutronMult(0)
+                              ,fPi0Mult(0)
+                              ,fPiPlusMult(0)
+                              ,fPiMinusMult(0)
+                              ,fKPlusMult(0)
+                              ,fKMinusMult(0)
+                              ,fK0sMult(0)
+                              ,fK0lMult(0)
+                              ,fMuMinusMult(0)
+                              ,fMuPlusMult(0)
+                              ,fEMinusMult(0)
+                              ,fEPlusMult(0)
+                              ,fGammaMult(0)
+                              ,fProtonRemovedMult(0)
+                              ,fAntiProtonRemovedMult(0)
+                              ,fNeutronRemovedMult(0)
+                              ,fAntiNeutronRemovedMult(0)
+                              ,fPi0RemovedMult(0)
+                              ,fPiPlusRemovedMult(0)
+                              ,fPiMinusRemovedMult(0)
+                              ,fKPlusRemovedMult(0)
+                              ,fKMinusRemovedMult(0)
+                              ,fK0sRemovedMult(0)
+                              ,fK0lRemovedMult(0)
+                              ,fMuMinusRemovedMult(0)
+                              ,fMuPlusRemovedMult(0)
+                              ,fEMinusRemovedMult(0)
+                              ,fEPlusRemovedMult(0)
+                              ,fGammaRemovedMult(0)
                               ,fEnergyDeposited(0)
-                              ,fEnergyTPC(0)
+                              ,fMomentumTPC(0)
                               ,fCharge(0)
                               ,fParticlePid(0)
                               ,fPidProb(0)
                               ,fTrackPassedCut(kFALSE)
                               ,fCentClass(0)
-                              ,fEtaCut(0)
-                              ,fEtaCutAcc(0)
-                              ,fPhiCutAccMin(0)
-                              ,fPhiCutAccMax(0)
                               ,fDetectorRadius(0)
-                              ,fClusterEnergyCut(0) 
                               ,fSingleCellEnergyCut(0)
-                              ,fTrackDistanceCut(0)
-                              ,fTrackDxCut(0)
-                              ,fTrackDzCut(0)
                               ,fChargedEnergyRemoved(0)
                               ,fNeutralEnergyRemoved(0)
                               ,fGammaEnergyAdded(0)
                               ,fHistEt(0)
-                              ,fHistChargedEt(0)
-                              ,fHistNeutralEt(0)
-                              ,fHistEtAcc(0)
-                              ,fHistChargedEtAcc(0)
-                              ,fHistNeutralEtAcc(0)
-                              ,fHistMult(0)
-                              ,fHistChargedMult(0)
                               ,fHistNeutralMult(0)
                               ,fHistPhivsPtPos(0)
                               ,fHistPhivsPtNeg(0)
-                              ,fHistBaryonEt(0)
-                              ,fHistAntiBaryonEt(0)
-                              ,fHistMesonEt(0)
-                              ,fHistProtonEt(0)
-                              ,fHistPionEt(0)
-                              ,fHistChargedKaonEt(0)
-                              ,fHistMuonEt(0)
-                              ,fHistElectronEt(0)
-                              ,fHistNeutronEt(0)
-                              ,fHistAntiNeutronEt(0)
-                              ,fHistGammaEt(0)
-                              ,fHistProtonEtAcc(0)
-                              ,fHistPionEtAcc(0)
-                              ,fHistChargedKaonEtAcc(0)
-                              ,fHistMuonEtAcc(0)
-                              ,fHistElectronEtAcc(0)
-                              ,fHistTMDeltaR(0)
-                              ,fHistTMDxDz(0)
-                              ,fTree(0)
-                              ,fTreeDeposit(0)
                               ,fCentrality(0)
-                              ,fDetector(0)
-                              ,fMakeSparse(kTRUE)
-                              ,fSparseHistTracks(0)
-                              ,fSparseHistClusters(0)
-                              ,fSparseHistEt(0)
-                                      ,fSparseTracks(0)
-                              ,fSparseClusters(0)
-                              ,fSparseEt(0)
-                              ,fClusterType(0)
-                              ,fMatrixInitialized(kFALSE)
+                              ,fMakeSparse(kFALSE)
                               ,fCutFlow(0)
                               ,fSelector(0)
-                              
+              
 {}
 
 AliAnalysisEt::~AliAnalysisEt()
 {//Destructor
-  if(fTreeDeposit){
-    fTreeDeposit->Clear();
-    delete fTreeDeposit; // optional TTree
+  if(fDepositTree){
+    fDepositTree->Clear();
+    delete fDepositTree; // optional TTree
   }
-  if(fTree){
-    fTree->Clear();
-    delete fTree; // optional TTree
+  if(fEventSummaryTree)
+  {
+    fEventSummaryTree->Clear();
+    delete fEventSummaryTree;
   }
   delete fHistEt; //Et spectrum
-  delete fHistChargedEt; //Charged Et spectrum 
-  delete fHistNeutralEt; //Neutral Et spectrum
-  delete fHistEtAcc; //Et in acceptance
-  delete fHistChargedEtAcc; //Charged Et in acceptance
-  delete fHistNeutralEtAcc; //Et in acceptance
-  delete fHistMult; //Multiplicity
-  delete fHistChargedMult; //Charged multiplicity
   delete fHistNeutralMult; //Neutral multiplicity
   delete fHistPhivsPtPos; //phi vs pT plot for positive tracks
   delete fHistPhivsPtNeg; //phi vs pT plot for negative tracks
-  delete fHistBaryonEt; /** Et of identified baryons */    
-  delete fHistAntiBaryonEt; /** Et of identified anti-baryons */
-  delete fHistMesonEt; /** Et of identified mesons */
-  delete fHistProtonEt; /** Et of identified protons */
-  delete fHistPionEt; /** Et of identified protons */
-  delete fHistChargedKaonEt; /** Et of identified charged kaons */
-  delete fHistMuonEt; /** Et of identified muons */
-  delete fHistElectronEt; /** Et of identified electrons */
-  delete fHistNeutronEt; /** Et of neutrons (MC only for now) */
-  delete fHistAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
-  delete fHistGammaEt; /** Et of gammas (MC only for now) */
-  delete fHistProtonEtAcc; /** Et of identified protons in calorimeter acceptance */    
-  delete fHistPionEtAcc; /** Et of identified protons in calorimeter acceptance */    
-  delete fHistChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */    
-  delete fHistMuonEtAcc; /** Et of identified muons in calorimeter acceptance */
-  delete fHistElectronEtAcc; /** Et of identified electrons in calorimeter acceptance */
-  delete fHistTMDeltaR; /* Track matching plots; Rec only for now */
-  delete fHistTMDxDz; /* Track matching plots; Rec only for now */
-  //arrays for axes were not dynamically created so don't need to be deleted
-  delete fTree;
-  delete fTreeDeposit;
   //delete fCentrality;//this code does not actually own AliCentrality so we don't have to worry about deleting it...  we just borrow it...
-  delete fSparseHistTracks;
-  delete fSparseHistClusters;
-  delete fSparseHistEt;
-  delete [] fSparseTracks;
-  delete [] fSparseClusters;
-  delete [] fSparseEt;
 }
 
 void AliAnalysisEt::FillOutputList(TList *list)
 { // histograms to be added to output
     list->Add(fHistEt);
-    list->Add(fHistChargedEt);
-    list->Add(fHistNeutralEt);
-
-    list->Add(fHistEtAcc);
-    list->Add(fHistChargedEtAcc);
-    list->Add(fHistNeutralEtAcc);
-
-    list->Add(fHistMult);
-    list->Add(fHistChargedMult);
     list->Add(fHistNeutralMult);
 
     list->Add(fHistPhivsPtPos);
     list->Add(fHistPhivsPtNeg);
 
-    list->Add(fHistBaryonEt);
-    list->Add(fHistAntiBaryonEt);
-    list->Add(fHistMesonEt);
-
-    list->Add(fHistProtonEt);
-    list->Add(fHistPionEt);
-    list->Add(fHistChargedKaonEt);
-    list->Add(fHistMuonEt);
-    list->Add(fHistElectronEt);
-    
-    list->Add(fHistNeutronEt);
-    list->Add(fHistAntiNeutronEt);
-    list->Add(fHistGammaEt);
-    
-    list->Add(fHistProtonEtAcc);
-    list->Add(fHistPionEtAcc);
-    list->Add(fHistChargedKaonEtAcc);
-    list->Add(fHistMuonEtAcc);
-    list->Add(fHistElectronEtAcc);
-
-    list->Add(fHistTMDeltaR);
-    list->Add(fHistTMDxDz);
-       
     if (fCuts) {
       if (fCuts->GetHistMakeTree()) {
-       list->Add(fTree);
+       //list->Add(fTree);
+       list->Add(fEventSummaryTree);
       }
       if (fCuts->GetHistMakeTreeDeposit()) {
-       list->Add(fTreeDeposit);
+       list->Add(fDepositTree);
       }
     }
     
-    if(fMakeSparse){
-      list->Add(fSparseHistTracks);
-      list->Add(fSparseHistClusters);
-      list->Add(fSparseHistEt);
-    }
-    
     list->Add(fCutFlow);
 
 }
@@ -301,42 +235,6 @@ void AliAnalysisEt::CreateHistograms()
     fHistEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
     fHistEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
 
-    histname = "fHistChargedEt" + fHistogramNameSuffix;
-    fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", nbinsEt, minEt, maxEt);
-    fHistChargedEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
-    fHistChargedEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
-
-    histname = "fHistNeutralEt" + fHistogramNameSuffix;
-    fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", nbinsEt, minEt, maxEt);
-    fHistNeutralEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
-    fHistNeutralEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
-
-    histname = "fHistEtAcc" + fHistogramNameSuffix;
-    fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
-    fHistEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
-    fHistEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
-
-    histname = "fHistChargedEtAcc" + fHistogramNameSuffix;
-    fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
-    fHistChargedEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
-    fHistChargedEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
-
-    histname = "fHistNeutralEtAcc" + fHistogramNameSuffix;
-    fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
-    fHistNeutralEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
-    fHistNeutralEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
-    std::cout << histname << std::endl;
-
-       histname = "fHistMult" + fHistogramNameSuffix;
-    fHistMult = new TH1F(histname.Data(), "Total Multiplicity", nbinsMult, minMult, maxMult);
-    fHistMult->GetXaxis()->SetTitle("N");
-    fHistMult->GetYaxis()->SetTitle("Multiplicity");
-
-    histname = "fHistChargedMult" + fHistogramNameSuffix;
-    fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", nbinsMult, minMult, maxMult);
-    fHistChargedMult->GetXaxis()->SetTitle("N");
-    fHistChargedMult->GetYaxis()->SetTitle("Multiplicity");
-
     histname = "fHistNeutralMult" + fHistogramNameSuffix;
     fHistNeutralMult = new TH1F(histname.Data(), "Neutral Multiplicity", nbinsMult, minMult, maxMult);
     fHistNeutralMult->GetXaxis()->SetTitle("N");
@@ -348,116 +246,6 @@ void AliAnalysisEt::CreateHistograms()
     histname = "fHistPhivsPtNeg" + fHistogramNameSuffix;
     fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter",      200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
 
-    histname = "fHistBaryonEt" + fHistogramNameSuffix;
-    fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons",  nbinsEt, minEt, maxEt);
-
-    histname = "fHistAntiBaryonEt" + fHistogramNameSuffix;
-    fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons",  nbinsEt, minEt, maxEt);
-
-    histname = "fHistMesonEt" + fHistogramNameSuffix;
-    fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons",  nbinsEt, minEt, maxEt);
-
-    histname = "fHistProtonEt" + fHistogramNameSuffix;
-    fHistProtonEt = new TH1F(histname.Data(), "E_{T} for (anti-)protons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistPionEt" + fHistogramNameSuffix;
-    fHistPionEt = new TH1F(histname.Data(), "E_{T} for #pi^+/#pi^-", nbinsEt, minEt, maxEt);
-
-    histname = "fHistKaonEt" + fHistogramNameSuffix;
-    fHistChargedKaonEt = new TH1F(histname.Data(), "E_{T} for charged kaons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistMuonEt" + fHistogramNameSuffix;
-    fHistMuonEt = new TH1F(histname.Data(), "E_{T} for muons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistElectronEt" + fHistogramNameSuffix;
-    fHistElectronEt = new TH1F(histname.Data(), "E_{T} for electrons/positrons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistNeutronEt" + fHistogramNameSuffix;
-    fHistNeutronEt = new TH1F(histname.Data(), "E_{T} for neutrons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistAntiNeutronEt" + fHistogramNameSuffix;
-    fHistAntiNeutronEt = new TH1F(histname.Data(), "E_{T} for anti-neutrons", nbinsEt, minEt, maxEt);
-
-    histname = "fHistGammaEt" + fHistogramNameSuffix;
-    fHistGammaEt = new TH1F(histname.Data(), "E_{T} for gammas", nbinsEt, minEt, maxEt);
-
-    histname = "fHistProtonEtAcc" + fHistogramNameSuffix;
-    fHistProtonEtAcc = new TH1F(histname.Data(), "E_{T} for (anti-)protons in calorimeter acceptance", nbinsEt, minEt, maxEt);
-
-    histname = "fHistPionEtAcc" + fHistogramNameSuffix;
-    fHistPionEtAcc = new TH1F(histname.Data(), "E_{T} for #pi^+/#pi^- in calorimeter acceptance", nbinsEt, minEt, maxEt);
-
-    histname = "fHistKaonEtAcc" + fHistogramNameSuffix;
-    fHistChargedKaonEtAcc = new TH1F(histname.Data(), "E_{T} for charged kaons in calorimeter acceptance", nbinsEt, minEt, maxEt);
-
-    histname = "fHistMuonEtAcc" + fHistogramNameSuffix;
-    fHistMuonEtAcc = new TH1F(histname.Data(), "E_{T} for muons in calorimeter acceptance", nbinsEt, minEt, maxEt);
-
-    histname = "fHistElectronEtAcc" + fHistogramNameSuffix;
-    fHistElectronEtAcc = new TH1F(histname.Data(), "E_{T} for electrons/positrons in calorimeter acceptance", nbinsEt, minEt, maxEt);
-
-    //
-    histname = "fHistTMDeltaR" + fHistogramNameSuffix;
-    fHistTMDeltaR = new TH1F(histname.Data(), "#Delta R for calorimeter clusters", 200, 0, 50);
-    
-    histname = "fHistTMDxDz" + fHistogramNameSuffix;
-    fHistTMDxDz = new TH2F(histname.Data(), "#Delta x vs #Delta z for calorimeter clusters", 800, -200, 200, 800, -200, 200);
-    
-    if(fMakeSparse){
-      histname = "fSparseHistTracks" + fHistogramNameSuffix;
-      const Int_t stsize = 7;
-      Int_t binsHist[stsize]   = {  1001,    7, 200000, 10000, 10000, 100,   11};
-      Double_t minHist[stsize] = {-500.5, -3.5,    0.0,   0.0,   0.0, -1.5, -0.5};
-      Double_t maxHist[stsize] = { 499.5,  3.5,  200.0, 100.0, 100.0,  1.5, 10.5};
-      fSparseHistTracks = new THnSparseF(histname.Data(), "pid:charge:mass:et:pt:rap:cent", stsize, binsHist, minHist, maxHist);
-      fSparseTracks = new Double_t[stsize];
-    
-      fSparseHistTracks->GetAxis(0)->SetTitle("pid");
-      fSparseHistTracks->GetAxis(1)->SetTitle("charge");
-      fSparseHistTracks->GetAxis(2)->SetTitle("mass");
-      fSparseHistTracks->GetAxis(3)->SetTitle("et");
-      fSparseHistTracks->GetAxis(4)->SetTitle("pt");
-      fSparseHistTracks->GetAxis(5)->SetTitle("rap");
-      fSparseHistTracks->GetAxis(6)->SetTitle("cent");
-
-      histname = "fSparseHistClusters" + fHistogramNameSuffix;
-      const Int_t scsize = 11;
-      //                            pid     ch    mass     et     pt   eta   et_t   pt_t  eta_t  cent   dist
-      Int_t scbinsHist[scsize]   = {  1001,    7, 200000, 10000, 10000,  100, 10000, 10000,  100,   11,   4000};
-      Double_t scminHist[scsize] = {-500.5, -3.5,    0.0,   0.0,   0.0, -1.5,   0.0,   0.0, -1.5, -0.5, -200.0};
-      Double_t scmaxHist[scsize] = { 499.5,  3.5,  200.0, 100.0, 100.0,  1.5, 100.0, 100.0,  1.5, 10.5,  200.0};
-      fSparseClusters = new Double_t[scsize];
-      fSparseHistClusters = new THnSparseF(histname.Data(), "pid:charge:mass:et:pt:rap:et_track:pt_track:eta_track:cent:dist_matched", scsize, scbinsHist, scminHist, scmaxHist);
-    
-      fSparseHistClusters->GetAxis(0)->SetTitle("pid");
-      fSparseHistClusters->GetAxis(1)->SetTitle("charge");
-      fSparseHistClusters->GetAxis(2)->SetTitle("mass");
-      fSparseHistClusters->GetAxis(3)->SetTitle("et");
-      fSparseHistClusters->GetAxis(4)->SetTitle("pt");
-      fSparseHistClusters->GetAxis(5)->SetTitle("rap");
-      fSparseHistClusters->GetAxis(6)->SetTitle("et_track");
-      fSparseHistClusters->GetAxis(7)->SetTitle("pt_track");
-      fSparseHistClusters->GetAxis(8)->SetTitle("rap_track");
-      fSparseHistClusters->GetAxis(9)->SetTitle("cent");
-      fSparseHistClusters->GetAxis(10)->SetTitle("dist_matched");
-
-      histname = "fSparseHistEt" + fHistogramNameSuffix;
-      const Int_t etsize = 7;
-      Int_t etbinsHist[etsize]   = { 10000, 10000,  10000,   3000,   500,   30000,   11};
-      Double_t etminHist[etsize] = {   0.0,   0.0,    0.0,   -0.5,  -0.5,    -0.5, -0.5};
-      Double_t etmaxHist[etsize] = { 200.0, 200.0,  200.0, 2999.5, 499.5,  2999.5, 10.5};
-      fSparseEt = new Double_t[etsize];
-      fSparseHistEt = new THnSparseF(histname.Data(), "tot_et:neutral_et:charged_et:tot_mult:neutral_mult:charged_mult:cent", etsize, etbinsHist, etminHist, etmaxHist);
-    
-      fSparseHistEt->GetAxis(0)->SetTitle("tot_et");
-      fSparseHistEt->GetAxis(1)->SetTitle("neutral_et");
-      fSparseHistEt->GetAxis(2)->SetTitle("charged_et");
-      fSparseHistEt->GetAxis(3)->SetTitle("tot_mult");
-      fSparseHistEt->GetAxis(4)->SetTitle("netral_mult");
-      fSparseHistEt->GetAxis(5)->SetTitle("charged_mult");
-      fSparseHistEt->GetAxis(6)->SetTitle("cent");
-    }
-    
     histname = "fCutFlow" + fHistogramNameSuffix;
     fCutFlow = new TH1I(histname, histname, 20, -0.5, 19.5);
 }
@@ -591,50 +379,130 @@ THnSparseF* AliAnalysisEt::CreateChargedPartHistoSparse(TString name, TString ti
 
 void AliAnalysisEt::CreateTrees()
 { // create tree..
-  TString treename = "fTree" + fHistogramNameSuffix;
+  TString treename = "fEventSummaryTree" + fHistogramNameSuffix;
   if(fCuts->GetHistMakeTree())
   {
   
-    fTree = new TTree(treename, treename);
-    fTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
-    fTree->Branch("fTotEtAcc",&fTotEtAcc,"fTotEtAcc/D");
-    fTree->Branch("fTotNeutralEt",&fTotNeutralEt,"fTotNeutralEt/D");
-    fTree->Branch("fTotNeutralEtAcc",&fTotNeutralEtAcc,"fTotNeutralEtAcc/D");
-    fTree->Branch("fTotChargedEt",&fTotChargedEt,"fTotChargedEt/D");
-    fTree->Branch("fTotChargedEtAcc",&fTotChargedEtAcc,"fTotChargedEtAcc/D");
-    fTree->Branch("fMultiplicity",&fMultiplicity,"fMultiplicity/I");
-    fTree->Branch("fChargedMultiplicity",&fChargedMultiplicity,"fChargedMultiplicity/I");
-    fTree->Branch("fNeutralMultiplicity",&fNeutralMultiplicity,"fNeutralMultiplicity/I");
-    fTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
-    fTree->Branch("fChargedEnergyRemoved", &fChargedEnergyRemoved, "fChargedEnergyRemoved/D");
-    fTree->Branch("fNeutralEnergyRemoved", &fNeutralEnergyRemoved, "fNeutralEnergyRemoved/D");
-    fTree->Branch("fGammaEnergyAdded", &fGammaEnergyAdded, "fGammaEnergyAdded/D");
+    fEventSummaryTree = new TTree(treename, treename);
+    fEventSummaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
+    fEventSummaryTree->Branch("fTotNeutralEt",&fTotNeutralEt,"fTotNeutralEt/D");
+    fEventSummaryTree->Branch("fTotChargedEt",&fTotChargedEt,"fTotChargedEt/D");
+    fEventSummaryTree->Branch("fMultiplicity",&fMultiplicity,"fMultiplicity/I");
+    fEventSummaryTree->Branch("fChargedMultiplicity",&fChargedMultiplicity,"fChargedMultiplicity/I");
+    fEventSummaryTree->Branch("fNeutralMultiplicity",&fNeutralMultiplicity,"fNeutralMultiplicity/I");
+    fEventSummaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
+    fEventSummaryTree->Branch("fChargedEnergyRemoved", &fChargedEnergyRemoved, "fChargedEnergyRemoved/D");
+    fEventSummaryTree->Branch("fNeutralEnergyRemoved", &fNeutralEnergyRemoved, "fNeutralEnergyRemoved/D");
+    fEventSummaryTree->Branch("fGammaEnergyAdded", &fGammaEnergyAdded, "fGammaEnergyAdded/D");
+    
+
+    fEventSummaryTree->Branch("fProtonEt",&fProtonEt,"fProtonEt/D");
+    fEventSummaryTree->Branch("fAntiProtonEt",&fAntiProtonEt,"fAntiProtonEt/D");
+
+    fEventSummaryTree->Branch("fNeutronEt",&fNeutronEt,"fNeutronEt/D");
+    fEventSummaryTree->Branch("fAntiNeutronEt",&fAntiNeutronEt,"fAntiNeutronEt/D");
+
+    fEventSummaryTree->Branch("fPi0Et",&fPi0Et,"fPi0Et/D");
+    fEventSummaryTree->Branch("fPiPlusEt",&fPiPlusEt,"fPiPlusEt/D");
+    fEventSummaryTree->Branch("fPiMinusEt",&fPiMinusEt,"fPiMinusEt/D");
+
+    fEventSummaryTree->Branch("fKPlusEt",&fKPlusEt,"fKPlusEt/D");
+    fEventSummaryTree->Branch("fKMinusEt",&fKMinusEt,"fKMinusEt/D");
+    fEventSummaryTree->Branch("fK0sEt",&fK0sEt,"fK0sEt/D");
+    fEventSummaryTree->Branch("fK0lEt",&fK0lEt,"fK0lEt/D");
+
+    fEventSummaryTree->Branch("fMuMinusEt",&fMuMinusEt,"fMuMinusEt/D");
+    fEventSummaryTree->Branch("fMuPlusEt",&fMuPlusEt,"fMuPlusEt/D");
+    
+    fEventSummaryTree->Branch("fEMinusEt",&fEMinusEt,"fEMinusEt/D");
+    fEventSummaryTree->Branch("fEPlusEt",&fEPlusEt,"fEPlusEt/D");
+    
+    fEventSummaryTree->Branch("fGammaEt", &fGammaEt, "fGammaEt/D");
+    
+    fEventSummaryTree->Branch("fProtonRemovedEt",&fProtonRemovedEt,"fProtonRemovedEt/D");
+    fEventSummaryTree->Branch("fAntiProtonRemovedEt",&fAntiProtonRemovedEt,"fAntiProtonRemovedEt/D");
+
+    fEventSummaryTree->Branch("fNeutronRemovedEt",&fNeutronRemovedEt,"fNeutronRemovedEt/D");
+    fEventSummaryTree->Branch("fAntiNeutronRemovedEt",&fAntiNeutronRemovedEt,"fAntiNeutronRemovedEt/D");
+
+    fEventSummaryTree->Branch("fPi0RemovedEt",&fPi0RemovedEt,"fPi0RemovedEt/D");
+    fEventSummaryTree->Branch("fPiPlusRemovedEt",&fPiPlusRemovedEt,"fPiPlusRemovedEt/D");
+    fEventSummaryTree->Branch("fPiMinusRemovedEt",&fPiMinusRemovedEt,"fPiMinusRemovedEt/D");
+
+    fEventSummaryTree->Branch("fKPlusRemovedEt",&fKPlusRemovedEt,"fKPlusRemovedEt/D");
+    fEventSummaryTree->Branch("fKMinusRemovedEt",&fKMinusEt,"fKMinusRemovedEt/D");
+    fEventSummaryTree->Branch("fK0sRemovedEt",&fK0sEt,"fK0sRemovedEt/D");
+    fEventSummaryTree->Branch("fK0lRemovedEt",&fK0lRemovedEt,"fK0lRemovedEt/D");
+
+    fEventSummaryTree->Branch("fMuMinusRemovedEt",&fMuMinusRemovedEt,"fMuMinusRemovedEt/D");
+    fEventSummaryTree->Branch("fMuPlusRemovedEt",&fMuPlusRemovedEt,"fMuPlusRemovedEt/D");
+    
+    fEventSummaryTree->Branch("fEMinusRemovedEt",&fEMinusRemovedEt,"fEMinusRemovedEt/D");
+    fEventSummaryTree->Branch("fEPlusRemovedEt",&fEPlusRemovedEt,"fEPlusRemovedEt/D");
+    
+    fEventSummaryTree->Branch("fGammaRemovedEt", &fGammaRemovedEt, "fGammaEtRemoved/D");
+
+    fEventSummaryTree->Branch("fProtonMult",&fProtonMult,"fProtonMult/D");
+    fEventSummaryTree->Branch("fAntiProtonMult",&fAntiProtonMult,"fAntiProtonMult/D");
+
+    fEventSummaryTree->Branch("fNeutronMult",&fNeutronMult,"fNeutronMult/D");
+    fEventSummaryTree->Branch("fAntiNeutronMult",&fAntiNeutronMult,"fAntiNeutronMult/D");
+
+    fEventSummaryTree->Branch("fPi0Mult",&fPi0Mult,"fPi0Mult/D");
+    fEventSummaryTree->Branch("fPiPlusMult",&fPiPlusMult,"fPiPlusMult/D");
+    fEventSummaryTree->Branch("fPiMinusMult",&fPiMinusMult,"fPiMinusMult/D");
+
+    fEventSummaryTree->Branch("fKPlusMult",&fKPlusMult,"fKPlusMult/D");
+    fEventSummaryTree->Branch("fKMinusMult",&fKMinusMult,"fKMinusMult/D");
+    fEventSummaryTree->Branch("fK0sMult",&fK0sMult,"fK0sMult/D");
+    fEventSummaryTree->Branch("fK0lMult",&fK0lMult,"fK0lMult/D");
+
+    fEventSummaryTree->Branch("fMuMinusMult",&fMuMinusMult,"fMuMinusMult/D");
+    fEventSummaryTree->Branch("fMuPlusMult",&fMuPlusMult,"fMuPlusMult/D");
+    
+    fEventSummaryTree->Branch("fEMinusMult",&fEMinusMult,"fEMinusMult/D");
+    fEventSummaryTree->Branch("fEPlusMult",&fEPlusMult,"fEPlusMult/D");
+    
+    fEventSummaryTree->Branch("fGammaMult", &fGammaMult, "fGammaMult/D");
+    
+    fEventSummaryTree->Branch("fProtonRemovedMult",&fProtonRemovedMult,"fProtonRemovedMult/D");
+    fEventSummaryTree->Branch("fAntiProtonRemovedMult",&fAntiProtonRemovedMult,"fAntiProtonRemovedMult/D");
+
+    fEventSummaryTree->Branch("fNeutronRemovedMult",&fNeutronRemovedMult,"fNeutronRemovedMult/D");
+    fEventSummaryTree->Branch("fAntiNeutronRemovedMult",&fAntiNeutronRemovedMult,"fAntiNeutronRemovedMult/D");
+
+    fEventSummaryTree->Branch("fPi0RemovedMult",&fPi0RemovedMult,"fPi0RemovedMult/D");
+    fEventSummaryTree->Branch("fPiPlusRemovedMult",&fPiPlusRemovedMult,"fPiPlusRemovedMult/D");
+    fEventSummaryTree->Branch("fPiMinusRemovedMult",&fPiMinusRemovedMult,"fPiMinusRemovedMult/D");
+
+    fEventSummaryTree->Branch("fKPlusRemovedMult",&fKPlusRemovedMult,"fKPlusRemovedMult/D");
+    fEventSummaryTree->Branch("fKMinusRemovedMult",&fKMinusMult,"fKMinusRemovedMult/D");
+    fEventSummaryTree->Branch("fK0sRemovedMult",&fK0sMult,"fK0sRemovedMult/D");
+    fEventSummaryTree->Branch("fK0lRemovedMult",&fK0lRemovedMult,"fK0lRemovedMult/D");
+
+    fEventSummaryTree->Branch("fMuMinusRemovedMult",&fMuMinusRemovedMult,"fMuMinusRemovedMult/D");
+    fEventSummaryTree->Branch("fMuPlusRemovedMult",&fMuPlusRemovedMult,"fMuPlusRemovedMult/D");
+    
+    fEventSummaryTree->Branch("fEMinusRemovedMult",&fEMinusRemovedMult,"fEMinusRemovedMult/D");
+    fEventSummaryTree->Branch("fEPlusRemovedMult",&fEPlusRemovedMult,"fEPlusRemovedMult/D");
+    
+    fEventSummaryTree->Branch("fGammaRemovedMult", &fGammaRemovedMult, "fGammaMultRemoved/D");
+
+    
     
-    fTree->Branch("fBaryonEt",&fBaryonEt,"fBaryonEt/D");
-    fTree->Branch("fAntiBaryonEt",&fAntiBaryonEt,"fAntiBaryonEt/D");
-    fTree->Branch("fMesonEt",&fMesonEt,"fMesonEt/D");
-    fTree->Branch("fProtonEt",&fProtonEt,"fProtonEt/D");
-    fTree->Branch("fChargedKaonEt",&fChargedKaonEt,"fChargedKaonEt/D");
-    fTree->Branch("fMuonEt",&fMuonEt,"fMuonEt/D");
-    fTree->Branch("fElectronEt",&fElectronEt,"fElectronEt/D");
-    fTree->Branch("fProtonEtAcc",&fProtonEtAcc,"fProtonEtAcc/D");
-    fTree->Branch("fChargedKaonEtAcc",&fChargedKaonEtAcc,"fChargedKaonEtAcc/D");
-    fTree->Branch("fMuonEtAcc",&fMuonEtAcc,"fMuonEtAcc/D");
-    fTree->Branch("fElectronEtAcc",&fElectronEtAcc,"fElectronEtAcc/D");
   }
   
   if(fCuts->GetHistMakeTreeDeposit())
   {
     treename = "fTreeDeposit" + fHistogramNameSuffix;
-    fTreeDeposit = new TTree(treename, treename);
-  
-    fTreeDeposit->Branch("fEnergyDeposited", &fEnergyDeposited, "fEnergyDeposited/F");
-    fTreeDeposit->Branch("fEnergyTPC", &fEnergyTPC, "fEnergyTPC/F");
-    fTreeDeposit->Branch("fCharge", &fCharge, "fCharge/S");
-    fTreeDeposit->Branch("fParticlePid", &fParticlePid, "fParticlePid/S");
-    fTreeDeposit->Branch("fPidProb", &fPidProb, "fPidProb/F");
-    fTreeDeposit->Branch("fTrackPassedCut", &fTrackPassedCut, "fTrackPassedCut/B");
+    fDepositTree = new TTree(treename, treename);
+  
+    fDepositTree->Branch("fEnergyDeposited", &fEnergyDeposited, "fEnergyDeposited/F");
+    fDepositTree->Branch("fMomentumTPC", &fMomentumTPC, "fMomentumTPC/F");
+    fDepositTree->Branch("fCharge", &fCharge, "fCharge/S");
+    fDepositTree->Branch("fParticlePid", &fParticlePid, "fParticlePid/S");
+    fDepositTree->Branch("fPidProb", &fPidProb, "fPidProb/F");
+    fDepositTree->Branch("fTrackPassedCut", &fTrackPassedCut, "fTrackPassedCut/B");
   }
 
   return;
@@ -642,50 +510,28 @@ void AliAnalysisEt::CreateTrees()
 void AliAnalysisEt::FillHistograms()
 { // fill histograms..
     fHistEt->Fill(fTotEt);
-    fHistChargedEt->Fill(fTotChargedEt);
-    fHistNeutralEt->Fill(fTotNeutralEt);
-
-    fHistEtAcc->Fill(fTotEtAcc);
-    fHistChargedEtAcc->Fill(fTotChargedEtAcc);
-    fHistNeutralEtAcc->Fill(fTotNeutralEtAcc);
 
-    fHistMult->Fill(fMultiplicity);
-    fHistChargedMult->Fill(fChargedMultiplicity);
     fHistNeutralMult->Fill(fNeutralMultiplicity);
 
-    fHistBaryonEt->Fill(fBaryonEt);
-    fHistAntiBaryonEt->Fill(fAntiBaryonEt);
-    fHistMesonEt->Fill(fMesonEt);
-
-    fHistProtonEt->Fill(fProtonEt);
-    fHistPionEt->Fill(fPionEt);
-    fHistChargedKaonEt->Fill(fChargedKaonEt);
-    fHistMuonEt->Fill(fMuonEt);
-    fHistElectronEt->Fill(fElectronEt);
-    fHistNeutronEt->Fill(fNeutronEt);
-    fHistAntiNeutronEt->Fill(fAntiNeutronEt);
-    fHistGammaEt->Fill(fGammaEt);
-    
-    fHistProtonEtAcc->Fill(fProtonEtAcc);
-    fHistPionEtAcc->Fill(fPionEtAcc);
-    fHistChargedKaonEtAcc->Fill(fChargedKaonEtAcc);
-    fHistMuonEtAcc->Fill(fMuonEtAcc);
-    fHistElectronEtAcc->Fill(fElectronEtAcc);
-
     if (fCuts) {
       if (fCuts->GetHistMakeTree()) {
-       fTree->Fill();
+       fEventSummaryTree->Fill();
       }
     }
-    
-    
-    if(fMakeSparse){fSparseHistEt->Fill(fSparseEt);}
+    if(fCuts)
+    {
+      if(fCuts->GetHistMakeTreeDeposit())
+      {
+       fDepositTree->Fill();
+      }
+    }
+      
 }
 
 Int_t AliAnalysisEt::AnalyseEvent(AliVEvent *event)
 { //this line is basically here to eliminate a compiler warning that event is not used.  Making it a virtual function did not work with the plugin.
-  fSelector->SetEvent(event);
   AliAnalysisEtCommon::AnalyseEvent(event);
+  //fSelector->SetEvent(event);
   ResetEventValues();
   return 0;
 }
@@ -694,30 +540,105 @@ void AliAnalysisEt::ResetEventValues()
 { // clear
   AliAnalysisEtCommon::ResetEventValues();
   fTotEt = 0;
-  fTotEtAcc = 0;
   fTotNeutralEt = 0;
-  fTotNeutralEtAcc = 0;
   fTotChargedEt  = 0;
-  fTotChargedEtAcc = 0;
   fMultiplicity = 0;
   fChargedMultiplicity = 0;
   fNeutralMultiplicity = 0;
-  fBaryonEt = 0;
-  fAntiBaryonEt = 0;
-  fMesonEt = 0;
+  
   fProtonEt = 0;
-  fPionEt = 0;
-  fChargedKaonEt = 0;
-  fMuonEt = 0;
-  fElectronEt = 0;
+  fAntiProtonEt = 0;
+  
   fNeutronEt = 0;
   fAntiNeutronEt = 0;
+  
+  fPi0Et = 0;
+  fPiPlusEt = 0;
+  fPiMinusEt = 0;
+  
+  fKPlusEt = 0;
+  fKMinusEt = 0;
+  fK0sEt = 0;
+  fK0lEt = 0;
+  
+  fMuMinusEt = 0;
+  fMuPlusEt = 0;
+  
+  fEMinusEt = 0;
+  fEPlusEt = 0;
+  
   fGammaEt = 0;
-  fProtonEtAcc = 0;
-  fPionEtAcc = 0;
-  fChargedKaonEtAcc = 0;
-  fMuonEtAcc = 0;
-  fElectronEtAcc = 0;
+  
+  fProtonRemovedEt = 0;
+  fAntiProtonRemovedEt = 0;
+  
+  fNeutronRemovedEt = 0;
+  fAntiNeutronRemovedEt = 0;
+  
+  fPi0RemovedEt = 0;
+  fPiPlusRemovedEt = 0;
+  fPiMinusRemovedEt = 0;
+  
+  fKPlusRemovedEt = 0;
+  fKMinusRemovedEt = 0;
+  fK0sRemovedEt = 0;
+  fK0lRemovedEt = 0;
+  
+  fMuMinusRemovedEt = 0;
+  fMuPlusRemovedEt = 0;
+  
+  fEMinusRemovedEt = 0;
+  fEPlusRemovedEt = 0;
+  
+  fGammaRemovedEt = 0;
+  
+  
+  fProtonMult = 0;
+  fAntiProtonMult = 0;
+  
+  fNeutronMult = 0;
+  fAntiNeutronMult = 0;
+  
+  fPi0Mult = 0;
+  fPiPlusMult = 0;
+  fPiMinusMult = 0;
+  
+  fKPlusMult = 0;
+  fKMinusMult = 0;
+  fK0sMult = 0;
+  fK0lMult = 0;
+  
+  fMuMinusMult = 0;
+  fMuPlusMult = 0;
+  
+  fEMinusMult = 0;
+  fEPlusMult = 0;
+  
+  fGammaMult = 0;
+  
+  fProtonRemovedMult = 0;
+  fAntiProtonRemovedMult = 0;
+  
+  fNeutronRemovedMult = 0;
+  fAntiNeutronRemovedMult = 0;
+  
+  fPi0RemovedMult = 0;
+  fPiPlusRemovedMult = 0;
+  fPiMinusRemovedMult = 0;
+  
+  fKPlusRemovedMult = 0;
+  fKMinusRemovedMult = 0;
+  fK0sRemovedMult = 0;
+  fK0lRemovedMult = 0;
+  
+  fMuMinusRemovedMult = 0;
+  fMuPlusRemovedMult = 0;
+  
+  fEMinusRemovedMult = 0;
+  fEPlusRemovedMult = 0;
+  
+  fGammaRemovedMult = 0;
+  
   return;
 }
 
@@ -742,35 +663,3 @@ Double_t AliAnalysisEt::CalculateTransverseEnergy(AliESDCaloCluster* cluster)
   return corrEnergy * TMath::Sin(cp.Theta());
 }
 
-//____________________________________________________________________________
-Double_t AliAnalysisEt::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, AliVEvent* e) const {
-  //Parameterization of LHC10h period
-  //_true if neutral_
-  
-  Double_t meanX=0;
-  Double_t meanZ=0.;
-  Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
-              6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
-  Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
-  
-  Double_t mf = e->GetMagneticField(); //Positive for ++ and negative for --
-
-  if(mf<0.){ //field --
-    meanZ = -0.468318 ;
-    if(charge>0)
-      meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
-    else
-      meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
-  }
-  else{ //Field ++
-    meanZ= -0.468318;
-    if(charge>0)
-      meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
-    else
-      meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;     
-  }
-
-  Double_t rz=(dz-meanZ)/sz ;
-  Double_t rx=(dx-meanX)/sx ;
-  return TMath::Sqrt(rx*rx+rz*rz) ;
-}
index 5ca27b1..61114f9 100644 (file)
@@ -76,31 +76,16 @@ public:
         return fTotEt;
     }
 
-    /** Total Et in the event within the acceptance cuts */
-    Double_t GetTotEtAcc() const {
-        return fTotEtAcc;
-    }
-
     /** Total neutral Et in the event (without acceptance cuts) */
     Double_t GetTotNeutralEt() const {
         return fTotNeutralEt;
     }
 
-    /** Total neutral Et in the event within the acceptance cuts */
-    Double_t GetTotNeutralEtAcc() const {
-        return fTotNeutralEtAcc;
-    }
-
     /** Total charged Et in the event (without acceptance cuts) */
     Double_t GetTotChargedEt() const {
         return fTotChargedEt;
     }
 
-    /** Total charged Et in the event within the acceptance cuts */
-    Double_t GetTotChargedEtAcc() const {
-        return fTotChargedEtAcc;
-    }
-
     void SetTPCOnlyTrackCuts(const AliESDtrackCuts *cuts) {
         fEsdtrackCutsTPC = (AliESDtrackCuts *) cuts;
     }
@@ -125,6 +110,11 @@ public:
         return 0;
     }
 
+    /** Get contribution from secondaries */
+    virtual Double_t GetSecondaryContribution(Int_t /*clusterMultiplicity*/) {
+        return 0;
+    }
+
     void MakeSparseHistograms() {
         fMakeSparse=kTRUE;
     }
@@ -135,44 +125,117 @@ protected:
 
     //AliAnalysisEtCuts *fCuts; // keeper of basic cuts
     Double_t CalculateTransverseEnergy(AliESDCaloCluster *cluster);
+    
+    TTree *fEventSummaryTree; // Contains event level information
 
-    Double_t TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, AliVEvent *e) const;
+    TTree *fAcceptedTree; // Tree for information about accepted particles
     
-    virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &/*cluster*/) const { return kFALSE; }
+    TTree *fDepositTree; // optional TTree for energy deposit measurements
     
     Double_t fTotEt;/** Total Et in the event (without acceptance cuts) */
-    Double_t fTotEtAcc;/** Total Et in the event within the acceptance cuts */
 
     Double_t fTotNeutralEt;/** Total neutral Et in the event */
-    Double_t fTotNeutralEtAcc;/** Total neutral Et in the event within the acceptance cuts */
+
     Double_t fTotChargedEt;/** Total charged Et in the event */
-    Double_t fTotChargedEtAcc;/** Total charged Et in the event within the acceptance cuts */
 
     Int_t fMultiplicity;/** Multiplicity of particles in the event */
     Int_t fChargedMultiplicity;/** Multiplicity of charged particles in the event */
     Int_t fNeutralMultiplicity; /** Multiplicity of neutral particles in the event */
 
-    Double_t fBaryonEt;     /** Et of identified baryons; calo based (Rec only for now) */
-    Double_t fAntiBaryonEt; /** Et of identified anti-baryons; calo based (Rec only for now) */
-    Double_t fMesonEt;     /** Et of identified mesons; calo based (Rec only for now) */
-
     Double_t fProtonEt; /** Et of identified protons */
-    Double_t fPionEt; /** Et of identified pions */
-    Double_t fChargedKaonEt; /** Et of identified charged kaons */
-    Double_t fMuonEt; /** Et of identified muons */
-    Double_t fElectronEt; /** Et of identified electrons */
+    Double_t fAntiProtonEt; /** Et of identified protons */
+
     Double_t fNeutronEt; /** Et of neutrons (MC only for now) */
     Double_t fAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
+    
+    Double_t fPi0Et; // Et of identified pi0
+    Double_t fPiPlusEt; // Et of identified pi+
+    Double_t fPiMinusEt; // Et of identified pi-
+    
+    Double_t fKPlusEt; // Et of identified K+ 
+    Double_t fKMinusEt; // Et of identified K- 
+    Double_t fK0sEt; // Et of identified K0 short
+    Double_t fK0lEt; // Et of identified K0 long
+    
+    Double_t fMuMinusEt; // Et of identified mu- 
+    Double_t fMuPlusEt; // Et of identified mu+ 
+
+    Double_t fEMinusEt; // Et of identified e-
+    Double_t fEPlusEt; // Et of identified e+
+    
     Double_t fGammaEt; /** Et of identified electrons (MC only for now) */
 
-    Double_t fProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
-    Double_t fPionEtAcc; /** Et of identified pions in calorimeter acceptance */
-    Double_t fChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
-    Double_t fMuonEtAcc; /** Et of identified muons in calorimeter acceptance */
-    Double_t fElectronEtAcc; /** Et of identified electrons in calorimeter acceptance */
+    Double_t fProtonRemovedEt; /** Et of identified protons */
+    Double_t fAntiProtonRemovedEt; /** Et of identified protons */
+
+    Double_t fNeutronRemovedEt; /** Et of neutrons (MC only for now) */
+    Double_t fAntiNeutronRemovedEt; /** Et of anti-neutrons (MC only for now) */
+    
+    Double_t fPi0RemovedEt; // Removed Et of identified pi0
+    Double_t fPiPlusRemovedEt; // Removed Et of identified pi+
+    Double_t fPiMinusRemovedEt; // Removed Et of identified pi-
+
+    Double_t fKPlusRemovedEt; // Removed Et of identified K+ 
+    Double_t fKMinusRemovedEt; // Removed Et of identified K- 
+    Double_t fK0sRemovedEt; // Removed Et of identified K0 short
+    Double_t fK0lRemovedEt; // Removed Et of identified K0 long
+    
+    Double_t fMuMinusRemovedEt; // Removed Et of identified mu- 
+    Double_t fMuPlusRemovedEt; // Removed Et of identified mu+ 
+    
+    Double_t fEMinusRemovedEt; // Removed Et of identified e-
+    Double_t fEPlusRemovedEt; // Removed Et of identified e+
+    
+    Double_t fGammaRemovedEt; /** Removed Et of identified electrons (MC only for now) */
+
+    Double_t fProtonMult; /** Mult of identified protons */
+    Double_t fAntiProtonMult; /** Mult of identified protons */
+
+    Double_t fNeutronMult; /** Mult of neutrons (MC only for now) */
+    Double_t fAntiNeutronMult; /** Mult of anti-neutrons (MC only for now) */
+    
+    Double_t fPi0Mult; // Mult of identified pi0
+    Double_t fPiPlusMult; // Mult of identified pi+
+    Double_t fPiMinusMult; // Mult of identified pi-
+    
+    Double_t fKPlusMult; // Mult of identified K+ 
+    Double_t fKMinusMult; // Mult of identified K- 
+    Double_t fK0sMult; // Mult of identified K0 short
+    Double_t fK0lMult; // Mult of identified K0 long
+    
+    Double_t fMuMinusMult; // Mult of identified mu- 
+    Double_t fMuPlusMult; // Mult of identified mu+ 
+
+    Double_t fEMinusMult; // Mult of identified e-
+    Double_t fEPlusMult; // Mult of identified e+
+    
+    Double_t fGammaMult; /** Mult of identified electrons (MC only for now) */
+
+    Double_t fProtonRemovedMult; /** Mult of identified protons */
+    Double_t fAntiProtonRemovedMult; /** Mult of identified protons */
+
+    Double_t fNeutronRemovedMult; /** Mult of neutrons (MC only for now) */
+    Double_t fAntiNeutronRemovedMult; /** Mult of anti-neutrons (MC only for now) */
+    
+    Double_t fPi0RemovedMult; // Removed Mult of identified pi0
+    Double_t fPiPlusRemovedMult; // Removed Mult of identified pi+
+    Double_t fPiMinusRemovedMult; // Removed Mult of identified pi-
+
+    Double_t fKPlusRemovedMult; // Removed Mult of identified K+ 
+    Double_t fKMinusRemovedMult; // Removed Mult of identified K- 
+    Double_t fK0sRemovedMult; // Removed Mult of identified K0 short
+    Double_t fK0lRemovedMult; // Removed Mult of identified K0 long
+    
+    Double_t fMuMinusRemovedMult; // Removed Mult of identified mu- 
+    Double_t fMuPlusRemovedMult; // Removed Mult of identified mu+ 
+    
+    Double_t fEMinusRemovedMult; // Removed Mult of identified e-
+    Double_t fEPlusRemovedMult; // Removed Mult of identified e+
+    
+    Double_t fGammaRemovedMult; /** Removed Mult of identified electrons (MC only for now) */
 
     Float_t fEnergyDeposited; /** Energy deposited in calorimeter */
-    Float_t fEnergyTPC; /** Energy measured in TPC */
+    Float_t fMomentumTPC; /** Momentum measured in TPC */
     Short_t fCharge; /** Charge of the particle */
     Short_t fParticlePid; /** Particle PID */
     Float_t fPidProb; /** Probability of PID */
@@ -180,94 +243,28 @@ protected:
 
     Int_t fCentClass; // centrality class
 
-    Double_t fEtaCut;/** Cut in eta (standard |eta| < 0.5 )*/
-
-    /** Eta cut for our acceptance */
-    Double_t fEtaCutAcc; // Eta cut for our acceptance
-
-    /** Min phi cut for our acceptance in radians */
-    Double_t fPhiCutAccMin; // Min phi cut for our acceptance in radians
-
-    /** Max phi cut for our acceptance in radians */
-    Double_t fPhiCutAccMax; // Max phi cut for our acceptance in radians
-
     /** Detector radius */
     Double_t fDetectorRadius; // Detector radius
 
-    /** Cut on the cluster energy */
-    Double_t fClusterEnergyCut; // Cut on the cluster energy
-
     /** Minimum energy to cut on single cell cluster */
     Double_t fSingleCellEnergyCut;  // Minimum energy to cut on single cell cluster
 
-    Double_t fTrackDistanceCut; // cut on track distance
-
-    Double_t fTrackDxCut; // cut on track distance in x
-
-    Double_t fTrackDzCut; // cut on track distance in z
-    
     Double_t fChargedEnergyRemoved;
     Double_t fNeutralEnergyRemoved;
     Double_t fGammaEnergyAdded;
-    
 
     // Declare the histograms
 
-    /** The full Et spectrum measured */
+    /** The EM Et spectrum measured */
     TH1F *fHistEt; //Et spectrum
 
-    /** The full charged Et spectrum measured */
-    TH1F *fHistChargedEt; //Charged Et spectrum
-
-    /** The full neutral Et spectrum measured */
-    TH1F *fHistNeutralEt; //Neutral Et spectrum
-
-    /** The Et spectrum within the calorimeter acceptance */
-    TH1F *fHistEtAcc; //Et in acceptance
-
-    /** The charged Et spectrum within the calorimeter acceptance */
-    TH1F *fHistChargedEtAcc; //Charged Et in acceptance
+    /** Multiplicity of neutral particles in the events */
+    TH1F *fHistNeutralMult; //Multiplicity
 
-    /** The neutral Et spectrum within the calorimeter acceptance */
-    TH1F *fHistNeutralEtAcc; //Et in acceptance
-
-    /** Multiplicity of particles in the events */
-    TH1F *fHistMult; //Multiplicity
-
-    /** Charged multiplicity of particles in the events */
-    TH1F *fHistChargedMult; //Charged multiplicity
-
-    /** Neutral multiplicity of particles in the events */
-    TH1F *fHistNeutralMult; //Neutral multiplicity
-
-    /* Acceptance plots */
+    // Acceptance plots 
     TH2F *fHistPhivsPtPos; //phi vs pT plot for positive tracks
     TH2F *fHistPhivsPtNeg; //phi vs pT plot for negative tracks
 
-    /* PID plots */
-    TH1F *fHistBaryonEt; /** Et of identified baryons */
-    TH1F *fHistAntiBaryonEt; /** Et of identified anti-baryons */
-    TH1F *fHistMesonEt; /** Et of identified mesons */
-
-    TH1F *fHistProtonEt; /** Et of identified protons */
-    TH1F *fHistPionEt; /** Et of identified protons */
-    TH1F *fHistChargedKaonEt; /** Et of identified charged kaons */
-    TH1F *fHistMuonEt; /** Et of identified muons */
-    TH1F *fHistElectronEt; /** Et of identified electrons */
-    TH1F *fHistNeutronEt; /** Et of neutrons (MC only for now) */
-    TH1F *fHistAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
-    TH1F *fHistGammaEt; /** Et of gammas (MC only for now) */
-
-    TH1F *fHistProtonEtAcc; /** Et of identified protons in calorimeter acceptance */
-    TH1F *fHistPionEtAcc; /** Et of identified protons in calorimeter acceptance */
-    TH1F *fHistChargedKaonEtAcc; /** Et of identified charged kaons in calorimeter acceptance */
-    TH1F *fHistMuonEtAcc; /** Et of identified muons in calorimeter acceptance */
-    TH1F *fHistElectronEtAcc; /** Et of identified electrons in calorimeter acceptance */
-
-    /* Correction plots */
-    TH1F *fHistTMDeltaR; /* Track matching plots; Rec only for now */
-    TH2F *fHistTMDxDz; /* Track matching plots; Rec only for now */
-
     /* Auxiliary Histogram variables */
     static Float_t fgEtaAxis[17];//bins for eta axis of histograms
     static Int_t fgnumOfEtaBins;//number of eta bins
@@ -278,47 +275,21 @@ protected:
     static Float_t fgRAxis[48];//bins for R axis
     static Int_t fgNumOfRBins;//number of R bins
 
-
-    TTree *fTree; // optional TTree
-    TTree *fTreeDeposit; // optional TTree for energy deposit measurements
-
     /** Centrality object */
     AliCentrality *fCentrality; //Centrality object
 
-    Short_t fDetector;     /** Which detector? (-1 -> PHOS, 1 -> EMCAL)*/
-
     Bool_t fMakeSparse;//Boolean for whether or not to make sparse histograms
 
-    THnSparseF *fSparseHistTracks;     /** THnSparse histograms */
-
-    THnSparseF *fSparseHistClusters;     /** THnSparse histograms */
-
-    /** ET sparse valuses */
-    THnSparseF *fSparseHistEt; //!
-
-    /** Values for sparse hists */
-    Double_t *fSparseTracks; //!
-
-    /** Values for sparse hists */
-    Double_t *fSparseClusters; //!
-
-    /** ET sparse valuses */
-    Double_t *fSparseEt; //!
-
-    Char_t fClusterType; // selection on cluster type
-    
-    Bool_t fMatrixInitialized;
-    
     TH1I *fCutFlow; // Cut flow
     
-    AliAnalysisEtSelector *fSelector;
+    AliAnalysisEtSelector *fSelector; // Selector class
 
 private:
     //Declare private to avoid compilation warning
     AliAnalysisEt & operator = (const AliAnalysisEt & g) ;//cpy assignment
     AliAnalysisEt(const AliAnalysisEt & g) ; // cpy ctor
 
-    ClassDef(AliAnalysisEt, 1);
+    ClassDef(AliAnalysisEt, 3);
 };
 
 #endif // ALIANALYSISET_H
index d816082..76e38d9 100644 (file)
@@ -36,7 +36,6 @@ public:
     /** Initialise the analysis, must be overloaded. */
     virtual void Init();
 
-
     /** Reset event specific values (Et etc.) */
     virtual void ResetEventValues();
 
index 379a825..d4d1a85 100644 (file)
@@ -9,6 +9,8 @@
 
 #include "AliAnalysisEtMonteCarlo.h"
 #include "AliAnalysisEtCuts.h"
+#include "AliAnalysisEtSelector.h"
+#include "AliAnalysisEtSelector.h"
 #include "AliESDtrack.h"
 #include "AliStack.h"
 #include "AliVEvent.h"
@@ -34,267 +36,270 @@ ClassImp(AliAnalysisEtMonteCarlo);
 
 // ctor
 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
-        ,fImpactParameter(0)
-        ,fNcoll(0)
-        ,fNpart(0)
-
-        ,fHistDecayVertexNonRemovedCharged(0)
-        ,fHistDecayVertexRemovedCharged(0)
-        ,fHistDecayVertexNonRemovedNeutral(0)
-        ,fHistDecayVertexRemovedNeutral(0)
-
-        ,fHistRemovedOrNot(0)
-        ,fHistEtNonRemovedProtons(0)
-        ,fHistEtNonRemovedAntiProtons(0)
-        ,fHistEtNonRemovedPiPlus(0)
-        ,fHistEtNonRemovedPiMinus(0)
-        ,fHistEtNonRemovedKaonPlus(0)
-        ,fHistEtNonRemovedKaonMinus(0)
-        ,fHistEtNonRemovedK0s(0)
-       ,fHistEtNonRemovedK0L(0)
-        ,fHistEtNonRemovedLambdas(0)
-        ,fHistEtNonRemovedElectrons(0)
-        ,fHistEtNonRemovedPositrons(0)
-        ,fHistEtNonRemovedMuPlus(0)
-        ,fHistEtNonRemovedMuMinus(0)
-        ,fHistEtNonRemovedNeutrons(0)
-        ,fHistEtNonRemovedAntiNeutrons(0)
-        ,fHistEtNonRemovedGammas(0)
-        ,fHistEtNonRemovedGammasFromPi0(0)
-        
-       ,fHistEtRemovedGammas(0)
-        ,fHistEtRemovedNeutrons(0)
-        ,fHistEtRemovedAntiNeutrons(0)
-       
-       ,fHistEtRemovedCharged(0)
-       ,fHistEtRemovedNeutrals(0)
-       ,fHistEtNonRemovedCharged(0)
-       ,fHistEtNonRemovedNeutrals(0)
-       
-        ,fHistMultNonRemovedProtons(0)
-        ,fHistMultNonRemovedAntiProtons(0)
-        ,fHistMultNonRemovedPiPlus(0)
-        ,fHistMultNonRemovedPiMinus(0)
-        ,fHistMultNonRemovedKaonPlus(0)
-        ,fHistMultNonRemovedKaonMinus(0)
-        ,fHistMultNonRemovedK0s(0)
-       ,fHistMultNonRemovedK0L(0)
-        ,fHistMultNonRemovedLambdas(0)
-        ,fHistMultNonRemovedElectrons(0)
-        ,fHistMultNonRemovedPositrons(0)
-        ,fHistMultNonRemovedMuPlus(0)
-        ,fHistMultNonRemovedMuMinus(0)
-        ,fHistMultNonRemovedNeutrons(0)
-        ,fHistMultNonRemovedAntiNeutrons(0)
-        ,fHistMultNonRemovedGammas(0)
-       
-        ,fHistMultRemovedGammas(0)
-        ,fHistMultRemovedNeutrons(0)
-        ,fHistMultRemovedAntiNeutrons(0)
-       
-       ,fHistMultRemovedCharged(0)
-       ,fHistMultRemovedNeutrals(0)
-       
-       ,fHistMultNonRemovedCharged(0)
-       ,fHistMultNonRemovedNeutrals(0)
-       
-        ,fHistTrackMultvsNonRemovedCharged(0)
-        ,fHistTrackMultvsNonRemovedNeutral(0)
-        ,fHistTrackMultvsRemovedGamma(0)
-        ,fHistClusterMultvsNonRemovedCharged(0)
-        ,fHistClusterMultvsNonRemovedNeutral(0)
-        ,fHistClusterMultvsRemovedGamma(0)
-        ,fHistMultvsNonRemovedChargedE(0)
-        ,fHistMultvsNonRemovedNeutralE(0)
-        ,fHistMultvsRemovedGammaE(0)
-        ,fEtNonRemovedProtons(0)
-        ,fEtNonRemovedAntiProtons(0)
-        ,fEtNonRemovedPiPlus(0)
-        ,fEtNonRemovedPiMinus(0)
-        ,fEtNonRemovedKaonPlus(0)
-        ,fEtNonRemovedKaonMinus(0)
-        ,fEtNonRemovedK0S(0)
-       ,fEtNonRemovedK0L(0)
-        ,fEtNonRemovedLambdas(0)
-        ,fEtNonRemovedElectrons(0)
-        ,fEtNonRemovedPositrons(0)
-        ,fEtNonRemovedMuMinus(0)
-        ,fEtNonRemovedMuPlus(0)
-        ,fEtNonRemovedGammas(0)
-        ,fEtNonRemovedGammasFromPi0(0)
-        ,fEtNonRemovedNeutrons(0)
-        ,fEtNonRemovedAntiNeutrons(0)
-       
-       ,fEtRemovedProtons(0)
-       ,fEtRemovedAntiProtons(0)
-       ,fEtRemovedPiPlus(0)
-       ,fEtRemovedPiMinus(0)
-       ,fEtRemovedKaonPlus(0)
-       ,fEtRemovedKaonMinus(0)
-       ,fEtRemovedK0s(0)
-       ,fEtRemovedK0L(0)
-       ,fEtRemovedLambdas(0)
-       ,fEtRemovedElectrons(0)
-       ,fEtRemovedPositrons(0)
-       ,fEtRemovedMuMinus(0)
-       ,fEtRemovedMuPlus(0)
-       
-       ,fEtRemovedGammasFromPi0(0)
-        ,fEtRemovedGammas(0)
-        ,fEtRemovedNeutrons(0)
-        ,fEtRemovedAntiNeutrons(0)
-        ,fMultNonRemovedProtons(0)
-        ,fMultNonRemovedAntiProtons(0)
-        ,fMultNonRemovedPiPlus(0)
-        ,fMultNonRemovedPiMinus(0)
-        ,fMultNonRemovedKaonPlus(0)
-        ,fMultNonRemovedKaonMinus(0)
-        ,fMultNonRemovedK0s(0)
-       ,fMultNonRemovedK0L(0)
-        ,fMultNonRemovedLambdas(0)
-        ,fMultNonRemovedElectrons(0)
-        ,fMultNonRemovedPositrons(0)
-        ,fMultNonRemovedMuMinus(0)
-        ,fMultNonRemovedMuPlus(0)
-        ,fMultNonRemovedGammas(0)
-        ,fMultNonRemovedNeutrons(0)
-        ,fMultNonRemovedAntiNeutrons(0)
-       
-       ,fMultRemovedProtons(0)
-       ,fMultRemovedAntiProtons(0)
-       ,fMultRemovedPiPlus(0)
-       ,fMultRemovedPiMinus(0)
-       ,fMultRemovedKaonPlus(0)
-       ,fMultRemovedKaonMinus(0)
-       ,fMultRemovedK0s(0)
-       ,fMultRemovedK0L(0)
-       ,fMultRemovedLambdas(0)
-       ,fMultRemovedElectrons(0)
-       ,fMultRemovedPositrons(0)
-       ,fMultRemovedMuMinus(0)
-       ,fMultRemovedMuPlus(0)
-       
-       
-        ,fMultRemovedGammas(0)
-        ,fMultRemovedNeutrons(0)
-        ,fMultRemovedAntiNeutrons(0)
-       
-        ,fTrackMultInAcc(0)
-        ,fHistDxDzNonRemovedCharged(0)
-        ,fHistDxDzRemovedCharged(0)
-        ,fHistDxDzNonRemovedNeutral(0)
-        ,fHistDxDzRemovedNeutral(0)
-       ,fHistPiPlusMult(0)
-       ,fHistPiMinusMult(0)
-       ,fHistPiZeroMult(0)
-       ,fHistPiPlusMultAcc(0)
-       ,fHistPiMinusMultAcc(0)
-       ,fHistPiZeroMultAcc(0)
-       ,fPiPlusMult(0)
-       ,fPiMinusMult(0)
-       ,fPiZeroMult(0)
-       ,fPiPlusMultAcc(0)
-       ,fPiMinusMultAcc(0)
-       ,fPiZeroMultAcc(0)
-        ,fNeutralRemoved(0)
-        ,fChargedRemoved(0)
-        ,fChargedNotRemoved(0)
-        ,fNeutralNotRemoved(0)
-        ,fEnergyNeutralRemoved(0)
-        ,fEnergyChargedRemoved(0)
-        ,fEnergyChargedNotRemoved(0)
-        ,fEnergyNeutralNotRemoved(0)
-       ,fNClusters(0)
-       ,fGeoUtils(0)
-       ,fBadMapM2(0)
-       ,fBadMapM3(0)
-       ,fBadMapM4(0)
+    ,fImpactParameter(0)
+    ,fNcoll(0)
+    ,fNpart(0)
+    ,fPrimaryTree(0)
+    ,fTotEtWithSecondaryRemoved(0)
+    ,fTotEtSecondaryFromEmEtPrimary(0)
+    ,fTotEtSecondary(0)
+    ,fPrimaryCode(0)
+    ,fPrimaryCharge(0)
+    ,fPrimaryE(0)
+    ,fPrimaryEt(0)
+    ,fPrimaryPx(0)
+    ,fPrimaryPy(0)
+    ,fPrimaryPz(0)
+    ,fPrimaryVx(0)
+    ,fPrimaryVy(0)
+    ,fPrimaryVz(0)
+    ,fPrimaryAccepted(0)
+    ,fDepositedCode(0)
+    ,fDepositedEt(0)
+    ,fDepositedCharge(0)
+    ,fDepositedVx(0)
+    ,fDepositedVy(0)
+    ,fDepositedVz(0)
+    ,fHistDecayVertexNonRemovedCharged(0)
+    ,fHistDecayVertexRemovedCharged(0)
+    ,fHistDecayVertexNonRemovedNeutral(0)
+    ,fHistDecayVertexRemovedNeutral(0)
+
+    ,fHistRemovedOrNot(0)
+    ,fHistEtNonRemovedProtons(0)
+    ,fHistEtNonRemovedAntiProtons(0)
+    ,fHistEtNonRemovedPiPlus(0)
+    ,fHistEtNonRemovedPiMinus(0)
+    ,fHistEtNonRemovedKaonPlus(0)
+    ,fHistEtNonRemovedKaonMinus(0)
+    ,fHistEtNonRemovedK0s(0)
+    ,fHistEtNonRemovedK0L(0)
+    ,fHistEtNonRemovedLambdas(0)
+    ,fHistEtNonRemovedElectrons(0)
+    ,fHistEtNonRemovedPositrons(0)
+    ,fHistEtNonRemovedMuPlus(0)
+    ,fHistEtNonRemovedMuMinus(0)
+    ,fHistEtNonRemovedNeutrons(0)
+    ,fHistEtNonRemovedAntiNeutrons(0)
+    ,fHistEtNonRemovedGammas(0)
+    ,fHistEtNonRemovedGammasFromPi0(0)
+    ,fHistEtRemovedGammas(0)
+    ,fHistEtRemovedNeutrons(0)
+    ,fHistEtRemovedAntiNeutrons(0)
+    ,fHistEtRemovedCharged(0)
+    ,fHistEtRemovedNeutrals(0)
+    ,fHistEtNonRemovedCharged(0)
+    ,fHistEtNonRemovedNeutrals(0)
+    ,fHistMultNonRemovedProtons(0)
+    ,fHistMultNonRemovedAntiProtons(0)
+    ,fHistMultNonRemovedPiPlus(0)
+    ,fHistMultNonRemovedPiMinus(0)
+    ,fHistMultNonRemovedKaonPlus(0)
+    ,fHistMultNonRemovedKaonMinus(0)
+    ,fHistMultNonRemovedK0s(0)
+    ,fHistMultNonRemovedK0L(0)
+    ,fHistMultNonRemovedLambdas(0)
+    ,fHistMultNonRemovedElectrons(0)
+    ,fHistMultNonRemovedPositrons(0)
+    ,fHistMultNonRemovedMuPlus(0)
+    ,fHistMultNonRemovedMuMinus(0)
+    ,fHistMultNonRemovedNeutrons(0)
+    ,fHistMultNonRemovedAntiNeutrons(0)
+    ,fHistMultNonRemovedGammas(0)
+    ,fHistMultRemovedGammas(0)
+    ,fHistMultRemovedNeutrons(0)
+    ,fHistMultRemovedAntiNeutrons(0)
+    ,fHistMultRemovedCharged(0)
+    ,fHistMultRemovedNeutrals(0)
+    ,fHistMultNonRemovedCharged(0)
+    ,fHistMultNonRemovedNeutrals(0)
+    ,fHistTrackMultvsNonRemovedCharged(0)
+    ,fHistTrackMultvsNonRemovedNeutral(0)
+    ,fHistTrackMultvsRemovedGamma(0)
+    ,fHistClusterMultvsNonRemovedCharged(0)
+    ,fHistClusterMultvsNonRemovedNeutral(0)
+    ,fHistClusterMultvsRemovedGamma(0)
+    ,fHistMultvsNonRemovedChargedE(0)
+    ,fHistMultvsNonRemovedNeutralE(0)
+    ,fHistMultvsRemovedGammaE(0)
+    ,fEtNonRemovedProtons(0)
+    ,fEtNonRemovedAntiProtons(0)
+    ,fEtNonRemovedPiPlus(0)
+    ,fEtNonRemovedPiMinus(0)
+    ,fEtNonRemovedKaonPlus(0)
+    ,fEtNonRemovedKaonMinus(0)
+    ,fEtNonRemovedK0S(0)
+    ,fEtNonRemovedK0L(0)
+    ,fEtNonRemovedLambdas(0)
+    ,fEtNonRemovedElectrons(0)
+    ,fEtNonRemovedPositrons(0)
+    ,fEtNonRemovedMuMinus(0)
+    ,fEtNonRemovedMuPlus(0)
+    ,fEtNonRemovedGammas(0)
+    ,fEtNonRemovedGammasFromPi0(0)
+    ,fEtNonRemovedNeutrons(0)
+    ,fEtNonRemovedAntiNeutrons(0)
+    ,fEtRemovedProtons(0)
+    ,fEtRemovedAntiProtons(0)
+    ,fEtRemovedPiPlus(0)
+    ,fEtRemovedPiMinus(0)
+    ,fEtRemovedKaonPlus(0)
+    ,fEtRemovedKaonMinus(0)
+    ,fEtRemovedK0s(0)
+    ,fEtRemovedK0L(0)
+    ,fEtRemovedLambdas(0)
+    ,fEtRemovedElectrons(0)
+    ,fEtRemovedPositrons(0)
+    ,fEtRemovedMuMinus(0)
+    ,fEtRemovedMuPlus(0)
+    ,fEtRemovedGammasFromPi0(0)
+    ,fEtRemovedGammas(0)
+    ,fEtRemovedNeutrons(0)
+    ,fEtRemovedAntiNeutrons(0)
+    ,fMultNonRemovedProtons(0)
+    ,fMultNonRemovedAntiProtons(0)
+    ,fMultNonRemovedPiPlus(0)
+    ,fMultNonRemovedPiMinus(0)
+    ,fMultNonRemovedKaonPlus(0)
+    ,fMultNonRemovedKaonMinus(0)
+    ,fMultNonRemovedK0s(0)
+    ,fMultNonRemovedK0L(0)
+    ,fMultNonRemovedLambdas(0)
+    ,fMultNonRemovedElectrons(0)
+    ,fMultNonRemovedPositrons(0)
+    ,fMultNonRemovedMuMinus(0)
+    ,fMultNonRemovedMuPlus(0)
+    ,fMultNonRemovedGammas(0)
+    ,fMultNonRemovedNeutrons(0)
+    ,fMultNonRemovedAntiNeutrons(0)
+    ,fMultRemovedProtons(0)
+    ,fMultRemovedAntiProtons(0)
+    ,fMultRemovedPiPlus(0)
+    ,fMultRemovedPiMinus(0)
+    ,fMultRemovedKaonPlus(0)
+    ,fMultRemovedKaonMinus(0)
+    ,fMultRemovedK0s(0)
+    ,fMultRemovedK0L(0)
+    ,fMultRemovedLambdas(0)
+    ,fMultRemovedElectrons(0)
+    ,fMultRemovedPositrons(0)
+    ,fMultRemovedMuMinus(0)
+    ,fMultRemovedMuPlus(0)
+    ,fMultRemovedGammas(0)
+    ,fMultRemovedNeutrons(0)
+    ,fMultRemovedAntiNeutrons(0)
+    ,fTrackMultInAcc(0)
+    ,fHistDxDzNonRemovedCharged(0)
+    ,fHistDxDzRemovedCharged(0)
+    ,fHistDxDzNonRemovedNeutral(0)
+    ,fHistDxDzRemovedNeutral(0)
+    ,fHistPiPlusMult(0)
+    ,fHistPiMinusMult(0)
+    ,fHistPiZeroMult(0)
+    ,fHistPiPlusMultAcc(0)
+    ,fHistPiMinusMultAcc(0)
+    ,fHistPiZeroMultAcc(0)
+    ,fPiPlusMult(0)
+    ,fPiMinusMult(0)
+    ,fPiZeroMult(0)
+    ,fPiPlusMultAcc(0)
+    ,fPiMinusMultAcc(0)
+    ,fPiZeroMultAcc(0)
+    ,fNeutralRemoved(0)
+    ,fChargedRemoved(0)
+    ,fChargedNotRemoved(0)
+    ,fNeutralNotRemoved(0)
+    ,fEnergyNeutralRemoved(0)
+    ,fEnergyChargedRemoved(0)
+    ,fEnergyChargedNotRemoved(0)
+    ,fEnergyNeutralNotRemoved(0)
+    ,fNClusters(0)
+    ,fTotNeutralEtAfterMinEnergyCut(0)
 
 {
-    fTrackDistanceCut = 7.0;
 }
 
 // dtor
 AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
-{//Destructor
-       delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
-       delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
-       delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
-       delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
-       
-       delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
-       
-       delete fHistEtNonRemovedProtons; // enter comment here
-       delete fHistEtNonRemovedAntiProtons; // enter comment here
-       delete fHistEtNonRemovedPiPlus; // enter comment here
-       delete fHistEtNonRemovedPiMinus; // enter comment here
-       delete fHistEtNonRemovedKaonPlus; // enter comment here
-       delete fHistEtNonRemovedKaonMinus; // enter comment here
-       delete fHistEtNonRemovedK0s; // enter comment here
-       delete fHistEtNonRemovedK0L; // enter comment here
-       delete fHistEtNonRemovedLambdas; // enter comment here
-       delete fHistEtNonRemovedElectrons; // enter comment here
-       delete fHistEtNonRemovedPositrons; // enter comment here
-       delete fHistEtNonRemovedMuPlus; // enter comment here
-       delete fHistEtNonRemovedMuMinus; // enter comment here
-       delete fHistEtNonRemovedNeutrons; // enter comment here
-       delete fHistEtNonRemovedAntiNeutrons; // enter comment here
-       delete fHistEtNonRemovedGammas; // enter comment here
-       delete fHistEtNonRemovedGammasFromPi0; // enter comment here
-       
-       delete fHistEtRemovedGammas; // enter comment here
-       delete fHistEtRemovedNeutrons; // enter comment here
-       delete fHistEtRemovedAntiNeutrons; // enter comment here
-       
-       
-       delete fHistMultNonRemovedProtons; // enter comment here 
-       delete fHistMultNonRemovedAntiProtons; // enter comment here 
-       delete fHistMultNonRemovedPiPlus; // enter comment here 
-       delete fHistMultNonRemovedPiMinus; // enter comment here 
-       delete fHistMultNonRemovedKaonPlus; // enter comment here 
-       delete fHistMultNonRemovedKaonMinus; // enter comment here 
-       delete fHistMultNonRemovedK0s; // enter comment here 
-       delete fHistMultNonRemovedK0L; // enter comment here 
-       delete fHistMultNonRemovedLambdas; // enter comment here 
-       delete fHistMultNonRemovedElectrons; // enter comment here 
-       delete fHistMultNonRemovedPositrons; // enter comment here 
-       delete fHistMultNonRemovedMuPlus; // enter comment here
-       delete fHistMultNonRemovedMuMinus; // enter comment here
-       delete fHistMultNonRemovedNeutrons; // enter comment here
-       delete fHistMultNonRemovedAntiNeutrons; // enter comment here
-       delete fHistMultNonRemovedGammas; // enter comment here
-       
-       delete fHistMultRemovedGammas; // enter comment here
-       delete fHistMultRemovedNeutrons; // enter comment here
-       delete fHistMultRemovedAntiNeutrons; // enter comment here
-       
-       delete fHistTrackMultvsNonRemovedCharged; // enter comment here
-       delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
-       delete fHistTrackMultvsRemovedGamma; // enter comment here
-       
-       delete fHistClusterMultvsNonRemovedCharged; // enter comment here
-       delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
-       delete fHistClusterMultvsRemovedGamma; // enter comment here
-       
-       delete fHistMultvsNonRemovedChargedE; // enter comment here
-       delete fHistMultvsNonRemovedNeutralE; // enter comment here
-       delete fHistMultvsRemovedGammaE; // enter comment here
-       delete fHistDxDzNonRemovedCharged; // enter comment here
-       delete fHistDxDzRemovedCharged; // enter comment here
-       delete fHistDxDzNonRemovedNeutral; // enter comment here
-       delete fHistDxDzRemovedNeutral; // enter comment here
-       
-       delete fHistPiPlusMult; // enter comment here
-       delete fHistPiMinusMult; // enter comment here
-       delete fHistPiZeroMult; // enter comment here
-
-       delete fHistPiPlusMultAcc; // enter comment here
-       delete fHistPiMinusMultAcc; // enter comment here
-       delete fHistPiZeroMultAcc; // enter comment here
+{   //Destructor
+    delete fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
+    delete fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
+    delete fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
+    delete fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
+
+    delete fHistRemovedOrNot; // If charged/neutral particles were removed or not
+
+    delete fHistEtNonRemovedProtons; // enter comment here
+    delete fHistEtNonRemovedAntiProtons; // enter comment here
+    delete fHistEtNonRemovedPiPlus; // enter comment here
+    delete fHistEtNonRemovedPiMinus; // enter comment here
+    delete fHistEtNonRemovedKaonPlus; // enter comment here
+    delete fHistEtNonRemovedKaonMinus; // enter comment here
+    delete fHistEtNonRemovedK0s; // enter comment here
+    delete fHistEtNonRemovedK0L; // enter comment here
+    delete fHistEtNonRemovedLambdas; // enter comment here
+    delete fHistEtNonRemovedElectrons; // enter comment here
+    delete fHistEtNonRemovedPositrons; // enter comment here
+    delete fHistEtNonRemovedMuPlus; // enter comment here
+    delete fHistEtNonRemovedMuMinus; // enter comment here
+    delete fHistEtNonRemovedNeutrons; // enter comment here
+    delete fHistEtNonRemovedAntiNeutrons; // enter comment here
+    delete fHistEtNonRemovedGammas; // enter comment here
+    delete fHistEtNonRemovedGammasFromPi0; // enter comment here
+
+    delete fHistEtRemovedGammas; // enter comment here
+    delete fHistEtRemovedNeutrons; // enter comment here
+    delete fHistEtRemovedAntiNeutrons; // enter comment here
+
+
+    delete fHistMultNonRemovedProtons; // enter comment here
+    delete fHistMultNonRemovedAntiProtons; // enter comment here
+    delete fHistMultNonRemovedPiPlus; // enter comment here
+    delete fHistMultNonRemovedPiMinus; // enter comment here
+    delete fHistMultNonRemovedKaonPlus; // enter comment here
+    delete fHistMultNonRemovedKaonMinus; // enter comment here
+    delete fHistMultNonRemovedK0s; // enter comment here
+    delete fHistMultNonRemovedK0L; // enter comment here
+    delete fHistMultNonRemovedLambdas; // enter comment here
+    delete fHistMultNonRemovedElectrons; // enter comment here
+    delete fHistMultNonRemovedPositrons; // enter comment here
+    delete fHistMultNonRemovedMuPlus; // enter comment here
+    delete fHistMultNonRemovedMuMinus; // enter comment here
+    delete fHistMultNonRemovedNeutrons; // enter comment here
+    delete fHistMultNonRemovedAntiNeutrons; // enter comment here
+    delete fHistMultNonRemovedGammas; // enter comment here
+
+    delete fHistMultRemovedGammas; // enter comment here
+    delete fHistMultRemovedNeutrons; // enter comment here
+    delete fHistMultRemovedAntiNeutrons; // enter comment here
+
+    delete fHistTrackMultvsNonRemovedCharged; // enter comment here
+    delete fHistTrackMultvsNonRemovedNeutral; // enter comment here
+    delete fHistTrackMultvsRemovedGamma; // enter comment here
+
+    delete fHistClusterMultvsNonRemovedCharged; // enter comment here
+    delete fHistClusterMultvsNonRemovedNeutral; // enter comment here
+    delete fHistClusterMultvsRemovedGamma; // enter comment here
+
+    delete fHistMultvsNonRemovedChargedE; // enter comment here
+    delete fHistMultvsNonRemovedNeutralE; // enter comment here
+    delete fHistMultvsRemovedGammaE; // enter comment here
+    delete fHistDxDzNonRemovedCharged; // enter comment here
+    delete fHistDxDzRemovedCharged; // enter comment here
+    delete fHistDxDzNonRemovedNeutral; // enter comment here
+    delete fHistDxDzRemovedNeutral; // enter comment here
+
+    delete fHistPiPlusMult; // enter comment here
+    delete fHistPiMinusMult; // enter comment here
+    delete fHistPiZeroMult; // enter comment here
+
+    delete fHistPiPlusMultAcc; // enter comment here
+    delete fHistPiMinusMultAcc; // enter comment here
+    delete fHistPiZeroMultAcc; // enter comment here
 }
 
 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
-{ // analyse MC event
+{   // analyse MC event
     ResetEventValues();
 
     fPiPlusMult = 0;
@@ -305,14 +310,6 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
         fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
 
     }
-    fSparseTracks[0] = 0;
-    fSparseTracks[1] = 0;
-    fSparseTracks[2] = 0;
-    fSparseTracks[3] = 0;
-    fSparseTracks[4] = 0;
-    fSparseTracks[5] = 0;
-    fSparseTracks[6] = fCentClass;
-
 
     // Get us an mc event
     if (!ev) {
@@ -338,15 +335,8 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
         fImpactParameter = hijingGenHeader->ImpactParameter();
         fNcoll = hijingGenHeader->HardScatters(); // or should this be some combination of NN() NNw() NwN() NwNw() ?
         fNpart = hijingGenHeader->ProjectileParticipants() + hijingGenHeader->TargetParticipants();
-      }
-
-    /* // placeholder if we want to get some Pythia info later
-       AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
-       if (pythiaGenHeader) { // not Hijing; try with Pythia
-       printf("Pythia: ProcessType %d  GetPtHard %g \n",
-       pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
-       }
-    */
+    }
+
     // Let's play with the stack!
     AliStack *stack = event->Stack();
 
@@ -357,8 +347,6 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
     {
 
         TParticle *part = stack->Particle(iPart);
-//         TParticle *partMom = 0;
-//         TParticle *partMomLast = 0;
 
         if (!part)
         {
@@ -366,29 +354,28 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
             continue;
         }
 
-//         Int_t iPartMom = part->GetMother(0);
-//         Int_t iPartLastMom = part->GetMother(1);
-
         TParticlePDG *pdg = part->GetPDG(0);
 
-       if (!pdg)
+        if (!pdg)
         {
-            //Printf("ERROR: Could not get particle PDG %d", iPart);
+            Printf("ERROR: Could not get particle PDG %d", iPart);
             continue;
         }
 
         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
-       Int_t code = pdg->PdgCode();
+        Int_t code = pdg->PdgCode();
 
         // Check for reasonable (for now neutral and singly charged) charge on the particle
         //TODO:Maybe not only singly charged?
-        if ((stack->IsPhysicalPrimary(iPart))&&( TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 || TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3))
+        if (fSelector->CutNeutralMcParticle(iPart,*stack,*pdg))
         {
 
             fMultiplicity++;
-
+            //PrintFamilyTree(iPart, stack);
+//
             if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
             {
+                //Printf("Particle with eta: %f, pid: %d", part->Eta(), code);
                 // calculate E_T
                 if (
                     TMath::Abs(code) == fgProtonCode ||
@@ -408,120 +395,102 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                 }
                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 
-                fSparseTracks[0] = code;
-                fSparseTracks[1] = pdg->Charge()*3;
-                fSparseTracks[2] = pdg->Mass();
-                fSparseTracks[3] = et;
-                fSparseTracks[4] = part->Pt();
-                fSparseTracks[5] = part->Eta();
-                fSparseHistTracks->Fill(fSparseTracks);
 
                 // Fill up total E_T counters for each particle species
                 if (code == fgProtonCode || code == fgAntiProtonCode)
                 {
-                    fProtonEt += et;
                 }
                 if (code == fgPiPlusCode || code == fgPiMinusCode)
                 {
-                    fPionEt += et;
                     if (code == fgPiPlusCode)
                     {
-                        fPiPlusMult++;
                     }
                     else
                     {
-                        fPiMinusMult++;
                     }
                 }
                 if (code == fgGammaCode)
                 {
-                    fPiZeroMult++;
                 }
-                if (code == fgKPlusCode || code == fgKMinusCode)
+                if (code == fgKPlusCode)
+                {
+                }
+                if(code == fgKMinusCode)
                 {
-                    fChargedKaonEt += et;
                 }
                 if (code == fgMuPlusCode || code == fgMuMinusCode)
                 {
-                    fMuonEt += et;
                 }
                 if (code == fgEPlusCode || code == fgEMinusCode)
                 {
-                    fElectronEt += et;
                 }
                 // some neutrals also
                 if (code == fgNeutronCode)
                 {
-                    fNeutronEt += et;
                 }
                 if (code == fgAntiNeutronCode)
                 {
-                    fAntiNeutronEt += et;
                 }
                 if (code == fgGammaCode)
                 {
-                    fGammaEt += et;
                 }
 
                 // Neutral particles
                 //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
-                
+
                 if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
                 {
 
-                    if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
+                    //Printf("Gamma, phi: %f, eta: %f, phi cut min: %f, phi cut max: %f, eta cut: %f", part->Phi(), part->Eta(), fPhiMinCutAcc, fPhiMaxCutAcc, fEtaCutAcc);
+                    //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
 
                     // inside EMCal acceptance
-                    if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
+
+                    //if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiMaxCutAcc && part->Phi() > fPhiMinCutAcc)
+
+                    if(fSelector->CutGeometricalAcceptance(*part))
                     {
                         fNeutralMultiplicity++;
-                        //std::cout << code << ", " << iPart << ", " <<  part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
-
-                        //if(part->GetDaughter(0) > 0) std::cout << stack->Particle(part->GetDaughter(0))->GetPdgCode() << " " << stack->Particle(part->GetDaughter(1))->GetPdgCode() << std::endl;
-                       fTotNeutralEtAcc += et;
-                       fTotEtAcc += et;
-                        //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
-                        //if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
-                       
+                        fTotNeutralEt += et;
+                        if(fSelector->CutMinEnergy(*part))
+                        {
+                            fTotNeutralEtAfterMinEnergyCut += et;
+                        }
                         if (part->Energy() > 0.05) partCount++;
                     }
                 }
                 //Charged particles
                 else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
                 {
-                    fChargedMultiplicity++;
-                    fTotChargedEt += et;
 
                     // inside EMCal acceptance
-                    if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
+                    if (fSelector->CutGeometricalAcceptance(*part))
                     {
-                        fTotChargedEtAcc += et;
-                        fTotEtAcc += et;
+
+                        fChargedMultiplicity++;
+
+                        fTotChargedEt += et;
 
                         if (code == fgProtonCode || code == fgAntiProtonCode)
                         {
-                            fProtonEtAcc += et;
                         }
                         if (code == fgPiPlusCode || code == fgPiMinusCode)
                         {
-                            fPionEtAcc += et;
                         }
                         if (code == fgKPlusCode || code == fgKMinusCode)
                         {
-                            fChargedKaonEtAcc += et;
                         }
                         if (code == fgMuPlusCode || code == fgMuMinusCode)
                         {
-                            fMuonEtAcc += et;
                         }
-                        
+
                         if (code == fgEPlusCode || code == fgEMinusCode)
                         {
-                          fElectronEtAcc += et;
+                            fTotNeutralEt += et; // calling electrons neutral
+                            fTotChargedEt -= et;
                         }
                     } // inside EMCal acceptance
 
-                    //   if (TrackHitsCalorimeter(part, event->GetMagneticField()))
                     if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
                     {
                         if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
@@ -534,43 +503,19 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
     }
 
     fTotEt = fTotChargedEt + fTotNeutralEt;
-    fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
+    //fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;//
     //std::cout << "Event done! # of particles: " << partCount << std::endl;
     return 0;
 }
 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
-{ // analyse MC and real event info
+{   // analyse MC and real event info
     //if(!mcEvent || !realEvent){
     if (!ev || !ev2) {
         AliFatal("ERROR: Event does not exist");
         return 0;
     }
-    //if (!fMatrixInitialized)
-    if(0)
-    {
-        for (Int_t mod=0; mod<5; mod++) {
-            if (!ev->GetPHOSMatrix(mod)) continue;
-            fGeoUtils->SetMisalMatrix(ev->GetPHOSMatrix(mod),mod) ;
-            Printf("PHOS geo matrix %p for module # %d is set\n", ev->GetPHOSMatrix(mod), mod);
-
-        }
-        fMatrixInitialized = kTRUE;
-    }
-    fSparseClusters[0] = 0;
-    fSparseClusters[1] = 0;
-    fSparseClusters[2] = 0;
-    fSparseClusters[3] = 0;
-    fSparseClusters[4] = 0;
-    fSparseClusters[5] = 0;
-    fSparseClusters[6] = 0;
-    fSparseClusters[7] = 0;
-    fSparseClusters[8] = 0;
-    fSparseClusters[9] = fCentClass;
-    fSparseClusters[10] = 0;
-
-    //AnalyseEvent(mcEvent);
-    AnalyseEvent(ev);
+    AliAnalysisEt::AnalyseEvent(ev);
     AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
     AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
     if (!mcEvent || !realEvent) {
@@ -578,50 +523,38 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
         return 0;
     }
 
+    fSelector->SetEvent(realEvent);
+
+    AnalyseEvent(ev);
+
     AliStack *stack = mcEvent->Stack();
 
     // get all detector clusters
-    TRefArray* caloClusters = new TRefArray();
-    if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
-    else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
-    else {
-        AliFatal("Detector ID has not been specified");
-        return -1;
-    }
-if(0)
-{
-    // Track loop to fill a pT spectrum
-    for (Int_t iTracks = 0; iTracks < realEvent->GetNumberOfTracks(); iTracks++)
-    {
-        AliESDtrack* track = realEvent->GetTrack(iTracks);
-        if (!track)
-        {
-            Printf("ERROR: Could not receive track %d", iTracks);
-            continue;
-        }
+    //  TRefArray* caloClusters = new TRefArray();
+
+//    if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
+    //else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
+    //else {
+    //AliFatal("Detector ID has not been specified");
+    //return -1;
+//    }
+
+    TRefArray *caloClusters = fSelector->GetClusters();
 
-        if ( track->Phi() < fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180. && track->Phi() > fCuts->GetGeometryPhosPhiAccMinCut()  * TMath::Pi()/180. && TMath::Abs(track->Eta()) < fCuts->GetGeometryPhosEtaAccCut())
-        {
-            fTrackMultInAcc++;
-        }
-    }
-}
     Int_t nCluster = caloClusters->GetEntries();
     fNClusters = 0;
     // loop the clusters
     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
     {
-       
+        Int_t cf = 0;
         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
         //Float_t caloE = caloCluster->E();
-       
-       if(caloCluster->GetType() != fClusterType) continue;
-       if(AliAnalysisEtMonteCarlo::TooCloseToBadChannel(*caloCluster))continue;
+        fCutFlow->Fill(cf++);
+        if(!fSelector->CutDistanceToBadChannel(*caloCluster)) continue;
 
-       fNClusters++;
+        fNClusters++;
         UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
         TParticle *part  =  stack->Particle(iPart);
-        //TParticle *partMom = 0;
 
         if (!part)
         {
@@ -629,74 +562,64 @@ if(0)
             continue;
         }
 
-        // compare MC and Rec energies for all particles
-        //fHistAllERecEMC->Fill(part->Energy(),caloE);
-
-        //TParticlePDG *pdgMom = 0;
-
-/*
-        Int_t iPartMom = part->GetMother(0);
-
-        if (iPartMom>0)
-        {
-            partMom = stack->Particle(iPartMom);
-            pdgMom = partMom->GetPDG(0);
-        }*/
-       //int mcClusterCharge = stack->Particle(iPart)->GetPDG()->Charge();
+        int primIdx = iPart;
         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
         {
-   //         if (!((code == fgEPlusCode) || (code == fgEMinusCode) || (code == fgGammaCode)))
-               int p = GetPrimMother(iPart, stack);
-       //      std::cout << p << std::endl;
-               if(p > 0) 
-               { 
-                 part = stack->Particle(p);
-               }
-               else 
-               {
-                 std::cout << "What!? No primary?" << std::endl;
-                 continue;
-               }
-               
+            primIdx = GetPrimMother(iPart, stack);
+            if(primIdx < 0)
+            {
+                std::cout << "What!? No primary?" << std::endl;
+                continue;
+            }
+
         } // end of primary particle check
+        const int primCode = stack->Particle(primIdx)->GetPdgCode();
         TParticlePDG *pdg = part->GetPDG();
-       Int_t code = pdg->PdgCode();
-       if (!pdg)
+        if (!pdg)
         {
             Printf("ERROR: Could not get particle PDG %d", iPart);
             continue;
-        }      
-        //if(code == fgK0SCode)PrintFamilyTree(iPart, stack);
-       //pdgMom = stack->Particle(GetPrimMother(iPart, stack))->GetPDG();
-       
-        //-if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
+        }
 
+        Int_t code = pdg->PdgCode();
         Double_t clEt = CalculateTransverseEnergy(caloCluster);
 //     if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
-        if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
+        if(!fSelector->CutMinEnergy(*caloCluster)) continue;
+        fCutFlow->Fill(cf++);
         Float_t pos[3];
 
         caloCluster->GetPosition(pos);
         TVector3 cp(pos);
 
-        fSparseClusters[0] = code;
-        fSparseClusters[1] = pdg->Charge()/3;
-        fSparseClusters[2] = pdg->Mass();
-        fSparseClusters[3] = clEt;
-        fSparseClusters[4] = caloCluster->E();
-        fSparseClusters[5] = cp.Eta();
-        fSparseClusters[6] = part->Energy() * TMath::Sin(part->Theta());
-        fSparseClusters[7] = part->Pt();
-        fSparseClusters[8] = part->Eta();
-        fSparseClusters[9] = fCentClass;
-        fSparseClusters[10] = 0;
-        fSparseHistClusters->Fill(fSparseClusters);
-
-       
+        TParticle *primPart = stack->Particle(primIdx);
+        fPrimaryCode = primPart->GetPdgCode();
+        fPrimaryCharge = primPart->GetPDG()->Charge();
+
+        fPrimaryE = primPart->Energy();
+        fPrimaryEt = primPart->Energy()*TMath::Sin(primPart->Theta());
+
+        fPrimaryPx = primPart->Px();
+        fPrimaryPy = primPart->Py();
+        fPrimaryPz = primPart->Pz();
+
+        fPrimaryVx = primPart->Vx();
+        fPrimaryVy = primPart->Vy();
+        fPrimaryVz = primPart->Vz();
+
+        fPrimaryAccepted = false;
+
+        fDepositedCode = part->GetPdgCode();
+        fDepositedEt = caloCluster->E() * TMath::Sin(cp.Eta());
+        fDepositedCharge = part->GetPDG()->Charge();
+
+        fDepositedVx = part->Vx();
+        fDepositedVy = part->Vy();
+        fDepositedVz = part->Vz();
         //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
-       pdg = part->GetPDG(0);
+        pdg = part->GetPDG(0);
         //if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
-       if(caloCluster->GetTrackMatchedIndex() > -1 && (fCuts->GetPhosTrackRCut() > TestCPV(caloCluster->GetTrackDx(), caloCluster->GetTrackDz(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Pt(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Charge(), ev)))
+        if(!fSelector->CutTrackMatching(*caloCluster))
+            //if(caloCluster->GetTrackMatchedIndex() > -1 && (fCuts->GetPhosTrackRCut() > TestCPV(caloCluster->GetTrackDx(), caloCluster->GetTrackDz(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Pt(), ev->GetTrack(caloCluster->GetTrackMatchedIndex())->Charge(), ev)))
         {
 
             if (pdg->Charge() != 0)
@@ -704,12 +627,91 @@ if(0)
                 //std::cout << "Removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
                 fChargedRemoved++;
                 //fEnergyChargedRemoved += caloCluster->E();
-               fEnergyChargedRemoved += clEt;
+                fEnergyChargedRemoved += clEt;
                 fHistRemovedOrNot->Fill(0.0, fCentClass);
                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
                 //fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
                 fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+                if(code == fgProtonCode)
+                {
+                    fProtonRemovedEt += clEt;
+                    fProtonRemovedMult++;
 
+                }
+                else if(code == fgAntiProtonCode)
+                {
+                    fAntiProtonRemovedEt += clEt;
+                    fAntiProtonRemovedMult++;
+                }
+                else if(code == fgPiPlusCode)
+                {
+                    fPiPlusRemovedEt += clEt;
+                    fPiPlusRemovedMult++;
+                }
+                else if(code == fgPiMinusCode)
+                {
+                    fPiMinusRemovedEt += clEt;
+                    fPiMinusRemovedMult++;
+                }
+                else if(code == fgKPlusCode)
+                {
+                    fKPlusRemovedEt += clEt;
+                    fKPlusRemovedMult++;
+                }
+                else if(code == fgKMinusCode)
+                {
+                    fKMinusRemovedEt += clEt;
+                    fKMinusRemovedMult++;
+                }
+                else if(code == fgMuMinusCode)
+                {
+                    fMuMinusRemovedEt += clEt;
+                    fMuMinusRemovedMult++;
+                }
+                else if(code == fgMuPlusCode)
+                {
+                    fMuPlusRemovedEt += clEt;
+                    fMuPlusRemovedMult++;
+                }
+                else if(code == fgEMinusCode)
+                {
+                    fEMinusRemovedEt += clEt;
+                    fEMinusRemovedMult++;
+                }
+                else if(code == fgEPlusCode)
+                {
+                    fEPlusRemovedEt += clEt;
+                    fEPlusRemovedMult++;
+                }
+                else
+                {
+                    Printf("Removed: %d, with E_T: %f", code, clEt);
+                }
+                if (primCode == fgGammaCode)
+                {
+                    fGammaRemovedEt += clEt;
+                    fGammaRemovedMult++;
+                }
+                else if (primCode == fgPi0Code)
+                {
+                    fPi0RemovedEt += clEt;
+                    fPi0RemovedMult++;
+                }
+                else if (primCode == fgEtaCode)
+                {
+                    fPi0RemovedEt += clEt;
+                    fPi0RemovedMult++;
+                }
+                else if(primCode == fgK0LCode)
+                {
+                    fK0lRemovedEt += clEt;
+                    fK0lRemovedMult++;
+                }
+                else if(primCode == fgK0SCode)
+                {
+                    fK0sRemovedEt += clEt;
+                    fK0sRemovedMult++;
+                }
 
             }
             else
@@ -718,7 +720,7 @@ if(0)
                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
                 fNeutralRemoved++;
                 //fEnergyNeutralRemoved += caloCluster->E();
-               fEnergyNeutralRemoved += clEt;
+                fEnergyNeutralRemoved += clEt;
                 fHistRemovedOrNot->Fill(1.0, fCentClass);
                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
                 //fHistDecayVertexRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
@@ -733,26 +735,53 @@ if(0)
                 {
                     fEtRemovedNeutrons += clEt;
                     fMultRemovedNeutrons++;
+                    fNeutronRemovedEt += clEt;
+                    fNeutronRemovedMult++;
                 }
                 else if (code == fgAntiNeutronCode)
                 {
                     fEtRemovedAntiNeutrons += clEt;
                     fMultRemovedAntiNeutrons++;
+                    fAntiNeutronRemovedEt += clEt;
+                    fAntiNeutronRemovedMult++;
+                }
+                if (primCode == fgGammaCode)
+                {
+                    fGammaRemovedEt += clEt;
+                    fGammaRemovedMult++;
+                }
+                else if(primCode == fgPi0Code)
+                {
+                    fPi0RemovedEt += clEt;
+                    fPi0RemovedMult++;
+                }
+                else if(primCode == fgEtaCode)
+                {
+                    fPi0RemovedEt += clEt;
+                    fPi0RemovedMult++;
+                }
+                else if(primCode == fgK0LCode)
+                {
+                    fK0lRemovedEt += clEt;
+                    fK0lRemovedMult++;
+                }
+                else if(primCode == fgK0SCode)
+                {
+                    fK0sRemovedEt += clEt;
+                    fK0sRemovedMult++;
                 }
-                
-                //else std::cout << "Hmm, what is this (neutral: " << code << std::endl;
             }
         }
         else
         {
-
+            fCutFlow->Fill(cf++);
             if (pdg->Charge() != 0)
             {
-                std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                //std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
                 fChargedNotRemoved++;
                 //fEnergyChargedNotRemoved += caloCluster->E();
-               fEnergyChargedNotRemoved += clEt;
+                fEnergyChargedNotRemoved += clEt;
                 fHistRemovedOrNot->Fill(2.0, fCentClass);
                 //std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
@@ -763,104 +792,162 @@ if(0)
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedProtons += clEt;
                     fMultNonRemovedProtons++;
+                    fProtonEt += clEt;
+                    fProtonMult++;
                 }
                 else if (code == fgAntiProtonCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedAntiProtons += clEt;
                     fMultNonRemovedAntiProtons++;
+                    fAntiProtonEt += clEt;
+                    fAntiProtonMult++;
                 }
                 else if (code == fgPiPlusCode)
                 {
                     //std::cout << "PI+" <<  clEt << std::endl;
                     fEtNonRemovedPiPlus += clEt;
                     fMultNonRemovedPiPlus++;
+                    fPiPlusEt += clEt;
+                    fPiPlusMult++;
                 }
                 else if (code == fgPiMinusCode)
                 {
                     // std::cout << "PI-"  << clEt << std::endl;
                     fEtNonRemovedPiMinus += clEt;
                     fMultNonRemovedPiMinus++;
+                    fPiMinusEt += clEt;
+                    fPiMinusMult++;
                 }
                 else if (code == fgKPlusCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedKaonPlus += clEt;
                     fMultNonRemovedKaonPlus++;
+                    fKPlusEt += clEt;
+                    fKPlusMult++;
                 }
                 else if (code == fgKMinusCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedKaonMinus += clEt;
                     fMultNonRemovedKaonMinus++;
+                    fKMinusEt += clEt;
+                    fKMinusMult++;
                 }
                 else if (code == fgEPlusCode)
                 {
                     //std::cout << clEt << std::endl;
-                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
+                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
                     {
                         fEtNonRemovedPositrons += clEt;
                         fMultNonRemovedPositrons++;
+                        fEPlusEt += clEt;
+                        fEPlusMult++;
                     }
                 }
                 else if (code == fgEMinusCode)
                 {
                     //std::cout << clEt << std::endl;
-                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
+                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 440)
                     {
                         fEtNonRemovedElectrons += clEt;
                         fMultNonRemovedElectrons++;
+                        fEMinusEt += clEt;
+                        fEMinusMult++;
                     }
                 }
                 else if (code == fgMuPlusCode)
                 {
                     fEtNonRemovedMuPlus += clEt;
                     fMultNonRemovedMuPlus++;
+                    fMuPlusEt += clEt;
+                    fMuPlusMult++;
                 }
                 else if (code == fgMuMinusCode)
                 {
                     fEtNonRemovedMuMinus += clEt;
                     fMultNonRemovedMuMinus++;
+                    fMuMinusEt += clEt;
+                    fMuMinusMult++;
+                }
+                if (primCode == fgGammaCode)
+                {
+                    fGammaEt += clEt;
+                    fGammaMult++;
+                }
+                else if(primCode == fgPi0Code)
+                {
+                    fPi0Et += clEt;
+                    fPi0Mult++;
+                }
+                else if(primCode == fgEtaCode)
+                {
+                    fPi0Et += clEt;
+                    fPi0Mult++;
+                }
+                else if(primCode == fgK0LCode)
+                {
+                    fK0lEt += clEt;
+                    fK0lMult++;
+                }
+                else if(primCode == fgK0SCode)
+                {
+                    fK0sEt += clEt;
+                    fK0sMult++;
                 }
 
             }
             else
             {
-                //std::cout << "Accepted: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                fPrimaryAccepted = true;
                 //std::cout << "Not removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
                 //std::cout << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
                 fNeutralNotRemoved++;
-               fEnergyNeutralNotRemoved += clEt;
+                fEnergyNeutralNotRemoved += clEt;
                 fHistRemovedOrNot->Fill(3.0, fCentClass);
                 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+
+                if(TMath::Abs(part->Vx()) < 1.0 && TMath::Abs(part->Vy()) < 1.0 && TMath::Abs(part->Vz()) < 20 && fSelector->IsEmEtParticle(primCode))
+                {
+                    fTotEtWithSecondaryRemoved += clEt;
+                }
+                else if(fSelector->IsEmEtParticle(primCode))
+                {
+                    fTotEtSecondaryFromEmEtPrimary += clEt;
+                }
+                else           {
+                    fTotEtSecondary += clEt;
+                }
+                //code = stack->Particle(primIdx)->GetPdgCode();
                 if (code == fgGammaCode)
                 {
                     fEtNonRemovedGammas += clEt;
                     fMultNonRemovedGammas++;
 //                    if (pdgMom)
-  //                  {
-    //                    if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
-      //                  {
+                    //                  {
+                    //                    if (TMath::Abs(pdgMom->PdgCode()) == fgPi0Code || TMath::Abs(pdgMom->PdgCode()) == fgEtaCode || TMath::Abs(pdgMom->PdgCode()) == 331)
+                    //                  {
 //                     std::cout << "Mother of gamma: " << pdgMom->PdgCode() << " " << pdgMom->GetName() <<  ", E: " << part->Energy() << std::endl;
-        //                    fEtNonRemovedGammasFromPi0 += clEt;
-          //              }
-            //        }
+                    //                    fEtNonRemovedGammasFromPi0 += clEt;
+                    //              }
+                    //        }
                 }
                 else if(TMath::Abs(code) == fgPi0Code)
-               {
-                 fEtNonRemovedGammasFromPi0 += clEt;
-                 fMultNonRemovedGammas++;
-               }
+                {
+                    fEtNonRemovedGammasFromPi0 += clEt;
+                    fMultNonRemovedGammas++;
+                }
                 else if(TMath::Abs(code) == fgEtaCode)
-               {
-                 fEtNonRemovedGammasFromPi0 += clEt;
-                 fMultNonRemovedGammas++;
-               }
+                {
+                    fEtNonRemovedGammasFromPi0 += clEt;
+                    fMultNonRemovedGammas++;
+                }
                 else if(TMath::Abs(code) == 331)
-               {
-                 fEtNonRemovedGammasFromPi0 += clEt;
-                 fMultNonRemovedGammas++;
-               }
+                {
+                    fEtNonRemovedGammasFromPi0 += clEt;
+                    fMultNonRemovedGammas++;
+                }
                 else if (code == fgNeutronCode)
                 {
                     fEtNonRemovedNeutrons += clEt;
@@ -874,30 +961,53 @@ if(0)
                 //else if (code == fgK0LCode || pdgMom->PdgCode() == fgK0SCode)
                 else if (code == fgK0SCode)
                 {
-                   //std::cout << "K0 with energy: " << clEt << std::endl;
+                    //std::cout << "K0 with energy: " << clEt << std::endl;
                     fEtNonRemovedK0S += clEt;
                     fMultNonRemovedK0s++;
-
                 }
                 else if(TMath::Abs(code) == fgK0LCode)
-               {
-                 
-                 fEtNonRemovedK0L += clEt;
-                 fMultNonRemovedK0L++;
-               }
-               
+                {
+
+                    fEtNonRemovedK0L += clEt;
+                    fMultNonRemovedK0L++;
+                }
+
                 else if (TMath::Abs(code) == fgLambdaCode)
                 {
                     fEtNonRemovedLambdas += clEt;
                     fMultNonRemovedLambdas++;
                 }
-                  
                 else std::cout << "Hmm, what is this (neutral not removed): " << code << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
+                if (primCode == fgGammaCode)
+                {
+                    fGammaEt += clEt;
+                    fGammaMult++;
+                }
+                else if(primCode == fgPi0Code)
+                {
+                    fPi0Et += clEt;
+                    fPi0Mult++;
+                }
+                else if(primCode == fgEtaCode)
+                {
+                    fPi0Et += clEt;
+                    fPi0Mult++;
+                }
+                else if(primCode == fgK0LCode)
+                {
+                    fK0lEt += clEt;
+                    fK0lMult++;
+                }
+                else if(primCode == fgK0SCode)
+                {
+                    fK0sEt += clEt;
+                    fK0sMult++;
+                }
             }
         }
+        fPrimaryTree->Fill();
     } // end of loop over clusters
 
-    delete caloClusters;
     //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
     //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
     //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
@@ -914,20 +1024,18 @@ if(0)
 }
 
 void AliAnalysisEtMonteCarlo::Init()
-{ // init
+{   // init
     AliAnalysisEt::Init();
-    TFile *f = TFile::Open("badchannels.root", "READ");
-
-    fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
-    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
-    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
-
 }
 
 void AliAnalysisEtMonteCarlo::ResetEventValues()
-{ // reset event values
+{   // reset event values
     AliAnalysisEt::ResetEventValues();
 
+    fTotEtSecondary = 0;
+    fTotEtSecondaryFromEmEtPrimary = 0;
+    fTotEtWithSecondaryRemoved = 0;
+
     // collision geometry defaults for p+p:
     fImpactParameter = 0;
     fNcoll = 1;
@@ -1010,26 +1118,32 @@ void AliAnalysisEtMonteCarlo::ResetEventValues()
     fEnergyChargedRemoved = 0;
     fEnergyNeutralNotRemoved = 0;
     fEnergyNeutralRemoved = 0;
-    
+
     fChargedNotRemoved = 0;
     fChargedRemoved = 0;
     fNeutralNotRemoved = 0;
     fNeutralRemoved = 0;
-    
-    
+
+
     fTrackMultInAcc = 0;
 
+    fTotNeutralEtAfterMinEnergyCut = 0;
+
 }
 
 void AliAnalysisEtMonteCarlo::CreateHistograms()
-{ // histogram related additions
+{   // histogram related additions
     AliAnalysisEt::CreateHistograms();
-    if (fTree) {
-        fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
-        fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
-        fTree->Branch("fNpart",&fNpart,"fNpart/I");
+    if (fEventSummaryTree) {
+        fEventSummaryTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
+        fEventSummaryTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
+        fEventSummaryTree->Branch("fNpart",&fNpart,"fNpart/I");
+        fEventSummaryTree->Branch("fTotEtWithSecondaryRemoved", &fTotEtWithSecondaryRemoved, "fTotEtWithSecondaryRemoved/D");
+        fEventSummaryTree->Branch("fTotEtSecondaryFromEmEtPrimary", &fTotEtSecondaryFromEmEtPrimary, "fTotEtSecondaryFromEmEtPrimary/D");
+        fEventSummaryTree->Branch("fTotEtSecondary", &fTotEtSecondary, "fTotEtSecondary/D");
+        fEventSummaryTree->Branch("fTotNeutralEtAfterMinEnergyCut", &fTotNeutralEtAfterMinEnergyCut, "fTotNeutralEtAfterMinEnergyCut/D");
     }
-       
+
     //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
     //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
     //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
@@ -1059,13 +1173,13 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistEtRemovedGammas = new  TH2F("fHistEtRemovedGammas", "fHistEtRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtRemovedNeutrons = new  TH2F("fHistEtRemovedNeutrons", "fHistEtRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtRemovedAntiNeutrons = new  TH2F("fHistEtRemovedAntiNeutrons", "fHistEtRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
-  
+
     fHistEtRemovedCharged = new  TH2F("fHistEtRemovedCharged", "fHistEtRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtRemovedNeutrals = new  TH2F("fHistEtRemovedNeutrals", "fHistEtRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
 
     fHistEtNonRemovedCharged = new  TH2F("fHistEtNonRemovedCharged", "fHistEtNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedNeutrals = new  TH2F("fHistEtNonRemovedNeutrals", "fHistEtNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
-  
+
     fHistMultNonRemovedProtons = new TH2F("fHistMultNonRemovedProtons", "fHistMultNonRemovedProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedAntiProtons = new TH2F("fHistMultNonRemovedAntiProtons", "fHistMultNonRemovedAntiProtons", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedPiPlus = new TH2F("fHistMultNonRemovedPiPlus", "fHistMultNonRemovedPiPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
@@ -1087,14 +1201,14 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
     fHistMultRemovedNeutrons = new  TH2F("fHistMultRemovedNeutrons", "fHistMultRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultRemovedAntiNeutrons = new  TH2F("fHistMultRemovedAntiNeutrons", "fHistMultRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
-  /*
-    fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
-    fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
+    /*
+      fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+      fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);
+
+      fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
+      fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
 
-    fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 1500, 0, 30, 10, -0.5, 9.5);
-    fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 1500, 0, 30, 10, -0.5, 9.5);*/
 
-  
     fHistMultRemovedCharged = new  TH2F("fHistMultRemovedCharged", "fHistMultRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
     fHistMultRemovedNeutrals = new  TH2F("fHistMultRemovedNeutrals", "fHistMultRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
 
@@ -1122,12 +1236,53 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
     fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
 
+    if(fCuts->GetHistMakeTree())
+    {
+        TString treename = "fPrimaryTree" + fHistogramNameSuffix;
+        fPrimaryTree = new TTree(treename, treename);
+
+        fPrimaryTree->Branch("fTotEt",&fTotEt,"fTotEt/D");
+        fPrimaryTree->Branch("fNeutralMultiplicity", &fNeutralMultiplicity, "fNeutralMultiplicity/I");
+        fPrimaryTree->Branch("fCentClass",&fCentClass,"fCentClass/I");
+
+        fPrimaryTree->Branch("fPrimaryCode", &fPrimaryCode, "fPrimaryCode/I");
+        fPrimaryTree->Branch("fPrimaryCharge", &fPrimaryCharge, "fPrimaryCharge/I");
+
+        fPrimaryTree->Branch("fPrimaryE", &fPrimaryE, "fPrimaryE/D");
+        fPrimaryTree->Branch("fPrimaryEt", &fPrimaryEt, "fPrimaryEt/D");
+
+        fPrimaryTree->Branch("fPrimaryPx", &fPrimaryPx, "fPrimaryPx/D");
+        fPrimaryTree->Branch("fPrimaryPy", &fPrimaryPy, "fPrimaryPy/D");
+        fPrimaryTree->Branch("fPrimaryPz", &fPrimaryPz, "fPrimaryPz/D");
+
+        fPrimaryTree->Branch("fPrimaryVx", &fPrimaryVx, "fPrimaryVx/D");
+        fPrimaryTree->Branch("fPrimaryVy", &fPrimaryVy, "fPrimaryVy/D");
+        fPrimaryTree->Branch("fPrimaryVz", &fPrimaryVz, "fPrimaryVz/D");
+
+        fPrimaryTree->Branch("fPrimaryAccepted", &fPrimaryAccepted, "fPrimaryAccepted/B");
+
+
+        fPrimaryTree->Branch("fDepositedCode", &fDepositedCode, "fDepositedCode/I");
+        fPrimaryTree->Branch("fDepositedCharge", &fDepositedCharge, "fDepositedCharge/I");
+        fPrimaryTree->Branch("fDepositedEt", &fDepositedEt, "fDepositedEt/D");
+
+        fPrimaryTree->Branch("fDepositedVx", &fDepositedVx, "fDepositedVx/D");
+        fPrimaryTree->Branch("fDepositedVy", &fDepositedVy, "fDepositedVy/D");
+        fPrimaryTree->Branch("fDepositedVz", &fDepositedVz, "fDepositedVz/D");
+
+    }
+
 }
 
 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
-{//fill the output list
+{   //fill the output list
     AliAnalysisEt::FillOutputList(list);
 
+    if(fCuts->GetHistMakeTree())
+    {
+        list->Add(fPrimaryTree);
+    }
+
     list->Add(fHistRemovedOrNot);
 
     list->Add(fHistEtNonRemovedProtons);
@@ -1223,17 +1378,14 @@ bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t mag
 
     // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
 
-    bool status = prop &&
-                  TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
-                  esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
-                  esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
+    bool status = prop && fSelector->CutGeometricalAcceptance(*esdTrack);
     delete esdTrack;
 
     return status;
 }
 
 void AliAnalysisEtMonteCarlo::FillHistograms()
-{ // let base class fill its histograms, and us fill the local ones
+{   // let base class fill its histograms, and us fill the local ones
     AliAnalysisEt::FillHistograms();
     //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
 
@@ -1264,7 +1416,7 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
 //                                             +fEtRemovedProtons.
 //                             fCentClass);
 //     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
-// 
+//
 //     fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
 //                                             +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
 //                                             +fEtNonRemovedProtons,
@@ -1275,13 +1427,13 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
     fHistEtRemovedNeutrals->Fill(fEnergyNeutralRemoved, fNClusters);
     fHistEtNonRemovedCharged->Fill(fEnergyChargedNotRemoved, fNClusters);
     fHistEtNonRemovedNeutrals->Fill(fEnergyNeutralNotRemoved, fNClusters);
-    
+
     fHistMultRemovedCharged->Fill(fChargedRemoved, fNClusters);
     fHistMultRemovedNeutrals->Fill(fNeutralRemoved, fNClusters);
     fHistMultNonRemovedCharged->Fill(fChargedNotRemoved, fNClusters);
     fHistMultNonRemovedNeutrals->Fill(fNeutralNotRemoved, fNClusters);
-    
-    
+
+
     fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
     fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
     fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
@@ -1332,210 +1484,104 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
 
 Int_t AliAnalysisEtMonteCarlo::PrintFamilyTree(Int_t partIdx, AliStack* stack)
 {
-  TParticle *part = stack->Particle(partIdx);
-  std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
-  std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
-  std::cout << "Energy: " << part->Energy() << std::endl;
-  std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
-  PrintMothers(partIdx, stack, 1);
-  return 0;
+    TParticle *part = stack->Particle(partIdx);
+    if(part->GetPdgCode() == fgK0SCode)
+    {
+        std::cout << "This is index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") , is it primary: " << stack->IsPhysicalPrimary(partIdx)<< std::endl;
+        std::cout << "PID: " << part->GetPdgCode() << "/" << part->GetName() << std::endl;
+        std::cout << "Energy: " << part->Energy() << std::endl;
+        std::cout << "Vertex: " << part->Vx() << ", " << part->Vy() << ", " << part->Vz() << std::endl;
+    }
+    return PrintMothers(partIdx, stack, 1);
 }
 
 Int_t AliAnalysisEtMonteCarlo::PrintMothers(Int_t partIdx, AliStack* stack, Int_t gen)
 {
-  char *tabs = new char[gen+1];
-  for(Int_t i = 0; i < gen; ++i)
-  {
-    //std::cout << i << std::endl;
-    tabs[i] = '\t';
-  }
-  tabs[gen] = '\0';
-  Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
-  if(mothIdx < 0) return 0;
-  TParticle *mother = stack->Particle(mothIdx);
-  std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
-  std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
-  std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
-  std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
-  if(mother->GetPdgCode() == fgK0SCode) 
-      {        
+    char *tabs = new char[gen+1];
+    for(Int_t i = 0; i < gen; ++i)
+    {
+        //std::cout << i << std::endl;
+        tabs[i] = '\t';
+    }
+    tabs[gen] = '\0';
+    Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+    if(mothIdx < 0)
+    {
+      delete [] tabs;
+      return 0;
+    }
+    TParticle *mother = stack->Particle(mothIdx);
+    if(mother->GetPdgCode() == fgK0SCode)
+    {
+        std::cout << tabs << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << ", is it primary: " << stack->IsPhysicalPrimary(mothIdx)<< std::endl;
+        std::cout << tabs << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
+        std::cout << tabs << "Energy: " << mother->Energy() << std::endl;
+        std::cout << tabs << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
+    }
+    if(mother->GetPdgCode() == fgK0SCode)
+    {
 //     std::cout << "K0S!!!!!!!!!!!!!11111!!!!!" << std::endl;
-      }
+    }
 //  std::cout << "Mother of index: " << partIdx << " (" << stack->Particle(partIdx)->GetName() <<") is: " << mothIdx << std::endl;
 //  std::cout << "PID: " << mother->GetPdgCode() << "/" << mother->GetName() << std::endl;
 //  std::cout << "Energy: " << mother->Energy() << std::endl;
 //  std::cout << "Vertex: " << mother->Vx() << ", " << mother->Vy() << ", " << mother->Vz() << std::endl;
-  
-  delete [] tabs;
-  return PrintMothers(mothIdx, stack, gen+1) + 1;
+
+    delete [] tabs;
+    return PrintMothers(mothIdx, stack, gen+1) + 1;
 }
 
 Int_t AliAnalysisEtMonteCarlo::GetPrimMother(Int_t partIdx, AliStack *stack)
 {
-  if(partIdx >= 0) 
-  {
-    Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
-    if(mothIdx < 0) return -1;
-    TParticle *mother = stack->Particle(mothIdx);
-    if(mother)
+    if(partIdx >= 0)
     {
-     // if(mother->GetPdgCode() == fgK0SCode) 
-      //{      
+        Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+        if(mothIdx < 0) return -1;
+        TParticle *mother = stack->Particle(mothIdx);
+        if(mother)
+        {
+            // if(mother->GetPdgCode() == fgK0SCode)
+            //{
 //     std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
-       //return mothIdx;
-  //    }
-      //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx)) 
-      //{      
+            //return mothIdx;
+            //    }
+            //if(mother->GetPdgCode() == fgK0SCode&& stack->IsPhysicalPrimary(mothIdx))
+            //{
 //     std::cout << "!!!!!!!!!!!!!!!!! Primary K0S !!!!!!!!!!!!!!!!!!" << std::endl;
-       //return mothIdx;
-      //      }
-      if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
-      else return GetPrimMother(mothIdx, stack);
-    }
-    else 
-    {
-      std::cout << "WAT!" << std::endl;
-      return -1;
+            //return mothIdx;
+            //      }
+            if(stack->IsPhysicalPrimary(mothIdx)) return mothIdx;
+            else return GetPrimMother(mothIdx, stack);
+        }
+        else
+        {
+            return -1;
+        }
     }
-  }
-  return -1;
+    return -1;
 }
 
 Int_t AliAnalysisEtMonteCarlo::GetK0InFamily(Int_t partIdx, AliStack* stack)
 {
- if(partIdx >= 0) 
-  {
-    if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
-    Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
-    if(mothIdx < 0) return -1;
-    TParticle *mother = stack->Particle(mothIdx);
-    if(mother)
-    {
-      if(mother->GetPdgCode() == fgK0SCode) 
-      {        
-       return mothIdx;
-      }
-      return GetK0InFamily(mothIdx, stack);
-    }
-    else 
+    if(partIdx >= 0)
     {
-      std::cout << "WAT WAT!" << std::endl;
-      return -1;
-    }
-  }
-  return -1;
-}
-
-
-
-
-
-
-
-
-
-Bool_t AliAnalysisEtMonteCarlo::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
-{
-    Float_t gPos[3];
-    cluster.GetPosition(gPos);
-    Int_t relId[4];
-    TVector3 glVec(gPos);
-    fGeoUtils->GlobalPos2RelId(glVec, relId);
-
-    TVector3 locVec;
-    fGeoUtils->Global2Local(locVec, glVec, relId[0]);
-//    std::cout << fGeoUtils << std::endl;
-    //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
-    //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
-    for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
-    {
-        for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
+        if(stack->Particle(partIdx)->GetPdgCode() == fgK0SCode) return partIdx;
+        Int_t mothIdx = stack->Particle(partIdx)->GetMother(0);
+        if(mothIdx < 0) return -1;
+        TParticle *mother = stack->Particle(mothIdx);
+        if(mother)
         {
-            if (relId[0] == 3)
-            {
-                if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
-                {
-                   Int_t tmpRel[4];
-                   tmpRel[0] = 3;
-                   tmpRel[1] = 0;
-                   tmpRel[2] = x+1;
-                   tmpRel[3] = z+1;
-                   
-                   Float_t tmpX;
-                   Float_t tmpZ;
-                   fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-                    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-                   
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                   std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 2)
+            if(mother->GetPdgCode() == fgK0SCode)
             {
-                if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
-                {
-                   Int_t tmpRel[4];
-                   tmpRel[0] = 2;
-                   tmpRel[1] = 0;
-                   tmpRel[2] = x+1;
-                   tmpRel[3] = z+1;
-                   
-                   Float_t tmpX;
-                   Float_t tmpZ;
-                   fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-
-//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                   //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                   std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 1)
-            {
-                if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
-                {
-                   Int_t tmpRel[4];
-                   tmpRel[0] = 1;
-                   tmpRel[1] = 0;
-                   tmpRel[2] = x+1;
-                   tmpRel[3] = z+1;
-                   
-                   Float_t tmpX;
-                   Float_t tmpZ;
-                   fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-
-//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                   //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                     std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
+                return mothIdx;
             }
-
+            return GetK0InFamily(mothIdx, stack);
+        }
+        else
+        {
+            return -1;
         }
     }
-
-    return kFALSE;
-
-
+    return -1;
 }
 
-
-
-
-
-
index 6880868..a861fcb 100644 (file)
@@ -49,16 +49,42 @@ protected:
     Int_t PrintFamilyTree(Int_t partIdx, AliStack *stack);
     Int_t PrintMothers(Int_t partIdx, AliStack *stack, Int_t gen);
     
-    
-    virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
-
-
 protected:
 
     Double_t fImpactParameter; // b(fm), for Hijing; 0 otherwise
     Int_t fNcoll; // Ncoll, for Hijing; 1 otherwise
     Int_t fNpart; // Ncoll, for Hijing; 2 otherwise
+    
+    TTree *fPrimaryTree; // Tree holding info on primaries
+
+    Double_t fTotEtWithSecondaryRemoved;
+    Double_t fTotEtSecondaryFromEmEtPrimary;
+    Double_t fTotEtSecondary;
+    
+    Int_t fPrimaryCode;
+    Int_t fPrimaryCharge;
+
+    Double_t fPrimaryE;
+    Double_t fPrimaryEt;
 
+    Double_t fPrimaryPx;
+    Double_t fPrimaryPy;
+    Double_t fPrimaryPz;
+    
+    Double_t fPrimaryVx;
+    Double_t fPrimaryVy;
+    Double_t fPrimaryVz;
+    
+    Bool_t fPrimaryAccepted;
+    Int_t fDepositedCode;
+    Double_t fDepositedEt;
+    Int_t fDepositedCharge;
+
+    Double_t fDepositedVx;
+    Double_t fDepositedVy;
+    Double_t fDepositedVz;
+    
+    
     TH3F *fHistDecayVertexNonRemovedCharged; // Decay vertex for non-removed charged particles
     TH3F *fHistDecayVertexRemovedCharged; // Decay vertex for non-removed charged particles
     TH3F *fHistDecayVertexNonRemovedNeutral; // Decay vertex for non-removed charged particles
@@ -243,12 +269,8 @@ protected:
 
     Int_t fNClusters; // Number of clusters in event
 
-    AliPHOSGeometry *fGeoUtils;
-
-    TH2I *fBadMapM2; // Bad map
-    TH2I *fBadMapM3; // Bad map
-    TH2I *fBadMapM4; // Bad map
-
+    Double_t fTotNeutralEtAfterMinEnergyCut;
+    
 private:
 
     //Declare it private to avoid compilation warning
index 660de03..c834529 100644 (file)
@@ -29,14 +29,6 @@ AliAnalysisEtMonteCarloEmcal::~AliAnalysisEtMonteCarloEmcal()
 void AliAnalysisEtMonteCarloEmcal::Init()
 { // Init
   AliAnalysisEtMonteCarlo::Init();
-   fClusterType = fCuts->GetEmcalClusterType();
   fDetectorRadius = fCuts->GetGeometryEmcalDetectorRadius();
-  fEtaCutAcc = fCuts->GetGeometryEmcalEtaAccCut();
-  fPhiCutAccMax = fCuts->GetGeometryEmcalPhiAccMaxCut() * TMath::Pi()/180.;
-  fPhiCutAccMin = fCuts->GetGeometryEmcalPhiAccMinCut() * TMath::Pi()/180.;
-  fClusterEnergyCut = fCuts->GetReconstructedEmcalClusterEnergyCut();
   fSingleCellEnergyCut = fCuts->GetReconstructedEmcalSingleCellEnergyCut();
-  
-  fDetector = fCuts->GetDetectorEmcal();
-
 }
index 91c29f2..89c1eec 100644 (file)
@@ -25,6 +25,7 @@ AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCar
 ,fBadMapM2(0)
 ,fBadMapM3(0)
 ,fBadMapM4(0)
+,fGeoUtils(0)
 {
    fHistogramNameSuffix = TString("PhosMC");
 }
@@ -39,21 +40,9 @@ void AliAnalysisEtMonteCarloPhos::Init()
   AliAnalysisEtMonteCarlo::Init();
   fSelector = new AliAnalysisEtSelectorPhos(fCuts);
     
-    fSelector->Init(137366);
-   fClusterType = fCuts->GetPhosClusterType();    
   fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
-  fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
-  fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
-  fPhiCutAccMin = fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180.;
-  fClusterEnergyCut = fCuts->GetReconstructedPhosClusterEnergyCut();
   fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
-  
-  fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
-  fTrackDxCut = fCuts->GetPhosTrackDxCut();
-  fTrackDzCut = fCuts->GetPhosTrackDzCut();
-    
-  fDetector = fCuts->GetDetectorPhos();
-  
+
  // ifstream f("badchannels.txt", ios::in);
   TFile *f = TFile::Open("badchannels.root", "READ");
   
@@ -61,7 +50,6 @@ void AliAnalysisEtMonteCarloPhos::Init()
    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 // 
-//  fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
    fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
   
 }
index 09881ec..1ac18eb 100644 (file)
@@ -1,7 +1,8 @@
 #ifndef ALIANALYSISETMONTECARLOPHOS_H
 #define ALIANALYSISETMONTECARLOPHOS_H
 
-class AliPHOSGeoUtils;//_________________________________________________________________________
+class AliPHOSGeoUtils;
+//_________________________________________________________________________
 //  Utility Class for transverse energy studies
 //  Base class for MC analysis, for PHOS
 //  - MC output
@@ -10,6 +11,7 @@ class AliPHOSGeoUtils;//________________________________________________________
 //_________________________________________________________________________
 
 #include "AliAnalysisEtMonteCarlo.h"
+
 class TH2I;
 
 class AliAnalysisEtMonteCarloPhos : public AliAnalysisEtMonteCarlo
@@ -29,6 +31,8 @@ protected:
     TH2I *fBadMapM2; // Bad map
     TH2I *fBadMapM3; // Bad map
     TH2I *fBadMapM4; // Bad map
+
+    AliPHOSGeoUtils *fGeoUtils; // Geo utils
     
     // Prohibited
     AliAnalysisEtMonteCarloPhos & operator = (const AliAnalysisEtMonteCarloPhos&) ;//cpy assignment
index 3d73863..5272533 100644 (file)
@@ -41,32 +41,28 @@ ClassImp(AliAnalysisEtReconstructed);
 
 
 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
-    AliAnalysisEt()
-    ,fCorrections(0)
-    ,fPidCut(0)
-    ,fHistChargedPionEnergyDeposit(0)
-    ,fHistProtonEnergyDeposit(0)
-    ,fHistAntiProtonEnergyDeposit(0)
-    ,fHistChargedKaonEnergyDeposit(0)
-    ,fHistMuonEnergyDeposit(0)
-    ,fHistRemovedEnergy(0)
-    ,fGeomCorrection(1.0)
-    ,fEMinCorrection(1.0/0.89)
-    ,fRecEffCorrection(1.0)
-    ,fGeoUtils(0)
-    ,fBadMapM2(0)
-    ,fBadMapM3(0)
-    ,fBadMapM4(0)
-    ,fClusterPosition(0)
-    ,fHistChargedEnergyRemoved(0)
-    ,fHistNeutralEnergyRemoved(0)
-    ,fHistGammaEnergyAdded(0)
+        AliAnalysisEt()
+        ,fCorrections(0)
+        ,fPidCut(0)
+        ,fHistChargedPionEnergyDeposit(0)
+        ,fHistProtonEnergyDeposit(0)
+        ,fHistAntiProtonEnergyDeposit(0)
+        ,fHistChargedKaonEnergyDeposit(0)
+        ,fHistMuonEnergyDeposit(0)
+        ,fHistRemovedEnergy(0)
+        ,fGeomCorrection(1.0)
+        ,fEMinCorrection(1.0/0.89)
+       ,fRecEffCorrection(1.0)
+       ,fClusterPosition(0)
+       ,fHistChargedEnergyRemoved(0)
+       ,fHistNeutralEnergyRemoved(0)
+       ,fHistGammaEnergyAdded(0)
 {
 
 }
 
 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
-{   //destructor
+{//destructor
     delete fCorrections;
     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
@@ -80,9 +76,10 @@ AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
 
 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 {
+
+    //AliAnalysisEt::AnalyseEvent(ev);
     // analyse ESD event
     ResetEventValues();
-
     if (!ev) {
         AliFatal("ERROR: Event does not exist");
         return 0;
@@ -94,17 +91,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         return 0;
     }
     fSelector->SetEvent(event);
-    if (!fMatrixInitialized)
-    {
-        for (Int_t mod=0; mod<5; mod++) {
-            if (!event->GetPHOSMatrix(mod)) continue;
-            fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
-//         std::cout << event->GetPHOSMatrix(mod) << std::endl;
-            Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
-        }
-        fMatrixInitialized = kTRUE;
-    }
-
+    
     Int_t cent = -1;
     if (fCentrality)
     {
@@ -112,8 +99,6 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         fCentClass = fCentrality->GetCentralityClass10("V0M");
     }
 
-    //Double_t protonMass = fgProtonMass;
-
     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
     {
         AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
@@ -123,20 +108,20 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
             continue;
         }
         int x = 0;
-        fCutFlow->Fill(x++);
-        if(cluster->IsEMCAL()) continue;
-        fCutFlow->Fill(x++);
-        if(!fSelector->CutMinEnergy(*cluster)) continue;
-        fCutFlow->Fill(x++);
+       fCutFlow->Fill(x++);
+       if(cluster->IsEMCAL()) continue;
+       fCutFlow->Fill(x++);
+       if(!fSelector->CutMinEnergy(*cluster)) continue;
+       fCutFlow->Fill(x++);
         if (!fSelector->CutDistanceToBadChannel(*cluster)) continue;
-        fCutFlow->Fill(x++);
+       fCutFlow->Fill(x++);
 
         Float_t pos[3];
 
         cluster->GetPosition(pos);
         TVector3 cp(pos);
-
-        Double_t distance = cluster->GetEmcCpvDistance();
+       
+       Double_t distance = cluster->GetEmcCpvDistance();
         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
         if ( cluster->IsEMCAL() ) {
             distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
@@ -144,13 +129,12 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 
         Bool_t matched = false;
 
-        Double_t r = 0;
-
-        matched = !fSelector->CutTrackMatching(*cluster, r);
+       
+       matched = !fSelector->CutTrackMatching(*cluster);
 
         if (matched)
         {
-
+         
             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
             {
                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
@@ -173,10 +157,10 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
                                 maxpid = p;
                             }
                         }
-                        if (fCuts->GetHistMakeTreeDeposit() && fTreeDeposit)
+                        if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
                         {
                             fEnergyDeposited = cluster->E();
-                            fEnergyTPC = track->E();
+                            fMomentumTPC = track->P();
                             fCharge = track->Charge();
                             fParticlePid = maxpid;
                             fPidProb = maxpidweight;
@@ -186,39 +170,35 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
                             }
                             else {
                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
-                                fTreeDeposit->Fill();
+                                fDepositTree->Fill();
                             }
                         }
 
                         if (maxpidweight > fPidCut)
                         {
-                            Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+                            //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
 
-                            Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
+                            //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
 
-                            Float_t et = cluster->E() * TMath::Sin(theta);
+                            //Float_t et = cluster->E() * TMath::Sin(theta);
                             if (maxpid == AliPID::kProton)
                             {
 
                                 if (track->Charge() == 1)
                                 {
-                                    fBaryonEt += et;
                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
                                 }
                                 else if (track->Charge() == -1)
                                 {
-                                    fAntiBaryonEt += et;
                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
                                 }
                             }
                             else if (maxpid == AliPID::kPion)
                             {
-                                fMesonEt += et;
                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
                             }
                             else if (maxpid == AliPID::kKaon)
                             {
-                                fMesonEt += et;
                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
                             }
                             else if (maxpid == AliPID::kMuon)
@@ -233,57 +213,44 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         } // distance
         else
         {
-            fCutFlow->Fill(x++);
-            //std::cout << x++ << std::endl;
-            fSparseClusters[0] = AliPID::kPhoton;
-            fSparseClusters[1] = 0;
+         fCutFlow->Fill(x++);
+         //std::cout << x++ << std::endl;
 
             //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
             //if (cluster->E() < fClusterEnergyCut) continue;
             cluster->GetPosition(pos);
+           
+           TVector3 p2(pos);
+           
+           fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
 
-            TVector3 p2(pos);
-
-            fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
-
-            fTotNeutralEt += CalculateTransverseEnergy(cluster);
+           fTotNeutralEt += CalculateTransverseEnergy(cluster);
             fNeutralMultiplicity++;
         }
         fMultiplicity++;
     }
-
+    
     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
-
+    
     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 
-    Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity);
+    Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
     fHistRemovedEnergy->Fill(removedEnergy);
-
+    
     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
-    fTotNeutralEtAcc = fTotNeutralEt;
     fTotEt = fTotChargedEt + fTotNeutralEt;
-    fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
-    if(fMakeSparse) {
-        fSparseEt[0] = fTotEt;
-        fSparseEt[1] = fTotNeutralEt;
-        fSparseEt[2] = fTotChargedEtAcc;
-        fSparseEt[3] = fMultiplicity;
-        fSparseEt[4] = fNeutralMultiplicity;
-        fSparseEt[5] = fChargedMultiplicity;
-        fSparseEt[6] = cent;
-    }
-    // Fill the histograms...
+// Fill the histograms...0
     FillHistograms();
-
+   // std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
     return 0;
 }
 
 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
-{   // check vertex
+{ // check vertex
 
     Float_t bxy = 999.;
     Float_t bz = 999.;
@@ -309,25 +276,17 @@ bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
 }
 
 void AliAnalysisEtReconstructed::Init()
-{   // Init
+{ // Init
     AliAnalysisEt::Init();
     fPidCut = fCuts->GetReconstructedPidCut();
     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
     if (!fCorrections) {
         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
     }
-    //fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
-    fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
-    // ifstream f("badchannels.txt", ios::in);
-    TFile *f = TFile::Open("badchannels.root", "READ");
-
-    fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
-    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
-    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
 }
 
 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
-{   // propagate track to detector radius
+{ // propagate track to detector radius
 
     if (!track) {
         cout<<"Warning: track empty"<<endl;
@@ -343,14 +302,11 @@ bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Doubl
     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
 
     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
-    return prop &&
-           TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
-           esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
-           esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
+    return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
 }
 
 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
-{   // add some extra histograms to the ones from base class
+{ // add some extra histograms to the ones from base class
     AliAnalysisEt::FillOutputList(list);
 
     list->Add(fHistChargedPionEnergyDeposit);
@@ -361,14 +317,14 @@ void AliAnalysisEtReconstructed::FillOutputList(TList* list)
 
     list->Add(fHistRemovedEnergy);
     list->Add(fClusterPosition);
-
+    
     list->Add(fHistChargedEnergyRemoved);
     list->Add(fHistNeutralEnergyRemoved);
     list->Add(fHistGammaEnergyAdded);
 }
 
 void AliAnalysisEtReconstructed::CreateHistograms()
-{   // add some extra histograms to the ones from base class
+{ // add some extra histograms to the ones from base class
     AliAnalysisEt::CreateHistograms();
 
     Int_t nbinsEt = 1000;
@@ -435,7 +391,7 @@ Double_t
 AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
         Int_t *trkMatchId,
         const AliESDEvent *event)
-{   // calculate distance between cluster and closest track
+{ // calculate distance between cluster and closest track
 
     Double_t trkPos[3] = {0,0,0};
 
@@ -477,170 +433,3 @@ AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
     *trkMatchId = bestTrkMatchId;
     return distance;
 }
-
-
-Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
-{
-
-    Float_t gPos[3];
-    cluster.GetPosition(gPos);
-    Int_t relId[4];
-    TVector3 glVec(gPos);
-    fGeoUtils->GlobalPos2RelId(glVec, relId);
-
-    TVector3 locVec;
-    fGeoUtils->Global2Local(locVec, glVec, relId[0]);
-//    std::cout << fGeoUtils << std::endl;
-    //std::cout << relId[0] << " " << cluster.IsPHOS() <<  std::endl;
-    //std::cout << locVec[0] << " " << " " << locVec[1] << " " << locVec[2] << std::endl;
-    for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
-    {
-        for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
-        {
-            if (relId[0] == 3)
-            {
-                if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
-                {
-                    Int_t tmpRel[4];
-                    tmpRel[0] = 3;
-                    tmpRel[1] = 0;
-                    tmpRel[2] = x+1;
-                    tmpRel[3] = z+1;
-
-                    Float_t tmpX;
-                    Float_t tmpZ;
-                    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-                    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                   std::cout << "Module 2, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 2)
-            {
-                if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
-                {
-                    Int_t tmpRel[4];
-                    tmpRel[0] = 2;
-                    tmpRel[1] = 0;
-                    tmpRel[2] = x+1;
-                    tmpRel[3] = z+1;
-
-                    Float_t tmpX;
-                    Float_t tmpZ;
-                    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-
-//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                   std::cout << "Module 3, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 1)
-            {
-                if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
-                {
-                    Int_t tmpRel[4];
-                    tmpRel[0] = 1;
-                    tmpRel[1] = 0;
-                    tmpRel[2] = x+1;
-                    tmpRel[3] = z+1;
-
-                    Float_t tmpX;
-                    Float_t tmpZ;
-                    fGeoUtils->RelPosInModule(tmpRel, tmpX, tmpZ);
-
-                    Float_t distance = TMath::Sqrt((tmpX-locVec[0])*(tmpX-locVec[0]) + (tmpZ - locVec[2])*(tmpZ-locVec[2]));
-
-//                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                    //Float_t distance = TMath::Sqrt((x-relId[3])*(x-relId[3]) + (z - relId[2])*(z-relId[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-//                     std::cout << "Module 4, position: " << locVec[0] << ", " << locVec[2] << ", distance to bad channel: " << distance << ", number of cells: " << cluster.GetNCells() <<  std::endl;
-                        return kTRUE;
-                    }
-                }
-            }
-
-        }
-    }
-
-    return kFALSE;
-
-
-}
-
-
-
-
-/*
-Bool_t AliAnalysisEtReconstructed::TooCloseToBadChannel(const AliESDCaloCluster &cluster) const
-{
-
-    Float_t gPos[3];
-
-    cluster.GetPosition(gPos);
-    Int_t relId[4];
-    TVector3 glVec(gPos);
-    fGeoUtils->GlobalPos2RelId(glVec, relId);
-    TVector3 locVec;
-    fGeoUtils->Global2Local(locVec, glVec, relId[0]);
-
-    std::vector<Int_t>::const_iterator badIt;
-
-    for (Int_t x = 0; x < fBadMapM2->GetNbinsX(); x++)
-    {
-        for (Int_t z = 0; z < fBadMapM2->GetNbinsY(); z++)
-        {
-            if (relId[0] == 3)
-            {
-                if (fBadMapM2->GetBinContent(x+1, z+1) != 0)
-                {
-
-                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 2)
-            {
-                if (fBadMapM3->GetBinContent(x+1, z+1) != 0)
-                {
-
-                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-                        return kTRUE;
-                    }
-                }
-            }
-            if (relId[0] == 1)
-            {
-                if (fBadMapM4->GetBinContent(x+1, z+1) != 0)
-                {
-
-                    Float_t distance = TMath::Sqrt((x-locVec[0])*(x-locVec[0]) + (z - locVec[2])*(z-locVec[2]));
-                    if (distance < fCuts->GetPhosBadDistanceCut())
-                    {
-                        return kTRUE;
-                    }
-                }
-            }
-        }
-    }
-
-    return kFALSE;
-}
-*/
index 718487d..3d4070e 100644 (file)
@@ -3,8 +3,8 @@
 
 class TH2F;
 class TH2D;
-class AliPHOSGeometry;
-class AliPHOSGeoUtils;//_________________________________________________________________________
+
+//_________________________________________________________________________
 //  Utility Class for transverse energy studies
 //  Base class for ESD analysis
 //  - reconstruction output
@@ -46,9 +46,6 @@ protected:
 
     bool CheckGoodVertex(AliVParticle *track);
     virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
-    virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
-
-    //Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
 
     AliAnalysisHadEtCorrections *fCorrections;//corrections needed for hadronic et
 
@@ -68,13 +65,6 @@ protected:
     Double_t fRecEffCorrection; // Eff correction
     Double_t CalcTrackClusterDistance(const Float_t pos[3],Int_t *trkMatchId, const AliESDEvent *event);
     
-//    AliPHOSGeoUtils *fGeoUtils;
-    AliPHOSGeometry *fGeoUtils;
-    TH2I *fBadMapM2; // Bad map
-    TH2I *fBadMapM3; // Bad map
-    TH2I *fBadMapM4; // Bad map
-
-    
     TH2D *fClusterPosition; // Position of clusters
     
     TH2D *fHistChargedEnergyRemoved;
index 4814a70..98c6ce8 100644 (file)
@@ -31,17 +31,8 @@ void AliAnalysisEtReconstructedEmcal::Init()
   AliAnalysisEtReconstructed::Init();
     
   fDetectorRadius = fCuts->GetGeometryEmcalDetectorRadius();
-  fEtaCutAcc = fCuts->GetGeometryEmcalEtaAccCut();
-  fPhiCutAccMax = fCuts->GetGeometryEmcalPhiAccMaxCut() * TMath::Pi()/180.;
-  fPhiCutAccMin = fCuts->GetGeometryEmcalPhiAccMinCut() * TMath::Pi()/180.;
-  fClusterEnergyCut = fCuts->GetReconstructedEmcalClusterEnergyCut();
   fSingleCellEnergyCut = fCuts->GetReconstructedEmcalSingleCellEnergyCut();
 
-  fClusterType = fCuts->GetEmcalClusterType();
-  fTrackDistanceCut = fCuts->GetEmcalTrackDistanceCut();
-  
-  fDetector = fCuts->GetDetectorEmcal();
-        
 }
 
 bool AliAnalysisEtReconstructedEmcal::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
index be2bb69..0abe592 100644 (file)
@@ -43,17 +43,32 @@ AliAnalysisEtReconstructedPhos::AliAnalysisEtReconstructedPhos() :
 {
     fHistogramNameSuffix = TString("PhosRec");
 
-    fChargedContributionCorrectionParameters[0] = -0.017;
-    fChargedContributionCorrectionParameters[1] = 0.065;
-
-    fNeutralContributionCorrectionParameters[0] = -0.002;
-    fNeutralContributionCorrectionParameters[1] = 0.017;
-    fNeutralContributionCorrectionParameters[2] = -3.6e-5;
-
-    fRemovedGammaContributionCorrectionParameters[0] = 0.001;
-    fRemovedGammaContributionCorrectionParameters[1] = 0.37e-5;
-    fRemovedGammaContributionCorrectionParameters[2] = 0.0003;
-
+//     fChargedContributionCorrectionParameters[0] = -0.017;
+//     fChargedContributionCorrectionParameters[1] = 0.065;
+// 
+//     fNeutralContributionCorrectionParameters[0] = -0.002;
+//     fNeutralContributionCorrectionParameters[1] = 0.017;
+//     fNeutralContributionCorrectionParameters[2] = -3.6e-5;
+// 
+//     fRemovedGammaContributionCorrectionParameters[0] = 0.001;
+//     fRemovedGammaContributionCorrectionParameters[1] = 0.37e-5;
+//     fRemovedGammaContributionCorrectionParameters[2] = 0.0003;
+
+    fChargedContributionCorrectionParameters[0] = -0.0231864;
+    fChargedContributionCorrectionParameters[1] = 0.112661;
+    fChargedContributionCorrectionParameters[2] = 9.2836e-05;
+
+    fNeutralContributionCorrectionParameters[0] = -9.18740e-03;
+    fNeutralContributionCorrectionParameters[1] = 3.58584e-02;
+    fNeutralContributionCorrectionParameters[2] = 9.48208e-04;
+
+    fRemovedGammaContributionCorrectionParameters[0] = -2.74323e-03;
+    fRemovedGammaContributionCorrectionParameters[1] = 2.71104e-03;
+    fRemovedGammaContributionCorrectionParameters[2] = 3.52758e-04;
+
+    fSecondaryContributionCorrectionParameters[0] = 0.075;
+    fSecondaryContributionCorrectionParameters[1] = 0.150187;
+    fSecondaryContributionCorrectionParameters[2] = 0.00329361;
 }
 
 AliAnalysisEtReconstructedPhos::~AliAnalysisEtReconstructedPhos()
@@ -66,22 +81,9 @@ void AliAnalysisEtReconstructedPhos::Init()
 
     fSelector = new AliAnalysisEtSelectorPhos(fCuts);
     
-    fSelector->Init(137366);
-    
     fDetectorRadius = fCuts->GetGeometryPhosDetectorRadius();
-    fEtaCutAcc = fCuts->GetGeometryPhosEtaAccCut();
-    fPhiCutAccMax = fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180.;
-    fPhiCutAccMin = fCuts->GetGeometryPhosPhiAccMinCut() * TMath::Pi()/180.;
-    fClusterEnergyCut = fCuts->GetReconstructedPhosClusterEnergyCut();
     fSingleCellEnergyCut = fCuts->GetReconstructedPhosSingleCellEnergyCut();
 
-    fClusterType = fCuts->GetPhosClusterType();
-    fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
-    fTrackDxCut = fCuts->GetPhosTrackDxCut();
-    fTrackDzCut = fCuts->GetPhosTrackDzCut();
-
-    fDetector = fCuts->GetDetectorPhos();
-  
 
 }
 
@@ -94,7 +96,7 @@ Double_t AliAnalysisEtReconstructedPhos::GetChargedContribution(Int_t clusterMul
 { // Charged contrib
     if (clusterMult > 0)
     {
-        Double_t nPart = fChargedContributionCorrectionParameters[0] + fChargedContributionCorrectionParameters[1]*clusterMult;
+        Double_t nPart = fChargedContributionCorrectionParameters[0] + fChargedContributionCorrectionParameters[1]*clusterMult + fChargedContributionCorrectionParameters[2]*clusterMult*clusterMult;
 
         Double_t contr = nPart*kMEANCHARGED;
 
@@ -130,6 +132,15 @@ Double_t AliAnalysisEtReconstructedPhos::GetGammaContribution(Int_t clusterMult)
     return 0;
 }
 
+Double_t AliAnalysisEtReconstructedPhos::GetSecondaryContribution(Int_t clusterMultiplicity)
+{
+  if(clusterMultiplicity > 0)
+  {
+    return fSecondaryContributionCorrectionParameters[0] + fSecondaryContributionCorrectionParameters[1]*clusterMultiplicity + fSecondaryContributionCorrectionParameters[2]*clusterMultiplicity*clusterMultiplicity;
+  }
+  
+  return 0;
+}
 
 
 
index 0365e12..a7f32ee 100644 (file)
@@ -27,11 +27,14 @@ public:
     virtual Double_t GetNeutralContribution(Int_t clusterMultiplicity);
 
     virtual Double_t GetGammaContribution(Int_t clusterMultiplicity);
+    
+    virtual Double_t GetSecondaryContribution(Int_t clusterMultiplicity);
 
-    void SetChargedContributionParameters(Double_t par[2])
+    void SetChargedContributionParameters(Double_t par[3])
     {
         fChargedContributionCorrectionParameters[0] = par[0];
         fChargedContributionCorrectionParameters[1] = par[1];
+       fChargedContributionCorrectionParameters[1] = par[2];
     }
     void SetNeutralContributionParameters(Double_t par[3])
     {
@@ -45,6 +48,13 @@ public:
         fRemovedGammaContributionCorrectionParameters[1] = par[1];
        fRemovedGammaContributionCorrectionParameters[2] = par[2];
     }
+    void SetSecondaryContributionParameters(Double_t par[3])
+    {
+        fSecondaryContributionCorrectionParameters[0] = par[0];
+        fSecondaryContributionCorrectionParameters[1] = par[1];
+       fSecondaryContributionCorrectionParameters[2] = par[2];
+    }
+    
 
 protected:
 
@@ -52,10 +62,12 @@ protected:
     
 private:
 
-    Double_t fChargedContributionCorrectionParameters[2]; // Parametrization of the charged contribution as function of cluster multiplicity
+    Double_t fChargedContributionCorrectionParameters[3]; // Parametrization of the charged contribution as function of cluster multiplicity
     Double_t fNeutralContributionCorrectionParameters[3]; // Parametrization of the neutral contribution as function of cluster multiplicity
     Double_t fRemovedGammaContributionCorrectionParameters[3]; // Parametrization of the negative contribution from removed gammas as function of cluster multiplicity
 
+    Double_t fSecondaryContributionCorrectionParameters[3]; // Parametrization of the positive contribution of secondary particles
+
 
     ClassDef(AliAnalysisEtReconstructedPhos, 1);
 };
index d31abaa..91e1545 100644 (file)
@@ -1,19 +1,62 @@
 #include "AliAnalysisEtSelector.h"
 #include "AliAnalysisEtCuts.h"
-#include <AliESDCaloCluster.h>
+#include "AliESDCaloCluster.h"
+#include "AliStack.h"
+#include "TParticle.h"
+#include "TParticlePDG.h"
+#include "AliAnalysisEtCommon.h"
+#include <iostream>
 
-ClassImp(AliAnalysisEtSelector)
+ClassImp(AliAnalysisEtSelector);
 
-AliAnalysisEtSelector::AliAnalysisEtSelector(AliAnalysisEtCuts *cuts) :
-fEvent(0)
+AliAnalysisEtSelector::AliAnalysisEtSelector(AliAnalysisEtCuts *cuts) : AliAnalysisEtCommon()
+,fEvent(0)
 ,fCuts(cuts)
+,fClusterArray(0)
 ,fRunNumber(0)
 {
-
 }
 
 AliAnalysisEtSelector::~AliAnalysisEtSelector()
 {
+  if(fClusterArray)
+  {
+    delete fClusterArray;
+  }
+}
+
+Bool_t AliAnalysisEtSelector::CutNeutralMcParticle(Int_t pIdx, AliStack& s, const TParticlePDG& pdg) const
+{
+  return s.IsPhysicalPrimary(pIdx) &&(TMath::Abs(TMath::Abs(pdg.Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3);
+}
 
+Bool_t AliAnalysisEtSelector::IsEmEtParticle(const Int_t pdgCode) const
+{
+  return pdgCode == fgGammaCode || pdgCode == fgPi0Code || pdgCode == fgEtaCode || pdgCode == fgEPlusCode || pdgCode == fgEMinusCode;
 }
 
+
+Bool_t AliAnalysisEtSelector::PrimaryIsEmEtParticle(const Int_t pIdx, AliStack& stack) const
+{
+  return IsEmEtParticle(stack.Particle(GetPrimary(pIdx, stack))->GetPdgCode());
+}
+Int_t AliAnalysisEtSelector::GetPrimary(const Int_t partIdx, AliStack& stack) const
+{
+  if(partIdx >= 0) 
+  {
+    Int_t mothIdx = stack.Particle(partIdx)->GetMother(0);
+    if(mothIdx < 0) return -1;
+    TParticle *mother = stack.Particle(mothIdx);
+    if(mother)
+    {
+      if(stack.IsPhysicalPrimary(mothIdx)) return mothIdx;
+      else return GetPrimary(mothIdx, stack);
+    }
+    else 
+    {
+      std::cout << "WAT!" << std::endl;
+      return -1;
+    }
+  }
+  return -1;
+}
index b44660c..dc59d23 100644 (file)
@@ -1,12 +1,18 @@
 #ifndef ALIANALYSISETSELECTOR_H
 #define ALIANALYSISETSELECTOR_H
 #include <Rtypes.h>
+#include "AliAnalysisEtCommon.h"
+#include "AliESDEvent.h"
 
-class AliVEvent;
 class AliESDCaloCluster;
+class AliVTrack;
 class TRefArray;
 class AliAnalysisEtCuts;
-class AliAnalysisEtSelector
+class TParticle;
+class TParticlePDG;
+class AliStack;
+
+class AliAnalysisEtSelector : public AliAnalysisEtCommon
 {
 
 public:
@@ -18,31 +24,60 @@ public:
     virtual ~AliAnalysisEtSelector();
     
     // Set the current event
-    void SetEvent(const AliVEvent *event) { fEvent = event; }
+    virtual void SetEvent(const AliESDEvent *event) = 0;// { fEvent = event; }
+    
+    // Init
+    virtual void Init() {} 
     
-    // Init 
-    virtual Int_t Init(Int_t runNumber) { fRunNumber = runNumber; return 0; }
+    // Init with event
+    virtual Int_t Init(const AliESDEvent *event) { fRunNumber = event->GetRunNumber(); return 0; }
     
     // Return CaloClusters for the detector
-    virtual TRefArray* GetClusters() = 0;
+    virtual TRefArray* GetClusters() { return 0; }
+    
+    // Return true if cluster has energy > cut
+    virtual Bool_t CutMinEnergy(const AliESDCaloCluster &/*cluster*/) const { return true; }
     
     // Return true if cluster has energy > cut
-    virtual Bool_t CutMinEnergy(const AliESDCaloCluster & /*cluster*/) const { return true; }
+    virtual Bool_t CutMinEnergy(const TParticle &/*part*/) const { return true; }
     
     // Cut on distance to bad channel
-    virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster & /*cluster*/) const { return true; }
+    virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster &/*cluster*/) const { return true; }
     
     // Cut on track matching
-    virtual Bool_t CutTrackMatching(const AliESDCaloCluster& /*cluster*/, Double_t &/*r*/) const { return true; }
-  
+    virtual Bool_t CutTrackMatching(const AliESDCaloCluster &/*cluster*/) const { return true; }
+    
+    // Cut on neutral monte carlo particle
+    virtual Bool_t CutNeutralMcParticle(Int_t pIdx, AliStack& s, const TParticlePDG& pdg) const;
+    
+    // Is it an EM E_T particle
+    virtual Bool_t IsEmEtParticle(const Int_t pdgCode) const;
+    
+    // Does the particle come from an EM E_T primary ?
+    virtual Bool_t PrimaryIsEmEtParticle(const Int_t pIdx, AliStack &stack) const;
+
+    // Get the index of primary particle for the particle
+    Int_t GetPrimary(const Int_t partIdx, AliStack &stack) const;
+    
+    // Cut on geometrical acceptance 
+    virtual Bool_t CutGeometricalAcceptance(const TParticle &/*part*/) const { return true; }
+    
+    // Cut on geometrical acceptance 
+    virtual Bool_t CutGeometricalAcceptance(const AliVTrack &/*part*/) const { return true; }
+    
+    
 protected:
   
     const AliVEvent *fEvent; // Pointer to current event
 
     AliAnalysisEtCuts *fCuts; // Pointer to the cuts object
     
+    TRefArray *fClusterArray; // Array of clusters
+    
     Int_t fRunNumber;
     
+
+    
 private:
 
     AliAnalysisEtSelector(); // Prohibited
index 085dc91..3c6fe99 100644 (file)
@@ -1,20 +1,25 @@
 #include "AliAnalysisEtSelectorPhos.h"
 #include "AliAnalysisEtCuts.h"
 #include "AliESDCaloCluster.h"
-#include <AliVEvent.h>
+#include "AliESDEvent.h"
 #include "TRefArray.h"
 #include "AliPHOSGeometry.h"
 #include "TH2I.h"
 #include "TFile.h"
 #include "TMath.h"
+#include "TParticle.h"
 #include "AliLog.h"
 #include <iostream>
+
 ClassImp(AliAnalysisEtSelectorPhos)
+
 AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
 ,fGeoUtils(0)
 ,fBadMapM2(0)
 ,fBadMapM3(0)
 ,fBadMapM4(0)
+,fInitialized(kFALSE)
+,fMatrixInitialized(kFALSE)
 {
   
 }
@@ -25,17 +30,46 @@ AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
 }
 
 TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
-{  
-  return 0;
+{
+  if(!fClusterArray) fClusterArray = new TRefArray;
+  
+  if(fClusterArray)
+  {
+    fEvent->GetPHOSClusters(fClusterArray);
+  }
+  else
+  {
+    Printf("Could not initialize cluster array");
+  }
+  
+  return fClusterArray;
 }
 
-Int_t AliAnalysisEtSelectorPhos::Init(Int_t runNumber)
+Int_t AliAnalysisEtSelectorPhos::Init(const AliESDEvent* event)
 {
-  AliAnalysisEtSelector::Init(runNumber);
   
+  AliAnalysisEtSelector::Init(event);
+  Printf("Initializing selector for run: %d", event->GetRunNumber());
   int res = LoadGeometry();
   if(res) return -1;
-  return LoadBadMaps();  
+  if(LoadBadMaps()) return -1;
+  fInitialized = kTRUE;
+  if (!fMatrixInitialized)
+    {
+       Printf("INITIALIZING MISALIGNMENT MATRICES");
+        for (Int_t mod=0; mod<5; mod++) {
+           
+            if (!event->GetPHOSMatrix(mod))
+           {
+             Printf("Could not find geo matrix for module %d", mod);
+             continue;
+           }
+           fMatrixInitialized = kTRUE;
+            fGeoUtils->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
+            Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
+        }
+    }
+  return 0;
 }
 
 Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
@@ -43,8 +77,19 @@ Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster)
   return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
 }
 
+Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const TParticle& part) const
+{
+    return part.Energy() > fCuts->GetReconstructedPhosClusterEnergyCut();
+}
+
+
 Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const
 {
+  if(!fMatrixInitialized)
+  {
+    Printf("Misalignment matrices are not initialized");
+    return kFALSE;
+  }
     Float_t gPos[3];
     cluster.GetPosition(gPos);
     Int_t relId[4];
@@ -142,9 +187,14 @@ Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(const AliESDCaloCluste
 
 }
 
-Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const
+Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster) const
 {
 
+  if(!fMatrixInitialized)
+  {
+    Printf("Misalignment matrices are not initialized");
+    return kFALSE;
+  }
   
   // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
   
@@ -185,7 +235,7 @@ Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& clus
 
   Double_t rz=(dz-meanZ)/sz ;
   Double_t rx=(dx-meanX)/sx ;
-  r = TMath::Sqrt(rx*rx+rz*rz);
+  Double_t r = TMath::Sqrt(rx*rx+rz*rz);
   if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
   
   return kTRUE;
@@ -209,26 +259,46 @@ TFile *f = TFile::Open("badchannels.root", "READ");
       std::cout << "Could not open badchannels.root" << std::endl;
       return -1;
     }
-    
+
     fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
-    if(fBadMapM2) 
+    if(!fBadMapM2) 
     {
       std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
     }
     fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
-    if(fBadMapM2) 
+    if(!fBadMapM2) 
     {
       std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
     }
     
     fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
-    if(fBadMapM4) 
+    if(!fBadMapM4) 
     {
       std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
     }
-
-    return 0;
     
     
+    return 0;
     
 }
+
+void AliAnalysisEtSelectorPhos::SetEvent(const AliESDEvent* event)
+{
+    //AliAnalysisEtSelector::SetEvent(event);
+    fEvent = event;
+    if(!fInitialized) Init(event);
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const TParticle& part) const
+{
+  return TMath::Abs(part.Eta()) < fCuts->GetGeometryPhosEtaAccCut() 
+         && part.Phi() < fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180.
+         && part.Phi() > fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutGeometricalAcceptance(const AliVTrack& track) const
+{
+  return TMath::Abs(track.Eta()) < fCuts->GetGeometryPhosEtaAccCut() &&
+           track.Phi() > fCuts->GetGeometryPhosPhiAccMaxCut()*TMath::Pi()/180. &&
+           track.Phi() < fCuts->GetGeometryPhosPhiAccMinCut()*TMath::Pi()/180.;
+}
index 67921b3..86664b9 100644 (file)
@@ -4,7 +4,9 @@
 #include "AliAnalysisEtSelector.h"
 
 class TH2I;
+class TParticle;
 class AliPHOSGeometry;
+class AliESDEvent;
 
 class AliAnalysisEtSelectorPhos : public AliAnalysisEtSelector
 {
@@ -16,11 +18,16 @@ public:
     
     virtual TRefArray* GetClusters();
     virtual Bool_t CutMinEnergy(const AliESDCaloCluster& cluster) const;
+    virtual Bool_t CutMinEnergy(const TParticle& part) const;
     virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const;
-    virtual Bool_t CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const;
+    virtual Bool_t CutTrackMatching(const AliESDCaloCluster& cluster) const;
+    virtual Bool_t CutGeometricalAcceptance(const TParticle& part) const;    
+    virtual Bool_t CutGeometricalAcceptance(const AliVTrack& part) const;    
+    virtual void Init() {}
+    virtual Int_t Init(const AliESDEvent *ev);
+
+    virtual void SetEvent(const AliESDEvent* event);
     
-    virtual Int_t Init(int runNumber);
-  
 private:
 
 
@@ -32,6 +39,9 @@ private:
     TH2I *fBadMapM2; // Bad map
     TH2I *fBadMapM3; // Bad map
     TH2I *fBadMapM4; // Bad map
+
+    Bool_t fInitialized; // matrix initialized
+    Bool_t fMatrixInitialized; // matrix initialized
     
     AliAnalysisEtSelectorPhos();
     AliAnalysisEtSelectorPhos(const AliAnalysisEtSelectorPhos& other); // Prohibited
index 68bbd74..195f571 100644 (file)
@@ -208,9 +208,9 @@ void AliAnalysisTaskTotEt::UserExec(Option_t *)
     }
   if(fMCAnalysis)
     {
-      fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEtAcc(), fMCAnalysis->GetTotNeutralEtAcc());
-      if(fMCAnalysis->GetTotNeutralEtAcc()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEtAcc(), cent->GetCentralityClass10("V0M"));
-      if(fMCAnalysis->GetTotNeutralEtAcc()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
+      fHistEtRecvsEtMC->Fill(fRecAnalysis->GetTotNeutralEt(), fMCAnalysis->GetTotNeutralEt());
+      if(fMCAnalysis->GetTotNeutralEt()) fHistEtRecOverEtMC->Fill(fRecAnalysis->GetTotNeutralEt()/fMCAnalysis->GetTotNeutralEt(), cent->GetCentralityClass10("V0M"));
+      if(fMCAnalysis->GetTotNeutralEt()) fHistDiffEtRecEtMCOverEtMC->Fill(fMCAnalysis->GetTotNeutralEt(), (fRecAnalysis->GetTotNeutralEt()-fMCAnalysis->GetTotNeutralEt())/fMCAnalysis->GetTotNeutralEt());
     }
   //}
   // Post output data.