- restructuring the cuts in the PHOS E_T code.
authorodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jun 2012 23:23:37 +0000 (23:23 +0000)
committerodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jun 2012 23:23:37 +0000 (23:23 +0000)
this is only partially complete, and breaks EMCAL compatibility for
this part of the E_T code.

19 files changed:
PWGLF/totEt/AliAnalysisEt.cxx
PWGLF/totEt/AliAnalysisEt.h
PWGLF/totEt/AliAnalysisEtCuts.cxx
PWGLF/totEt/AliAnalysisEtCuts.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 [new file with mode: 0644]
PWGLF/totEt/AliAnalysisEtSelector.h [new file with mode: 0644]
PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx [new file with mode: 0644]
PWGLF/totEt/AliAnalysisEtSelectorPhos.h [new file with mode: 0644]
PWGLF/totEt/AliAnalysisTaskTransverseEnergy.cxx

index 0139f34..1936d2b 100644 (file)
@@ -12,6 +12,7 @@
 #include "TList.h"
 #include "TH1F.h"
 #include "TH2F.h"
+#include "TH1I.h" 
 #include "TTree.h"
 #include <iostream>
 #include "AliAnalysisEtCuts.h"
@@ -21,6 +22,8 @@
 #include "Rtypes.h"
 #include "TString.h"
 #include "AliCentrality.h"
+#include "AliAnalysisEtSelector.h"
+
               //#include "THnSparse.h"
 
 using namespace std;
@@ -95,6 +98,9 @@ AliAnalysisEt::AliAnalysisEt() : AliAnalysisEtCommon()
                               ,fTrackDistanceCut(0)
                               ,fTrackDxCut(0)
                               ,fTrackDzCut(0)
+                              ,fChargedEnergyRemoved(0)
+                              ,fNeutralEnergyRemoved(0)
+                              ,fGammaEnergyAdded(0)
                               ,fHistEt(0)
                               ,fHistChargedEt(0)
                               ,fHistNeutralEt(0)
@@ -128,13 +134,17 @@ AliAnalysisEt::AliAnalysisEt() : AliAnalysisEtCommon()
                               ,fTreeDeposit(0)
                               ,fCentrality(0)
                               ,fDetector(0)
-                              ,fMakeSparse(kFALSE)
+                              ,fMakeSparse(kTRUE)
                               ,fSparseHistTracks(0)
                               ,fSparseHistClusters(0)
                               ,fSparseHistEt(0)
-                              ,fSparseTracks(0)
+                                      ,fSparseTracks(0)
                               ,fSparseClusters(0)
                               ,fSparseEt(0)
+                              ,fClusterType(0)
+                              ,fMatrixInitialized(kFALSE)
+                              ,fCutFlow(0)
+                              ,fSelector(0)
                               
 {}
 
@@ -243,6 +253,8 @@ void AliAnalysisEt::FillOutputList(TList *list)
       list->Add(fSparseHistClusters);
       list->Add(fSparseHistEt);
     }
+    
+    list->Add(fCutFlow);
 
 }
 
@@ -445,6 +457,9 @@ void AliAnalysisEt::CreateHistograms()
       fSparseHistEt->GetAxis(5)->SetTitle("charged_mult");
       fSparseHistEt->GetAxis(6)->SetTitle("cent");
     }
+    
+    histname = "fCutFlow" + fHistogramNameSuffix;
+    fCutFlow = new TH1I(histname, histname, 20, -0.5, 19.5);
 }
 
 TH2F* AliAnalysisEt::CreateEtaEHisto2D(TString name, TString title, TString ztitle)
@@ -591,6 +606,9 @@ void AliAnalysisEt::CreateTrees()
     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");
     
     fTree->Branch("fBaryonEt",&fBaryonEt,"fBaryonEt/D");
     fTree->Branch("fAntiBaryonEt",&fAntiBaryonEt,"fAntiBaryonEt/D");
@@ -666,6 +684,7 @@ void AliAnalysisEt::FillHistograms()
 
 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);
   ResetEventValues();
   return 0;
@@ -708,7 +727,50 @@ Double_t AliAnalysisEt::CalculateTransverseEnergy(AliESDCaloCluster* cluster)
   Float_t pos[3];
   cluster->GetPosition(pos);
   TVector3 cp(pos);
+  Double_t corrEnergy = 0;
+  
+  if(cluster->E() < 1.5)
+  {
+    corrEnergy =cluster->E()/(0.51 + 0.02*cluster->E());
+  }    
+  else
+  {
+    corrEnergy =cluster->E()/(0.51 + 0.02*1.5);
+  }
+  //std::cout << "Original energy: " << cluster->E() << ", corrected energy: " << corrEnergy << std::endl;
   
-  return cluster->E() * TMath::Sin(cp.Theta());
+  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 81efbc5..5ca27b1 100644 (file)
@@ -9,14 +9,20 @@
 #ifndef ALIANALYSISET_H
 #define ALIANALYSISET_H
 
-class AliCentrality;
 #include "AliAnalysisEtCommon.h"
 #include "THnSparse.h"
+#include "AliESDCaloCluster.h"
+#include "AliAnalysisEtCuts.h"
+#include <vector>
+#include "Rtypes.h"
 
+class AliAnalysisEtSelector;
+class AliCentrality;
 class TString;
 class TTree;
 class TH2F;
 class TH1F;
+class TH1I;
 class AliVEvent;
 class TList;
 class TString;
@@ -28,12 +34,12 @@ class AliESDCaloCluster;
 class AliAnalysisEt : public AliAnalysisEtCommon
 {
 public:
-   
+
     AliAnalysisEt();
     virtual ~AliAnalysisEt();
-  
+
 public:
-  
+
     /** Analyse the event! */
 
     virtual Int_t AnalyseEvent(AliVEvent *event);
@@ -44,21 +50,21 @@ public:
     /** Initialise the analysis, must be overloaded. */
     virtual void Init();
 
-    /** 
-    * Creates the histograms, must be overloaded if you want to add your own. 
+    /**
+    * Creates the histograms, must be overloaded if you want to add your own.
     * Uses the fHistogramNameSuffix to create proper histogram names
     */
     virtual void CreateHistograms();
     virtual void CreateTrees();
     TH2F* CreateEtaEHisto2D(TString name, TString title, TString ztitle);
-       TH2F* CreateEtaPtHisto2D(TString name, TString title, TString ztitle);
-       TH2F* CreateEtaEtHisto2D(TString name, TString title, TString ztitle);
-       TH2F* CreateResEHisto2D(TString name, TString title, TString ztitle);
-       TH2F* CreateResPtHisto2D(TString name, TString title, TString ztitle);
+    TH2F* CreateEtaPtHisto2D(TString name, TString title, TString ztitle);
+    TH2F* CreateEtaEtHisto2D(TString name, TString title, TString ztitle);
+    TH2F* CreateResEHisto2D(TString name, TString title, TString ztitle);
+    TH2F* CreateResPtHisto2D(TString name, TString title, TString ztitle);
     THnSparseF* CreateClusterHistoSparse(TString name, TString title);
     THnSparseF* CreateNeutralPartHistoSparse(TString name, TString title);
     THnSparseF* CreateChargedPartHistoSparse(TString name, TString title);
-       
+
     /** Fills the histograms, must be overloaded if you want to add your own */
     virtual void FillHistograms();
 
@@ -66,57 +72,87 @@ public:
     virtual void ResetEventValues();
 
     /** Total Et in the event (without acceptance cuts) */
-    Double_t GetTotEt() const { return fTotEt; }
+    Double_t GetTotEt() const {
+        return fTotEt;
+    }
 
     /** Total Et in the event within the acceptance cuts */
-    Double_t GetTotEtAcc() const { return fTotEtAcc; }
+    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 (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; }
-    
+    Double_t GetTotNeutralEtAcc() const {
+        return fTotNeutralEtAcc;
+    }
+
     /** Total charged Et in the event (without acceptance cuts) */
-    Double_t GetTotChargedEt() const { return fTotChargedEt; }
+    Double_t GetTotChargedEt() const {
+        return fTotChargedEt;
+    }
 
     /** Total charged Et in the event within the acceptance cuts */
-    Double_t GetTotChargedEtAcc() const { return fTotChargedEtAcc; }
+    Double_t GetTotChargedEtAcc() const {
+        return fTotChargedEtAcc;
+    }
+
+    void SetTPCOnlyTrackCuts(const AliESDtrackCuts *cuts) {
+        fEsdtrackCutsTPC = (AliESDtrackCuts *) cuts;
+    }
 
-    void SetTPCOnlyTrackCuts(const AliESDtrackCuts *cuts){ fEsdtrackCutsTPC = (AliESDtrackCuts *) cuts;}
-    
     /** Set the centrality object */
-    void SetCentralityObject(AliCentrality *cent) { fCentrality = cent; }
-    
+    void SetCentralityObject(AliCentrality *cent) {
+        fCentrality = cent;
+    }
+
     /** Get contribution from non-removed charged particles */
-    virtual Double_t GetChargedContribution(Int_t /*clusterMultiplicity*/) {return 0;}
+    virtual Double_t GetChargedContribution(Int_t /*clusterMultiplicity*/) {
+        return 0;
+    }
 
     /** Get contribution from non-removed neutral particles */
-    virtual Double_t GetNeutralContribution(Int_t /*clusterMultiplicity*/) {return 0;}
-    
+    virtual Double_t GetNeutralContribution(Int_t /*clusterMultiplicity*/) {
+        return 0;
+    }
+
     /** Get contribution from removed gammas */
-    virtual Double_t GetGammaContribution(Int_t /*clusterMultiplicity*/) {return 0;}
+    virtual Double_t GetGammaContribution(Int_t /*clusterMultiplicity*/) {
+        return 0;
+    }
 
-    void MakeSparseHistograms(){fMakeSparse=kTRUE;}
+    void MakeSparseHistograms() {
+        fMakeSparse=kTRUE;
+    }
+    
+    AliAnalysisEtCuts * GetCuts() const { return fCuts; }
 
 protected:
 
     //AliAnalysisEtCuts *fCuts; // keeper of basic cuts
     Double_t CalculateTransverseEnergy(AliESDCaloCluster *cluster);
 
-    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 TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge, AliVEvent *e) const;
     
-    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 */    
+    virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &/*cluster*/) const { return kFALSE; }
+    
+    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 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 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) */
 
@@ -128,13 +164,13 @@ protected:
     Double_t fNeutronEt; /** Et of neutrons (MC only for now) */
     Double_t fAntiNeutronEt; /** Et of anti-neutrons (MC only for now) */
     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 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 */
-    
+
     Float_t fEnergyDeposited; /** Energy deposited in calorimeter */
     Float_t fEnergyTPC; /** Energy measured in TPC */
     Short_t fCharge; /** Charge of the particle */
@@ -143,33 +179,37 @@ protected:
     Bool_t fTrackPassedCut; /** The track is accepted by ESDTrackCuts */
 
     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     
-    
+
+    /** 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 
-    
+    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 
-    
+    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 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
 
@@ -177,7 +217,7 @@ protected:
     TH1F *fHistEt; //Et spectrum
 
     /** The full charged Et spectrum measured */
-    TH1F *fHistChargedEt; //Charged Et spectrum 
+    TH1F *fHistChargedEt; //Charged Et spectrum
 
     /** The full neutral Et spectrum measured */
     TH1F *fHistNeutralEt; //Neutral Et spectrum
@@ -205,7 +245,7 @@ protected:
     TH2F *fHistPhivsPtNeg; //phi vs pT plot for negative tracks
 
     /* PID plots */
-    TH1F *fHistBaryonEt; /** Et of identified baryons */    
+    TH1F *fHistBaryonEt; /** Et of identified baryons */
     TH1F *fHistAntiBaryonEt; /** Et of identified anti-baryons */
     TH1F *fHistMesonEt; /** Et of identified mesons */
 
@@ -217,17 +257,17 @@ protected:
     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 *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
@@ -237,36 +277,41 @@ protected:
     static Int_t fgNumOfEBins;//number of pt bins
     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 */
+    /** 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;
 
 private:
     //Declare private to avoid compilation warning
index 4ce2a71..0f5482a 100644 (file)
@@ -21,7 +21,7 @@ AliAnalysisEtCuts::AliAnalysisEtCuts() :
   TNamed()
                                    //
   ,fCommonEtaCut(0.5)
-  ,fCommonClusterEnergyCut(0.1)
+  ,fCommonClusterEnergyCut(0.15)
   ,fCommonTrackPtCut(0.0)
   ,fCommonSingleCell(1)
   ,fEmcalTrackDistanceCut(15.0)
@@ -30,6 +30,8 @@ AliAnalysisEtCuts::AliAnalysisEtCuts() :
   ,fPhosTrackDistanceCut(10.0)  
   ,fPhosTrackDxCut(8.0)
   ,fPhosTrackDzCut(3.0)
+  ,fPhosTrackRCut(2.0)
+  ,fPhosBadDistanceCut(3.0)
   
   ,fGeometryPhosEtaAccCut(0.12)
   ,fGeometryPhosPhiAccMinCut(260.0)
@@ -51,7 +53,7 @@ AliAnalysisEtCuts::AliAnalysisEtCuts() :
   ,fReconstructedPidCut(0.0)
                                    //
   ,fReconstructedPhosClusterType(-1)
-  ,fReconstructedPhosClusterEnergyCut(0.1)
+  ,fReconstructedPhosClusterEnergyCut(0.15)
   ,fReconstructedPhosSingleCellEnergyCut(0.5)
   ,fReconstructedPhosTrackDistanceTightCut(3.0)
   ,fReconstructedPhosTrackDistanceMediumCut(5.0)
index bb00ebf..e22a846 100644 (file)
@@ -46,14 +46,18 @@ class AliAnalysisEtCuts : public TNamed
   Int_t GetReconstructedNItsClustersCut() const { return fReconstructedNItsClustersCut; }
   Double_t GetReconstructedPidCut() const { return fReconstructedPidCut; }
   // ReconstructedPhos
-  Char_t GetReconstructedPhosClusterType() const { return fReconstructedPhosClusterType; }
+  Char_t GetPhosClusterType() const { return fReconstructedPhosClusterType; }
   Double_t GetReconstructedPhosClusterEnergyCut() const { return fReconstructedPhosClusterEnergyCut; }
   Double_t GetReconstructedPhosSingleCellEnergyCut() const { return fReconstructedPhosSingleCellEnergyCut; }
   Double_t GetPhosTrackDistanceCut() const { return fPhosTrackDistanceCut; }
   Double_t GetPhosTrackDxCut() const { return fPhosTrackDxCut; }
   Double_t GetPhosTrackDzCut() const { return fPhosTrackDzCut; }
+  Double_t GetPhosTrackRCut() const { return fPhosTrackRCut; }
+  
+  Double_t GetPhosBadDistanceCut() const { return fPhosBadDistanceCut; }
+  
   // ReconstructedEmcal
-  Char_t GetReconstructedEmcalClusterType() const { return fReconstructedEmcalClusterType; }
+  Char_t GetEmcalClusterType() const { return fReconstructedEmcalClusterType; }
   Double_t GetReconstructedEmcalClusterEnergyCut() const { return fReconstructedEmcalClusterEnergyCut; }
   Double_t GetReconstructedEmcalSingleCellEnergyCut() const { return fReconstructedEmcalSingleCellEnergyCut; }
   Double_t GetEmcalTrackDistanceCut() const { return fEmcalTrackDistanceCut; }
@@ -83,6 +87,8 @@ class AliAnalysisEtCuts : public TNamed
   Double_t GetHistMinParticlePt() const { return fHistMinParticlePt; }
   Double_t GetHistMaxParticlePt() const { return fHistMaxParticlePt; }
   
+  
+  
   Short_t GetDetectorPhos() const { return fgkDetectorPhos; }
   Short_t GetDetectorEmcal() const { return fgkDetectorEmcal; }
 
@@ -118,6 +124,9 @@ class AliAnalysisEtCuts : public TNamed
   void SetPhosTrackDistanceCut(const Double_t val) { fPhosTrackDistanceCut = val; }
   void SetPhosTrackDxCut(const Double_t val) { fPhosTrackDxCut = val; }
   void SetPhosTrackDzCut(const Double_t val) { fPhosTrackDzCut = val; }
+  void SetPhosTrackRCut(const Double_t val) { fPhosTrackRCut = val; }
+  
+  void SetPhosBadDistanceCut(const Double_t val) { fPhosBadDistanceCut = val; }
   
   // ReconstructedEmcal
   void SetReconstructedEmcalClusterType(const Char_t val) { fReconstructedEmcalClusterType = val; }
@@ -162,7 +171,10 @@ class AliAnalysisEtCuts : public TNamed
   Double_t fPhosTrackDistanceCut; // PHOS track distance  
   Double_t fPhosTrackDxCut; // PHOS track distance in x
   Double_t fPhosTrackDzCut; // PHOS track distance  in z
-  
+  Double_t fPhosTrackRCut; // PHOS track distance  in r (using the parametrized track distance)
+ Double_t fPhosBadDistanceCut; // PHOS distance to bad channel 
   // GeometryPhos
   Double_t fGeometryPhosEtaAccCut; // PHOS Eta Acc cut
   Double_t fGeometryPhosPhiAccMinCut; // PHOS Phi Acc Min cut
index 3408dcc..379a825 100644 (file)
@@ -24,7 +24,9 @@
 #include "AliLog.h"
 #include <iostream>
 #include <AliCentrality.h>
-
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSGeometry.h"
+#include "TFile.h"
 using namespace std;
 
 ClassImp(AliAnalysisEtMonteCarlo);
@@ -49,6 +51,7 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fHistEtNonRemovedKaonPlus(0)
         ,fHistEtNonRemovedKaonMinus(0)
         ,fHistEtNonRemovedK0s(0)
+       ,fHistEtNonRemovedK0L(0)
         ,fHistEtNonRemovedLambdas(0)
         ,fHistEtNonRemovedElectrons(0)
         ,fHistEtNonRemovedPositrons(0)
@@ -58,9 +61,16 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fHistEtNonRemovedAntiNeutrons(0)
         ,fHistEtNonRemovedGammas(0)
         ,fHistEtNonRemovedGammasFromPi0(0)
-        ,fHistEtRemovedGammas(0)
+        
+       ,fHistEtRemovedGammas(0)
         ,fHistEtRemovedNeutrons(0)
         ,fHistEtRemovedAntiNeutrons(0)
+       
+       ,fHistEtRemovedCharged(0)
+       ,fHistEtRemovedNeutrals(0)
+       ,fHistEtNonRemovedCharged(0)
+       ,fHistEtNonRemovedNeutrals(0)
+       
         ,fHistMultNonRemovedProtons(0)
         ,fHistMultNonRemovedAntiProtons(0)
         ,fHistMultNonRemovedPiPlus(0)
@@ -68,6 +78,7 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fHistMultNonRemovedKaonPlus(0)
         ,fHistMultNonRemovedKaonMinus(0)
         ,fHistMultNonRemovedK0s(0)
+       ,fHistMultNonRemovedK0L(0)
         ,fHistMultNonRemovedLambdas(0)
         ,fHistMultNonRemovedElectrons(0)
         ,fHistMultNonRemovedPositrons(0)
@@ -76,9 +87,17 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,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)
@@ -94,7 +113,8 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fEtNonRemovedPiMinus(0)
         ,fEtNonRemovedKaonPlus(0)
         ,fEtNonRemovedKaonMinus(0)
-        ,fEtNonRemovedK0s(0)
+        ,fEtNonRemovedK0S(0)
+       ,fEtNonRemovedK0L(0)
         ,fEtNonRemovedLambdas(0)
         ,fEtNonRemovedElectrons(0)
         ,fEtNonRemovedPositrons(0)
@@ -104,6 +124,22 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,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)
@@ -114,6 +150,7 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fMultNonRemovedKaonPlus(0)
         ,fMultNonRemovedKaonMinus(0)
         ,fMultNonRemovedK0s(0)
+       ,fMultNonRemovedK0L(0)
         ,fMultNonRemovedLambdas(0)
         ,fMultNonRemovedElectrons(0)
         ,fMultNonRemovedPositrons(0)
@@ -122,9 +159,26 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,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)
@@ -150,6 +204,11 @@ AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
         ,fEnergyChargedRemoved(0)
         ,fEnergyChargedNotRemoved(0)
         ,fEnergyNeutralNotRemoved(0)
+       ,fNClusters(0)
+       ,fGeoUtils(0)
+       ,fBadMapM2(0)
+       ,fBadMapM3(0)
+       ,fBadMapM4(0)
 
 {
     fTrackDistanceCut = 7.0;
@@ -172,6 +231,7 @@ AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
        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
@@ -194,6 +254,7 @@ AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
        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 
@@ -277,17 +338,7 @@ 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();
-        /*
-          printf("Hijing: ImpactParameter %g ReactionPlaneAngle %g \n",
-          hijingGenHeader->ImpactParameter(), hijingGenHeader->ReactionPlaneAngle());
-          printf("HardScatters %d ProjecileParticipants %d TargetParticipants %d\n",
-          hijingGenHeader->HardScatters(), hijingGenHeader->ProjectileParticipants(), hijingGenHeader->TargetParticipants());
-          printf("ProjSpectatorsn %d ProjSpectatorsp %d TargSpectatorsn %d TargSpectatorsp %d\n",
-          hijingGenHeader->ProjSpectatorsn(), hijingGenHeader->ProjSpectatorsp(), hijingGenHeader->TargSpectatorsn(), hijingGenHeader->TargSpectatorsp());
-          printf("NN %d NNw %d NwN %d, NwNw %d\n",
-          hijingGenHeader->NN(), hijingGenHeader->NNw(), hijingGenHeader->NwN(), hijingGenHeader->NwNw());
-        */
-    }
+      }
 
     /* // placeholder if we want to get some Pythia info later
        AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
@@ -296,7 +347,6 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
        pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
        }
     */
-
     // Let's play with the stack!
     AliStack *stack = event->Stack();
 
@@ -320,34 +370,15 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 //         Int_t iPartLastMom = part->GetMother(1);
 
         TParticlePDG *pdg = part->GetPDG(0);
-        //TParticlePDG *pdgMom = 0;
-//         TParticlePDG *pdgMomLast = 0;
 
-        if (!pdg)
+       if (!pdg)
         {
             //Printf("ERROR: Could not get particle PDG %d", iPart);
             continue;
         }
 
-//         if (iPartMom>0)
-//         {
-//             partMom = stack->Particle(iPartMom);
-//             //pdgMom = partMom->GetPDG(0);
-//         }
-
-//         if (iPartLastMom>0)
-//         {
-//             partMomLast = stack->Particle(iPartLastMom);
-//             pdgMomLast = partMomLast->GetPDG(0);
-//         }
-
-
         Double_t particleMassPart = 0; //The mass part in the Et calculation for this particle
-
-        // Check if it is a primary particle
-        //if (!stack->IsPhysicalPrimary(iPart)) continue;
-
-        //printf("MC: iPart %03d eta %4.3f phi %4.3f code %d charge %g \n", iPart, part->Eta(), part->Phi(), pdg->PdgCode(), pdg->Charge()); // tmp/debug printout
+       Int_t code = pdg->PdgCode();
 
         // Check for reasonable (for now neutral and singly charged) charge on the particle
         //TODO:Maybe not only singly charged?
@@ -360,24 +391,24 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
             {
                 // calculate E_T
                 if (
-                    TMath::Abs(pdg->PdgCode()) == fgProtonCode ||
-                    TMath::Abs(pdg->PdgCode()) == fgNeutronCode ||
-                    TMath::Abs(pdg->PdgCode()) == fgLambdaCode ||
-                    TMath::Abs(pdg->PdgCode()) == fgXiCode ||
-                    TMath::Abs(pdg->PdgCode()) == fgXi0Code ||
-                    TMath::Abs(pdg->PdgCode()) == fgOmegaCode
+                    TMath::Abs(code) == fgProtonCode ||
+                    TMath::Abs(code) == fgNeutronCode ||
+                    TMath::Abs(code) == fgLambdaCode ||
+                    TMath::Abs(code) == fgXiCode ||
+                    TMath::Abs(code) == fgXi0Code ||
+                    TMath::Abs(code) == fgOmegaCode
                 )
                 {
-                    if (pdg->PdgCode() > 0) {
+                    if (code > 0) {
                         particleMassPart = - protonMass;
                     }
-                    if (pdg->PdgCode() < 0) {
+                    if (code < 0) {
                         particleMassPart = protonMass;
                     }
                 }
                 Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
 
-                fSparseTracks[0] = pdg->PdgCode();
+                fSparseTracks[0] = code;
                 fSparseTracks[1] = pdg->Charge()*3;
                 fSparseTracks[2] = pdg->Mass();
                 fSparseTracks[3] = et;
@@ -386,14 +417,14 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                 fSparseHistTracks->Fill(fSparseTracks);
 
                 // Fill up total E_T counters for each particle species
-                if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+                if (code == fgProtonCode || code == fgAntiProtonCode)
                 {
                     fProtonEt += et;
                 }
-                if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+                if (code == fgPiPlusCode || code == fgPiMinusCode)
                 {
                     fPionEt += et;
-                    if (pdg->PdgCode() == fgPiPlusCode)
+                    if (code == fgPiPlusCode)
                     {
                         fPiPlusMult++;
                     }
@@ -402,39 +433,40 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                         fPiMinusMult++;
                     }
                 }
-                if (pdg->PdgCode() == fgGammaCode)
+                if (code == fgGammaCode)
                 {
                     fPiZeroMult++;
                 }
-                if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+                if (code == fgKPlusCode || code == fgKMinusCode)
                 {
                     fChargedKaonEt += et;
                 }
-                if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+                if (code == fgMuPlusCode || code == fgMuMinusCode)
                 {
                     fMuonEt += et;
                 }
-                if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+                if (code == fgEPlusCode || code == fgEMinusCode)
                 {
                     fElectronEt += et;
                 }
-
                 // some neutrals also
-                if (pdg->PdgCode() == fgNeutronCode)
+                if (code == fgNeutronCode)
                 {
                     fNeutronEt += et;
                 }
-                if (pdg->PdgCode() == fgAntiNeutronCode)
+                if (code == fgAntiNeutronCode)
                 {
                     fAntiNeutronEt += et;
                 }
-                if (pdg->PdgCode() == fgGammaCode)
+                if (code == fgGammaCode)
                 {
                     fGammaEt += et;
                 }
 
                 // Neutral particles
-                if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+                //if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+                
+                if(code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
                 {
 
                     if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
@@ -443,11 +475,14 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                     if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
                     {
                         fNeutralMultiplicity++;
-                        //std::cout << pdg->PdgCode() << ", " << iPart << ", " <<  part->GetMother(0) << ", " << stack->IsPhysicalPrimary(iPart) << ", " << part->GetNDaughters() << ", " << part->Energy() << ", " << part->Eta() << ", " << part->Phi() << std::endl;
+                        //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;
-                        if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
-                        if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
+                       fTotNeutralEtAcc += et;
+                       fTotEtAcc += et;
+                        //if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEtAcc += et;
+                        //if (et > fCuts->GetCommonClusterEnergyCut()) fTotEtAcc += et;
+                       
                         if (part->Energy() > 0.05) partCount++;
                     }
                 }
@@ -463,24 +498,24 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
                         fTotChargedEtAcc += et;
                         fTotEtAcc += et;
 
-                        if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+                        if (code == fgProtonCode || code == fgAntiProtonCode)
                         {
                             fProtonEtAcc += et;
                         }
-                        if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+                        if (code == fgPiPlusCode || code == fgPiMinusCode)
                         {
                             fPionEtAcc += et;
                         }
-                        if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+                        if (code == fgKPlusCode || code == fgKMinusCode)
                         {
                             fChargedKaonEtAcc += et;
                         }
-                        if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+                        if (code == fgMuPlusCode || code == fgMuMinusCode)
                         {
                             fMuonEtAcc += et;
                         }
                         
-                        if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+                        if (code == fgEPlusCode || code == fgEMinusCode)
                         {
                           fElectronEtAcc += et;
                         }
@@ -499,9 +534,7 @@ 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;
 }
@@ -513,7 +546,17 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* 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;
@@ -545,7 +588,8 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
         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++)
     {
@@ -561,18 +605,23 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
             fTrackMultInAcc++;
         }
     }
-
+}
     Int_t nCluster = caloClusters->GetEntries();
-
+    fNClusters = 0;
     // loop the clusters
     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
     {
+       
         AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
         //Float_t caloE = caloCluster->E();
+       
+       if(caloCluster->GetType() != fClusterType) continue;
+       if(AliAnalysisEtMonteCarlo::TooCloseToBadChannel(*caloCluster))continue;
 
+       fNClusters++;
         UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
-        TParticle *part  = stack->Particle(iPart);
-        TParticle *partMom = 0;
+        TParticle *part  =  stack->Particle(iPart);
+        //TParticle *partMom = 0;
 
         if (!part)
         {
@@ -583,41 +632,54 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
         // compare MC and Rec energies for all particles
         //fHistAllERecEMC->Fill(part->Energy(),caloE);
 
-        TParticlePDG *pdg = part->GetPDG(0);
-        TParticlePDG *pdgMom = 0;
-
-        if (!pdg)
-        {
-            Printf("ERROR: Could not get particle PDG %d", iPart);
-            continue;
-        }
+        //TParticlePDG *pdgMom = 0;
 
+/*
         Int_t iPartMom = part->GetMother(0);
 
         if (iPartMom>0)
         {
             partMom = stack->Particle(iPartMom);
             pdgMom = partMom->GetPDG(0);
-        }
-
-        // Check if it is a primary particle
-        //if (!stack->IsPhysicalPrimary(iPart)) continue;
+        }*/
+       //int mcClusterCharge = stack->Particle(iPart)->GetPDG()->Charge();
         if (!stack->IsPhysicalPrimary(iPart)) // check whether particle is primary. we keep secondary electron and gamma for testing.
         {
-            if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
-                continue;
+   //         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;
+               }
+               
         } // end of primary particle check
-
-        if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
+        TParticlePDG *pdg = part->GetPDG();
+       Int_t code = pdg->PdgCode();
+       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;
 
         Double_t clEt = CalculateTransverseEnergy(caloCluster);
+//     if(code == fgK0SCode) std::cout << "K0 energy: " << caloCluster->E() << std::endl;
         if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
         Float_t pos[3];
 
         caloCluster->GetPosition(pos);
         TVector3 cp(pos);
 
-        fSparseClusters[0] = pdg->PdgCode();
+        fSparseClusters[0] = code;
         fSparseClusters[1] = pdg->Charge()/3;
         fSparseClusters[2] = pdg->Mass();
         fSparseClusters[3] = clEt;
@@ -630,16 +692,19 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
         fSparseClusters[10] = 0;
         fSparseHistClusters->Fill(fSparseClusters);
 
+       
         //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
-
-        if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
+       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 (pdg->Charge() != 0)
             {
-                //std::cout << "Removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                //std::cout << "Removing charged: " << code << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
                 fChargedRemoved++;
-                fEnergyChargedRemoved += caloCluster->E();
+                //fEnergyChargedRemoved += caloCluster->E();
+               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());
@@ -649,31 +714,33 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
             }
             else
             {
-                //std::cout << "Removing neutral: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                //std::cout << "Removing neutral: " << 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;
                 fNeutralRemoved++;
-                fEnergyNeutralRemoved += caloCluster->E();
+                //fEnergyNeutralRemoved += caloCluster->E();
+               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());
                 fHistDxDzRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
 
-                if (pdg->PdgCode() == fgGammaCode)
+                if (code == fgGammaCode || code == fgPi0Code || code == fgEtaCode)
                 {
                     fEtRemovedGammas += clEt;
                     fMultRemovedGammas++;
                 }
-                else if (pdg->PdgCode() == fgNeutronCode)
+                else if (code == fgNeutronCode)
                 {
                     fEtRemovedNeutrons += clEt;
                     fMultRemovedNeutrons++;
                 }
-                else if (pdg->PdgCode() == fgAntiNeutronCode)
+                else if (code == fgAntiNeutronCode)
                 {
                     fEtRemovedAntiNeutrons += clEt;
                     fMultRemovedAntiNeutrons++;
                 }
-                //else std::cout << "Hmm, what is this (neutral: " << pdg->PdgCode() << std::endl;
+                
+                //else std::cout << "Hmm, what is this (neutral: " << code << std::endl;
             }
         }
         else
@@ -681,52 +748,53 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 
             if (pdg->Charge() != 0)
             {
-                //std::cout << "Not removing charged: " << pdg->PdgCode() << ", 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 += caloCluster->E();
+               fEnergyChargedNotRemoved += clEt;
                 fHistRemovedOrNot->Fill(2.0, fCentClass);
                 //std::cout << fHistDecayVertexNonRemovedCharged << std::endl;
                 //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
                 //fHistDecayVertexNonRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
                 fHistDxDzNonRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
-                if (pdg->PdgCode() == fgProtonCode)
+                if (code == fgProtonCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedProtons += clEt;
                     fMultNonRemovedProtons++;
                 }
-                else if (pdg->PdgCode() == fgAntiProtonCode)
+                else if (code == fgAntiProtonCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedAntiProtons += clEt;
                     fMultNonRemovedAntiProtons++;
                 }
-                else if (pdg->PdgCode() == fgPiPlusCode)
+                else if (code == fgPiPlusCode)
                 {
                     //std::cout << "PI+" <<  clEt << std::endl;
                     fEtNonRemovedPiPlus += clEt;
                     fMultNonRemovedPiPlus++;
                 }
-                else if (pdg->PdgCode() == fgPiMinusCode)
+                else if (code == fgPiMinusCode)
                 {
                     // std::cout << "PI-"  << clEt << std::endl;
                     fEtNonRemovedPiMinus += clEt;
                     fMultNonRemovedPiMinus++;
                 }
-                else if (pdg->PdgCode() == fgKPlusCode)
+                else if (code == fgKPlusCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedKaonPlus += clEt;
                     fMultNonRemovedKaonPlus++;
                 }
-                else if (pdg->PdgCode() == fgKMinusCode)
+                else if (code == fgKMinusCode)
                 {
                     //std::cout << clEt << std::endl;
                     fEtNonRemovedKaonMinus += clEt;
                     fMultNonRemovedKaonMinus++;
                 }
-                else if (pdg->PdgCode() == fgEPlusCode)
+                else if (code == fgEPlusCode)
                 {
                     //std::cout << clEt << std::endl;
                     if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
@@ -735,7 +803,7 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                         fMultNonRemovedPositrons++;
                     }
                 }
-                else if (pdg->PdgCode() == fgEMinusCode)
+                else if (code == fgEMinusCode)
                 {
                     //std::cout << clEt << std::endl;
                     if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
@@ -744,15 +812,13 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
                         fMultNonRemovedElectrons++;
                     }
                 }
-                else if (pdg->PdgCode() == fgMuPlusCode)
+                else if (code == fgMuPlusCode)
                 {
-                    std::cout << clEt << std::endl;
                     fEtNonRemovedMuPlus += clEt;
                     fMultNonRemovedMuPlus++;
                 }
-                else if (pdg->PdgCode() == fgMuMinusCode)
+                else if (code == fgMuMinusCode)
                 {
-                    std::cout << clEt << std::endl;
                     fEtNonRemovedMuMinus += clEt;
                     fMultNonRemovedMuMinus++;
                 }
@@ -760,50 +826,73 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
             }
             else
             {
-                //std::cout << "Accepted: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
-                //std::cout << "Not removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                //std::cout << "Accepted: " << 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() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
                 fNeutralNotRemoved++;
-                fEnergyNeutralNotRemoved += caloCluster->E();
+               fEnergyNeutralNotRemoved += clEt;
                 fHistRemovedOrNot->Fill(3.0, fCentClass);
-                //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
-                //fHistDecayVertexNonRemovedNeutral->Fill(part->Vy(), part->Vx(), part->Vz());
                 fHistDxDzNonRemovedNeutral->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
-                if (pdg->PdgCode() == fgGammaCode)
+                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 (pdgMom)
+  //                  {
+    //                    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 (pdg->PdgCode() == fgNeutronCode)
+                else if(TMath::Abs(code) == fgPi0Code)
+               {
+                 fEtNonRemovedGammasFromPi0 += clEt;
+                 fMultNonRemovedGammas++;
+               }
+                else if(TMath::Abs(code) == fgEtaCode)
+               {
+                 fEtNonRemovedGammasFromPi0 += clEt;
+                 fMultNonRemovedGammas++;
+               }
+                else if(TMath::Abs(code) == 331)
+               {
+                 fEtNonRemovedGammasFromPi0 += clEt;
+                 fMultNonRemovedGammas++;
+               }
+                else if (code == fgNeutronCode)
                 {
                     fEtNonRemovedNeutrons += clEt;
                     fMultNonRemovedNeutrons++;
                 }
-                else if (pdg->PdgCode() == fgAntiNeutronCode)
+                else if (code == fgAntiNeutronCode)
                 {
                     fEtNonRemovedAntiNeutrons += clEt;
                     fMultNonRemovedAntiNeutrons++;
                 }
-                else if (pdg->PdgCode() == fgK0LCode || pdg->PdgCode() == fgK0SCode)
+                //else if (code == fgK0LCode || pdgMom->PdgCode() == fgK0SCode)
+                else if (code == fgK0SCode)
                 {
-                    fEtNonRemovedK0s += clEt;
+                   //std::cout << "K0 with energy: " << clEt << std::endl;
+                    fEtNonRemovedK0S += clEt;
                     fMultNonRemovedK0s++;
 
                 }
-                else if (TMath::Abs(pdg->PdgCode()) == fgLambdaCode)
+                else if(TMath::Abs(code) == fgK0LCode)
+               {
+                 
+                 fEtNonRemovedK0L += clEt;
+                 fMultNonRemovedK0L++;
+               }
+               
+                else if (TMath::Abs(code) == fgLambdaCode)
                 {
                     fEtNonRemovedLambdas += clEt;
                     fMultNonRemovedLambdas++;
                 }
-                else std::cout << "Hmm, what is this (neutral not removed): " << pdg->PdgCode() << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
+                  
+                else std::cout << "Hmm, what is this (neutral not removed): " << code << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
             }
         }
     } // end of loop over clusters
@@ -827,6 +916,12 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 void AliAnalysisEtMonteCarlo::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()
@@ -844,7 +939,8 @@ void AliAnalysisEtMonteCarlo::ResetEventValues()
     fEtNonRemovedPiMinus = 0;
     fEtNonRemovedKaonPlus = 0;
     fEtNonRemovedKaonMinus = 0;
-    fEtNonRemovedK0s = 0;
+    fEtNonRemovedK0S = 0;
+    fEtNonRemovedK0L = 0;
     fEtNonRemovedLambdas = 0;
     fEtNonRemovedElectrons = 0;
     fEtNonRemovedPositrons = 0;
@@ -855,6 +951,22 @@ void AliAnalysisEtMonteCarlo::ResetEventValues()
     fEtNonRemovedGammas = 0;
     fEtNonRemovedGammasFromPi0 = 0;
 
+    fEtRemovedProtons = 0;
+    fEtRemovedAntiProtons = 0;
+    fEtRemovedPiPlus = 0;
+    fEtRemovedPiMinus = 0;
+    fEtRemovedKaonPlus = 0;
+    fEtRemovedKaonMinus = 0;
+    fEtRemovedK0s = 0;
+    fEtRemovedK0L = 0;
+    fEtRemovedLambdas = 0;
+    fEtRemovedElectrons = 0;
+    fEtRemovedPositrons = 0;
+    fEtRemovedMuPlus = 0;
+    fEtRemovedMuMinus = 0;
+    fEtRemovedNeutrons = 0;
+
+    fEtRemovedGammasFromPi0 = 0;
     fEtRemovedGammas = 0;
     fEtRemovedNeutrons = 0;
     fEtRemovedAntiNeutrons = 0;
@@ -866,6 +978,7 @@ void AliAnalysisEtMonteCarlo::ResetEventValues()
     fMultNonRemovedKaonPlus = 0;
     fMultNonRemovedKaonMinus = 0;
     fMultNonRemovedK0s = 0;
+    fMultNonRemovedK0L = 0;
     fMultNonRemovedLambdas = 0;
     fMultNonRemovedElectrons = 0;
     fMultNonRemovedPositrons = 0;
@@ -875,10 +988,35 @@ void AliAnalysisEtMonteCarlo::ResetEventValues()
     fMultNonRemovedAntiNeutrons = 0;
     fMultNonRemovedGammas = 0;
 
+    fMultRemovedProtons = 0;
+    fMultRemovedAntiProtons = 0;
+    fMultRemovedPiPlus = 0;
+    fMultRemovedPiMinus = 0;
+    fMultRemovedKaonPlus = 0;
+    fMultRemovedKaonMinus = 0;
+    fMultRemovedK0s = 0;
+    fMultRemovedK0L = 0;
+    fMultRemovedLambdas = 0;
+    fMultRemovedElectrons = 0;
+    fMultRemovedPositrons = 0;
+    fMultRemovedMuPlus = 0;
+    fMultRemovedMuMinus = 0;
+
     fMultRemovedGammas = 0;
     fMultRemovedNeutrons = 0;
     fMultRemovedAntiNeutrons = 0;
 
+    fEnergyChargedNotRemoved = 0;
+    fEnergyChargedRemoved = 0;
+    fEnergyNeutralNotRemoved = 0;
+    fEnergyNeutralRemoved = 0;
+    
+    fChargedNotRemoved = 0;
+    fChargedRemoved = 0;
+    fNeutralNotRemoved = 0;
+    fNeutralRemoved = 0;
+    
+    
     fTrackMultInAcc = 0;
 
 }
@@ -906,6 +1044,7 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistEtNonRemovedKaonPlus = new TH2F("fHistEtNonRemovedKaonPlus", "fHistEtNonRemovedKaonPlus", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedKaonMinus = new TH2F("fHistEtNonRemovedKaonMinus", "fHistEtNonRemovedKaonMinus", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedK0s = new TH2F("fHistEtNonRemovedK0s", "fHistEtNonRemovedK0s", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedK0L = new TH2F("fHistEtNonRemovedK0L", "fHistEtNonRemovedK0L", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedLambdas = new TH2F("fHistEtNonRemovedLambdas", "fHistEtNonRemovedLambdas", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedElectrons = new TH2F("fHistEtNonRemovedElectrons", "fHistEtNonRemovedElectrons", 1500, 0, 30, 10, -0.5, 9.5);
     fHistEtNonRemovedPositrons = new TH2F("fHistEtNonRemovedPositrons", "fHistEtNonRemovedPositrons", 1500, 0, 30, 10, -0.5, 9.5);
@@ -920,7 +1059,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);
@@ -928,6 +1073,7 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistMultNonRemovedKaonPlus = new TH2F("fHistMultNonRemovedKaonPlus", "fHistMultNonRemovedKaonPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedKaonMinus = new TH2F("fHistMultNonRemovedKaonMinus", "fHistMultNonRemovedKaonMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedK0s = new TH2F("fHistMultNonRemovedK0s", "fHistMultNonRemovedK0s", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    fHistMultNonRemovedK0L = new TH2F("fHistMultNonRemovedK0L", "fHistMultNonRemovedK0L", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedLambdas = new TH2F("fHistMultNonRemovedLambdas", "fHistMultNonRemovedLambdas", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedElectrons = new TH2F("fHistMultNonRemovedElectrons", "fHistMultNonRemovedElectrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedPositrons = new TH2F("fHistMultNonRemovedPositrons", "fHistMultNonRemovedPositrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
@@ -936,11 +1082,24 @@ void AliAnalysisEtMonteCarlo::CreateHistograms()
     fHistMultNonRemovedNeutrons = new TH2F("fHistMultNonRemovedNeutrons", "fHistMultNonRemovedNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
     fHistMultNonRemovedAntiNeutrons = new TH2F("fHistMultNonRemovedAntiNeutrons", "fHistMultNonRemovedAntiNeutrons", 100, -0.5, 99.5, 10, -0.5, 9.5);
 
-    fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    fHistMultNonRemovedGammas = new  TH2F("fHistMultNonRemovedGammas", "fHistMultNonRemovedGammas", 100, -0.5, 99.5, 100, -0.5, 99.5);
 
-    fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    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);
+
+    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);
+
+    fHistMultNonRemovedCharged = new  TH2F("fHistMultNonRemovedCharged", "fHistMultNonRemovedCharged", 100, -0.5, 99.5, 100, -0.5, 99.5);
+    fHistMultNonRemovedNeutrals = new  TH2F("fHistMultNonRemovedNeutrals", "fHistMultNonRemovedNeutrals", 100, -0.5, 99.5, 100, -0.5, 99.5);
 
     fHistTrackMultvsNonRemovedCharged = new TH2F("fHistTrackMultvsNonRemovedCharged", "fHistTrackMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
     fHistTrackMultvsNonRemovedNeutral = new TH2F("fHistTrackMultvsNonRemovedNeutral", "fHistTrackMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
@@ -978,6 +1137,7 @@ void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
     list->Add(fHistEtNonRemovedKaonPlus);
     list->Add(fHistEtNonRemovedKaonMinus);
     list->Add(fHistEtNonRemovedK0s);
+    list->Add(fHistEtNonRemovedK0L);
     list->Add(fHistEtNonRemovedLambdas);
     list->Add(fHistEtNonRemovedElectrons);
     list->Add(fHistEtNonRemovedPositrons);
@@ -992,6 +1152,11 @@ void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
     list->Add(fHistEtRemovedNeutrons);
     list->Add(fHistEtRemovedAntiNeutrons);
 
+    list->Add(fHistEtRemovedCharged);
+    list->Add(fHistEtRemovedNeutrals);
+
+    list->Add(fHistEtNonRemovedCharged);
+    list->Add(fHistEtNonRemovedNeutrals);
 
     list->Add(fHistMultNonRemovedProtons);
     list->Add(fHistMultNonRemovedAntiProtons);
@@ -1000,6 +1165,7 @@ void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
     list->Add(fHistMultNonRemovedKaonPlus);
     list->Add(fHistMultNonRemovedKaonMinus);
     list->Add(fHistMultNonRemovedK0s);
+    list->Add(fHistMultNonRemovedK0L);
     list->Add(fHistMultNonRemovedLambdas);
     list->Add(fHistMultNonRemovedElectrons);
     list->Add(fHistMultNonRemovedPositrons);
@@ -1013,6 +1179,12 @@ void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
     list->Add(fHistMultRemovedNeutrons);
     list->Add(fHistMultRemovedAntiNeutrons);
 
+    list->Add(fHistMultRemovedCharged);
+    list->Add(fHistMultRemovedNeutrals);
+
+    list->Add(fHistMultNonRemovedCharged);
+    list->Add(fHistMultNonRemovedNeutrals);
+
     list->Add(fHistTrackMultvsNonRemovedCharged);
     list->Add(fHistTrackMultvsNonRemovedNeutral);
     list->Add(fHistTrackMultvsRemovedGamma);
@@ -1069,7 +1241,8 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
     fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
     fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
     fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
-    fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0s, fCentClass);
+    fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0S, fCentClass);
+    fHistEtNonRemovedK0L->Fill(fEtNonRemovedK0L, fCentClass);
     fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
     fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
     fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
@@ -1082,15 +1255,39 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
     fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
     fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
 
-    fHistEtRemovedGammas->Fill(fEtRemovedGammas, fCentClass);
+    fHistEtRemovedGammas->Fill(fEtRemovedGammas, fNClusters);
     fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
     fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
 
+//     fHistEtRemovedCharged->Fill(fEtRemovedAntiProtons+fEtRemovedElectrons+fEtRemovedKaonMinus+fEtRemovedKaonPlus
+//                                             +fEtRemovedMuMinus+fEtRemovedMuPlus+fEtRemovedPiMinus+fEtRemovedPiPlus+fEtRemovedPositrons
+//                                             +fEtRemovedProtons.
+//                             fCentClass);
+//     fHistEtRemovedNeutrals->Fill(fEtRemovedNeutrons+fEtRemovedAntiNeutrons, fCentClass);
+// 
+//     fHistEtNonRemovedCharged->Fill(fEtNonRemovedAntiProtons+fEtNonRemovedElectrons+fEtNonRemovedKaonMinus+fEtNonRemovedKaonPlus
+//                                             +fEtNonRemovedMuMinus+fEtNonRemovedMuPlus+fEtNonRemovedPiMinus+fEtNonRemovedPiPlus+fEtNonRemovedPositrons
+//                                             +fEtNonRemovedProtons,
+//                             fCentClass);
+//     fHistEtRemovedNeutrals->Fill(fEtNonRemovedNeutrons+fEtNonRemovedAntiNeutrons, fCentClass);
+
+    fHistEtRemovedCharged->Fill(fEnergyChargedRemoved, fNClusters);
+    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);
     fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
     fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
+    fHistMultNonRemovedK0L->Fill(fMultNonRemovedK0L, fCentClass);
     fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
     fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
     fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
@@ -1112,20 +1309,20 @@ void AliAnalysisEtMonteCarlo::FillHistograms()
                                             +fMultNonRemovedProtons);
 
     fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
-                                            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+                                            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
 
     fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
                                        fMultRemovedGammas);
 
-    fHistClusterMultvsNonRemovedCharged->Fill(fNeutralMultiplicity,
+    fHistClusterMultvsNonRemovedCharged->Fill(fNClusters,
             fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
             +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
             +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
 
-    fHistClusterMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
-            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+    fHistClusterMultvsNonRemovedNeutral->Fill(fNClusters,
+            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedK0L+fMultNonRemovedLambdas);
 
-    fHistClusterMultvsRemovedGamma->Fill(fTrackMultInAcc,
+    fHistClusterMultvsRemovedGamma->Fill(fNClusters,
                                          fMultRemovedGammas);
 
 }
@@ -1133,12 +1330,209 @@ 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;
+}
+
+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) 
+      {        
+//     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;
+}
+
+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(mother->GetPdgCode() == fgK0SCode) 
+      //{      
+//     std::cout << "!!!!!!!!!!!!!!!!! K0S !!!!!!!!!!!!!!!!!!" << std::endl;
+       //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 -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 
+    {
+      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 (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;
+
+
+}
 
 
 
index 90e9203..6880868 100644 (file)
@@ -1,6 +1,7 @@
 #ifndef ALIANALYSISETMONTECARLO_H
 #define ALIANALYSISETMONTECARLO_H
-//_________________________________________________________________________
+
+class AliPHOSGeometry;//_________________________________________________________________________
 //  Utility Class for transverse energy studies
 //  Base class for MC analysis
 //  - MC output
@@ -11,6 +12,9 @@
 #include "AliAnalysisEt.h"
 class TParticle;
 class TH3F;
+class TH2I;
+class AliPHOSGeoUtils;
+class AliStack;
 //class AliMCEvent;
 //class AliESDEvent;
 
@@ -18,13 +22,13 @@ class AliAnalysisEtMonteCarlo : public AliAnalysisEt
 {
 
 public:
-   
-  AliAnalysisEtMonteCarlo();
-  virtual ~AliAnalysisEtMonteCarlo();
+
+    AliAnalysisEtMonteCarlo();
+    virtual ~AliAnalysisEtMonteCarlo();
 
     virtual Int_t AnalyseEvent(AliVEvent* event);
-       virtual Int_t AnalyseEvent(AliVEvent* event, AliVEvent* event2);
-       //virtual Int_t AnalyseEvent(AliMCEvent* event, AliESDEvent* event2);
+    virtual Int_t AnalyseEvent(AliVEvent* event, AliVEvent* event2);
+    //virtual Int_t AnalyseEvent(AliMCEvent* event, AliESDEvent* event2);
 
     virtual void Init();
     virtual void ResetEventValues();
@@ -36,6 +40,18 @@ public:
 protected:
 
     virtual bool TrackHitsCalorimeter(TParticle *part, Double_t magField=0.5);
+    
+    
+    Int_t GetPrimMother(Int_t partIdx, AliStack *stack);
+    
+    Int_t GetK0InFamily(Int_t partIdx, AliStack *stack);
+    
+    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:
 
@@ -43,145 +59,197 @@ protected:
     Int_t fNcoll; // Ncoll, for Hijing; 1 otherwise
     Int_t fNpart; // Ncoll, for Hijing; 2 otherwise
 
-       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
-       TH3F *fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
-       
-       TH2F *fHistRemovedOrNot; // If charged/neutral particles were removed or not
-       
-       TH2F *fHistEtNonRemovedProtons; // enter comment here
-       TH2F *fHistEtNonRemovedAntiProtons; // enter comment here
-       TH2F *fHistEtNonRemovedPiPlus; // enter comment here
-       TH2F *fHistEtNonRemovedPiMinus; // enter comment here
-       TH2F *fHistEtNonRemovedKaonPlus; // enter comment here
-       TH2F *fHistEtNonRemovedKaonMinus; // enter comment here
-       TH2F *fHistEtNonRemovedK0s; // enter comment here
-       TH2F *fHistEtNonRemovedLambdas; // enter comment here
-       TH2F *fHistEtNonRemovedElectrons; // enter comment here
-       TH2F *fHistEtNonRemovedPositrons; // enter comment here
-       TH2F *fHistEtNonRemovedMuPlus; // enter comment here
-       TH2F *fHistEtNonRemovedMuMinus; // enter comment here
-       TH2F *fHistEtNonRemovedNeutrons; // enter comment here
-       TH2F *fHistEtNonRemovedAntiNeutrons; // enter comment here
-       TH2F *fHistEtNonRemovedGammas; // enter comment here
-       TH2F *fHistEtNonRemovedGammasFromPi0; // enter comment here
-       
-       TH2F *fHistEtRemovedGammas; // enter comment here
-       TH2F *fHistEtRemovedNeutrons; // enter comment here
-       TH2F *fHistEtRemovedAntiNeutrons; // enter comment here
-       
-       
-       TH2F *fHistMultNonRemovedProtons; // enter comment here 
-       TH2F *fHistMultNonRemovedAntiProtons; // enter comment here 
-       TH2F *fHistMultNonRemovedPiPlus; // enter comment here 
-       TH2F *fHistMultNonRemovedPiMinus; // enter comment here 
-       TH2F *fHistMultNonRemovedKaonPlus; // enter comment here 
-       TH2F *fHistMultNonRemovedKaonMinus; // enter comment here 
-       TH2F *fHistMultNonRemovedK0s; // enter comment here 
-       TH2F *fHistMultNonRemovedLambdas; // enter comment here 
-       TH2F *fHistMultNonRemovedElectrons; // enter comment here 
-       TH2F *fHistMultNonRemovedPositrons; // enter comment here 
-       TH2F *fHistMultNonRemovedMuPlus; // enter comment here
-       TH2F *fHistMultNonRemovedMuMinus; // enter comment here
-       TH2F *fHistMultNonRemovedNeutrons; // enter comment here
-       TH2F *fHistMultNonRemovedAntiNeutrons; // enter comment here
-       TH2F *fHistMultNonRemovedGammas; // enter comment here
-       
-       TH2F *fHistMultRemovedGammas; // enter comment here
-       TH2F *fHistMultRemovedNeutrons; // enter comment here
-       TH2F *fHistMultRemovedAntiNeutrons; // enter comment here
-       
-       TH2F *fHistTrackMultvsNonRemovedCharged; // enter comment here
-       TH2F *fHistTrackMultvsNonRemovedNeutral; // enter comment here
-       TH2F *fHistTrackMultvsRemovedGamma; // enter comment here
-       
-       TH2F *fHistClusterMultvsNonRemovedCharged; // enter comment here
-       TH2F *fHistClusterMultvsNonRemovedNeutral; // enter comment here
-       TH2F *fHistClusterMultvsRemovedGamma; // enter comment here
-       
-       TH2F *fHistMultvsNonRemovedChargedE; // enter comment here
-       TH2F *fHistMultvsNonRemovedNeutralE; // enter comment here
-       TH2F *fHistMultvsRemovedGammaE; // enter comment here
-       
-       Float_t fEtNonRemovedProtons; // enter comment here 
-       Float_t fEtNonRemovedAntiProtons; // enter comment here 
-       Float_t fEtNonRemovedPiPlus; // enter comment here 
-       Float_t fEtNonRemovedPiMinus; // enter comment here 
-       Float_t fEtNonRemovedKaonPlus; // enter comment here 
-       Float_t fEtNonRemovedKaonMinus; // enter comment here 
-       Float_t fEtNonRemovedK0s; // enter comment here 
-       Float_t fEtNonRemovedLambdas; // enter comment here 
-       Float_t fEtNonRemovedElectrons; // enter comment here 
-       Float_t fEtNonRemovedPositrons; // enter comment here 
-       Float_t fEtNonRemovedMuMinus; // enter comment here
-       Float_t fEtNonRemovedMuPlus; // enter comment here
-       Float_t fEtNonRemovedGammas; // enter comment here
-       Float_t fEtNonRemovedGammasFromPi0; // enter comment here
-       Float_t fEtNonRemovedNeutrons; // enter comment here
-       Float_t fEtNonRemovedAntiNeutrons; // enter comment here
-       
-       Float_t fEtRemovedGammas; // enter comment here
-       Float_t fEtRemovedNeutrons; // enter comment here
-       Float_t fEtRemovedAntiNeutrons; // enter comment here
-               
-       Int_t fMultNonRemovedProtons; // enter comment here 
-       Int_t fMultNonRemovedAntiProtons; // enter comment here 
-       Int_t fMultNonRemovedPiPlus; // enter comment here 
-       Int_t fMultNonRemovedPiMinus; // enter comment here 
-       Int_t fMultNonRemovedKaonPlus; // enter comment here 
-       Int_t fMultNonRemovedKaonMinus; // enter comment here 
-       Int_t fMultNonRemovedK0s; // enter comment here 
-       Int_t fMultNonRemovedLambdas; // enter comment here 
-       Int_t fMultNonRemovedElectrons; // enter comment here 
-       Int_t fMultNonRemovedPositrons; // enter comment here 
-       Int_t fMultNonRemovedMuMinus; // enter comment here
-       Int_t fMultNonRemovedMuPlus; // enter comment here
-       Int_t fMultNonRemovedGammas; // enter comment here
-       Int_t fMultNonRemovedNeutrons; // enter comment here
-       Int_t fMultNonRemovedAntiNeutrons; // enter comment here
-       
-       Int_t fMultRemovedGammas; // enter comment here
-       Int_t fMultRemovedNeutrons; // enter comment here
-       Int_t fMultRemovedAntiNeutrons; // enter comment here
-       
-       Int_t fTrackMultInAcc; // enter comment here
-       
-       
-       TH2F *fHistDxDzNonRemovedCharged; // enter comment here
-       TH2F *fHistDxDzRemovedCharged; // enter comment here
-       TH2F *fHistDxDzNonRemovedNeutral; // enter comment here
-       TH2F *fHistDxDzRemovedNeutral; // enter comment here
-       
-       TH1F *fHistPiPlusMult; // enter comment here
-       TH1F *fHistPiMinusMult; // enter comment here
-       TH1F *fHistPiZeroMult; // enter comment here
-
-       TH1F *fHistPiPlusMultAcc; // enter comment here
-       TH1F *fHistPiMinusMultAcc; // enter comment here
-       TH1F *fHistPiZeroMultAcc; // enter comment here
-       
-       Int_t fPiPlusMult; // enter comment here
-       Int_t fPiMinusMult; // enter comment here
-       Int_t fPiZeroMult; // enter comment here
-
-       Int_t fPiPlusMultAcc; // enter comment here
-       Int_t fPiMinusMultAcc; // enter comment here
-       Int_t fPiZeroMultAcc; // enter comment here
-       
-       
-       Int_t fNeutralRemoved; // number of neutral particles that where removed by track matching
-       Int_t fChargedRemoved; // number of charged particles that where removed by track matching
-       Int_t fChargedNotRemoved; // number of charged particles that were not removed
-       Int_t fNeutralNotRemoved; // number of neutral particles that were not removed
-       
-       Double_t fEnergyNeutralRemoved; // energy of neutral particles that where removed by track matching
-       Double_t fEnergyChargedRemoved; // energy of charged particles that where removed by track matching
-       Double_t fEnergyChargedNotRemoved; // energy of charged particles that were not removed
-       Double_t fEnergyNeutralNotRemoved; // energy of neutral particles that were not removed
-       
-       
- private:
+    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
+    TH3F *fHistDecayVertexRemovedNeutral; // Decay vertex for non-removed charged particles
+
+    TH2F *fHistRemovedOrNot; // If charged/neutral particles were removed or not
+
+    TH2F *fHistEtNonRemovedProtons; // enter comment here
+    TH2F *fHistEtNonRemovedAntiProtons; // enter comment here
+    TH2F *fHistEtNonRemovedPiPlus; // enter comment here
+    TH2F *fHistEtNonRemovedPiMinus; // enter comment here
+    TH2F *fHistEtNonRemovedKaonPlus; // enter comment here
+    TH2F *fHistEtNonRemovedKaonMinus; // enter comment here
+    TH2F *fHistEtNonRemovedK0s; // enter comment here
+    TH2F *fHistEtNonRemovedK0L; // enter comment here
+    TH2F *fHistEtNonRemovedLambdas; // enter comment here
+    TH2F *fHistEtNonRemovedElectrons; // enter comment here
+    TH2F *fHistEtNonRemovedPositrons; // enter comment here
+    TH2F *fHistEtNonRemovedMuPlus; // enter comment here
+    TH2F *fHistEtNonRemovedMuMinus; // enter comment here
+    TH2F *fHistEtNonRemovedNeutrons; // enter comment here
+    TH2F *fHistEtNonRemovedAntiNeutrons; // enter comment here
+    TH2F *fHistEtNonRemovedGammas; // enter comment here
+    TH2F *fHistEtNonRemovedGammasFromPi0; // enter comment here
+
+    TH2F *fHistEtRemovedGammas; // enter comment here
+    TH2F *fHistEtRemovedNeutrons; // enter comment here
+    TH2F *fHistEtRemovedAntiNeutrons; // enter comment here
+
+    TH2F *fHistEtRemovedCharged; // enter comment here
+    TH2F *fHistEtRemovedNeutrals; // enter comment here
+
+    TH2F *fHistEtNonRemovedCharged; // enter comment here
+    TH2F *fHistEtNonRemovedNeutrals; // enter comment here
+
+    TH2F *fHistMultNonRemovedProtons; // enter comment here
+    TH2F *fHistMultNonRemovedAntiProtons; // enter comment here
+    TH2F *fHistMultNonRemovedPiPlus; // enter comment here
+    TH2F *fHistMultNonRemovedPiMinus; // enter comment here
+    TH2F *fHistMultNonRemovedKaonPlus; // enter comment here
+    TH2F *fHistMultNonRemovedKaonMinus; // enter comment here
+    TH2F *fHistMultNonRemovedK0s; // enter comment here
+    TH2F *fHistMultNonRemovedK0L; // enter comment here
+    TH2F *fHistMultNonRemovedLambdas; // enter comment here
+    TH2F *fHistMultNonRemovedElectrons; // enter comment here
+    TH2F *fHistMultNonRemovedPositrons; // enter comment here
+    TH2F *fHistMultNonRemovedMuPlus; // enter comment here
+    TH2F *fHistMultNonRemovedMuMinus; // enter comment here
+    TH2F *fHistMultNonRemovedNeutrons; // enter comment here
+    TH2F *fHistMultNonRemovedAntiNeutrons; // enter comment here
+    TH2F *fHistMultNonRemovedGammas; // enter comment here
+
+    TH2F *fHistMultRemovedGammas; // enter comment here
+    TH2F *fHistMultRemovedNeutrons; // enter comment here
+    TH2F *fHistMultRemovedAntiNeutrons; // enter comment here
+
+    TH2F *fHistMultRemovedCharged; // enter comment here
+    TH2F *fHistMultRemovedNeutrals; // enter comment here
+
+    TH2F *fHistMultNonRemovedCharged; // enter comment here
+    TH2F *fHistMultNonRemovedNeutrals; // enter comment here
+
+    TH2F *fHistTrackMultvsNonRemovedCharged; // enter comment here
+    TH2F *fHistTrackMultvsNonRemovedNeutral; // enter comment here
+    TH2F *fHistTrackMultvsRemovedGamma; // enter comment here
+
+    TH2F *fHistClusterMultvsNonRemovedCharged; // enter comment here
+    TH2F *fHistClusterMultvsNonRemovedNeutral; // enter comment here
+    TH2F *fHistClusterMultvsRemovedGamma; // enter comment here
+
+    TH2F *fHistMultvsNonRemovedChargedE; // enter comment here
+    TH2F *fHistMultvsNonRemovedNeutralE; // enter comment here
+    TH2F *fHistMultvsRemovedGammaE; // enter comment here
+
+    Float_t fEtNonRemovedProtons; // enter comment here
+    Float_t fEtNonRemovedAntiProtons; // enter comment here
+    Float_t fEtNonRemovedPiPlus; // enter comment here
+    Float_t fEtNonRemovedPiMinus; // enter comment here
+    Float_t fEtNonRemovedKaonPlus; // enter comment here
+    Float_t fEtNonRemovedKaonMinus; // enter comment here
+    Float_t fEtNonRemovedK0S; // enter comment here
+    Float_t fEtNonRemovedK0L; // enter comment here
+    Float_t fEtNonRemovedLambdas; // enter comment here
+    Float_t fEtNonRemovedElectrons; // enter comment here
+    Float_t fEtNonRemovedPositrons; // enter comment here
+    Float_t fEtNonRemovedMuMinus; // enter comment here
+    Float_t fEtNonRemovedMuPlus; // enter comment here
+    Float_t fEtNonRemovedGammas; // enter comment here
+    Float_t fEtNonRemovedGammasFromPi0; // enter comment here
+    Float_t fEtNonRemovedNeutrons; // enter comment here
+    Float_t fEtNonRemovedAntiNeutrons; // enter comment here
+
+    Float_t fEtRemovedProtons; // enter comment here
+    Float_t fEtRemovedAntiProtons; // enter comment here
+    Float_t fEtRemovedPiPlus; // enter comment here
+    Float_t fEtRemovedPiMinus; // enter comment here
+    Float_t fEtRemovedKaonPlus; // enter comment here
+    Float_t fEtRemovedKaonMinus; // enter comment here
+    Float_t fEtRemovedK0s; // enter comment here
+    Float_t fEtRemovedK0L; // enter comment here
+    Float_t fEtRemovedLambdas; // enter comment here
+    Float_t fEtRemovedElectrons; // enter comment here
+    Float_t fEtRemovedPositrons; // enter comment here
+    Float_t fEtRemovedMuMinus; // enter comment here
+    Float_t fEtRemovedMuPlus; // enter comment here
+
+    Float_t fEtRemovedGammasFromPi0; // enter comment here
+    Float_t fEtRemovedGammas; // enter comment here
+    Float_t fEtRemovedNeutrons; // enter comment here
+    Float_t fEtRemovedAntiNeutrons; // enter comment here
+
+    Int_t fMultNonRemovedProtons; // enter comment here
+    Int_t fMultNonRemovedAntiProtons; // enter comment here
+    Int_t fMultNonRemovedPiPlus; // enter comment here
+    Int_t fMultNonRemovedPiMinus; // enter comment here
+    Int_t fMultNonRemovedKaonPlus; // enter comment here
+    Int_t fMultNonRemovedKaonMinus; // enter comment here
+    Int_t fMultNonRemovedK0s; // enter comment here
+    Int_t fMultNonRemovedK0L; // enter comment here
+    Int_t fMultNonRemovedLambdas; // enter comment here
+    Int_t fMultNonRemovedElectrons; // enter comment here
+    Int_t fMultNonRemovedPositrons; // enter comment here
+    Int_t fMultNonRemovedMuMinus; // enter comment here
+    Int_t fMultNonRemovedMuPlus; // enter comment here
+    Int_t fMultNonRemovedGammas; // enter comment here
+    Int_t fMultNonRemovedNeutrons; // enter comment here
+    Int_t fMultNonRemovedAntiNeutrons; // enter comment here
+
+    Int_t fMultRemovedProtons; // enter comment here
+    Int_t fMultRemovedAntiProtons; // enter comment here
+    Int_t fMultRemovedPiPlus; // enter comment here
+    Int_t fMultRemovedPiMinus; // enter comment here
+    Int_t fMultRemovedKaonPlus; // enter comment here
+    Int_t fMultRemovedKaonMinus; // enter comment here
+    Int_t fMultRemovedK0s; // enter comment here
+    Int_t fMultRemovedK0L; // enter comment here
+    
+    Int_t fMultRemovedLambdas; // enter comment here
+    Int_t fMultRemovedElectrons; // enter comment here
+    Int_t fMultRemovedPositrons; // enter comment here
+    Int_t fMultRemovedMuMinus; // enter comment here
+    Int_t fMultRemovedMuPlus; // enter comment here
+
+    Int_t fMultRemovedGammas; // enter comment here
+    Int_t fMultRemovedNeutrons; // enter comment here
+    Int_t fMultRemovedAntiNeutrons; // enter comment here
+
+    Int_t fTrackMultInAcc; // enter comment here
+
+
+    TH2F *fHistDxDzNonRemovedCharged; // enter comment here
+    TH2F *fHistDxDzRemovedCharged; // enter comment here
+    TH2F *fHistDxDzNonRemovedNeutral; // enter comment here
+    TH2F *fHistDxDzRemovedNeutral; // enter comment here
+
+    TH1F *fHistPiPlusMult; // enter comment here
+    TH1F *fHistPiMinusMult; // enter comment here
+    TH1F *fHistPiZeroMult; // enter comment here
+
+    TH1F *fHistPiPlusMultAcc; // enter comment here
+    TH1F *fHistPiMinusMultAcc; // enter comment here
+    TH1F *fHistPiZeroMultAcc; // enter comment here
+
+    Int_t fPiPlusMult; // enter comment here
+    Int_t fPiMinusMult; // enter comment here
+    Int_t fPiZeroMult; // enter comment here
+
+    Int_t fPiPlusMultAcc; // enter comment here
+    Int_t fPiMinusMultAcc; // enter comment here
+    Int_t fPiZeroMultAcc; // enter comment here
+
+
+    Int_t fNeutralRemoved; // number of neutral particles that where removed by track matching
+    Int_t fChargedRemoved; // number of charged particles that where removed by track matching
+    Int_t fChargedNotRemoved; // number of charged particles that were not removed
+    Int_t fNeutralNotRemoved; // number of neutral particles that were not removed
+
+    Double_t fEnergyNeutralRemoved; // energy of neutral particles that where removed by track matching
+    Double_t fEnergyChargedRemoved; // energy of charged particles that where removed by track matching
+    Double_t fEnergyChargedNotRemoved; // energy of charged particles that were not removed
+    Double_t fEnergyNeutralNotRemoved; // energy of neutral particles that were not removed
+
+    Int_t fNClusters; // Number of clusters in event
+
+    AliPHOSGeometry *fGeoUtils;
+
+    TH2I *fBadMapM2; // Bad map
+    TH2I *fBadMapM3; // Bad map
+    TH2I *fBadMapM4; // Bad map
+
+private:
 
     //Declare it private to avoid compilation warning
     AliAnalysisEtMonteCarlo & operator = (const AliAnalysisEtMonteCarlo & g) ;//cpy assignment
index e6cdc6e..660de03 100644 (file)
@@ -18,6 +18,7 @@ ClassImp(AliAnalysisEtMonteCarloEmcal);
 AliAnalysisEtMonteCarloEmcal::AliAnalysisEtMonteCarloEmcal()
 {
    fHistogramNameSuffix = TString("EmcalMC");
+
 }
 
 AliAnalysisEtMonteCarloEmcal::~AliAnalysisEtMonteCarloEmcal()
@@ -28,7 +29,7 @@ AliAnalysisEtMonteCarloEmcal::~AliAnalysisEtMonteCarloEmcal()
 void AliAnalysisEtMonteCarloEmcal::Init()
 { // Init
   AliAnalysisEtMonteCarlo::Init();
-  
+   fClusterType = fCuts->GetEmcalClusterType();
   fDetectorRadius = fCuts->GetGeometryEmcalDetectorRadius();
   fEtaCutAcc = fCuts->GetGeometryEmcalEtaAccCut();
   fPhiCutAccMax = fCuts->GetGeometryEmcalPhiAccMaxCut() * TMath::Pi()/180.;
index e9c5eac..91c29f2 100644 (file)
@@ -7,15 +7,24 @@
 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
 //_________________________________________________________________________
 #include "AliAnalysisEtMonteCarloPhos.h"
+#include "AliAnalysisEtSelectorPhos.h"
 #include "AliAnalysisEtCuts.h"
 #include "AliESDtrack.h"
+#include <iostream>
+#include "AliPHOSGeoUtils.h"
+#include "TFile.h"
+#include "TH2I.h"
+#include <AliPHOSGeometry.h>
 
 using namespace std;
 
 ClassImp(AliAnalysisEtMonteCarloPhos);
 
 
-AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos()
+AliAnalysisEtMonteCarloPhos::AliAnalysisEtMonteCarloPhos():AliAnalysisEtMonteCarlo()
+,fBadMapM2(0)
+,fBadMapM3(0)
+,fBadMapM4(0)
 {
    fHistogramNameSuffix = TString("PhosMC");
 }
@@ -28,7 +37,10 @@ AliAnalysisEtMonteCarloPhos::~AliAnalysisEtMonteCarloPhos()
 void AliAnalysisEtMonteCarloPhos::Init()
 { // 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.;
@@ -41,5 +53,119 @@ void AliAnalysisEtMonteCarloPhos::Init()
   fTrackDzCut = fCuts->GetPhosTrackDzCut();
     
   fDetector = fCuts->GetDetectorPhos();
-    
+  
+ // 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");
+// 
+//  fGeoUtils = new AliPHOSGeoUtils("PHOS", "noCPV");
+   fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
+  
 }
+
+
+Bool_t AliAnalysisEtMonteCarloPhos::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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << 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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << 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 << "d: " << distance << ", cut: " << fCuts->GetPhosBadDistanceCut() << std::endl;
+                    
+                        return kTRUE;
+                    }
+                }
+            }
+
+        }
+    }
+
+    return kFALSE;
+}
+
+
index 8d0177c..09881ec 100644 (file)
@@ -1,6 +1,7 @@
 #ifndef ALIANALYSISETMONTECARLOPHOS_H
 #define ALIANALYSISETMONTECARLOPHOS_H
-//_________________________________________________________________________
+
+class AliPHOSGeoUtils;//_________________________________________________________________________
 //  Utility Class for transverse energy studies
 //  Base class for MC analysis, for PHOS
 //  - MC output
@@ -9,7 +10,7 @@
 //_________________________________________________________________________
 
 #include "AliAnalysisEtMonteCarlo.h"
-
+class TH2I;
 
 class AliAnalysisEtMonteCarloPhos : public AliAnalysisEtMonteCarlo
 {
@@ -20,9 +21,18 @@ public:
     virtual ~AliAnalysisEtMonteCarloPhos();
 
     virtual void Init();
+protected:
+      virtual Bool_t TooCloseToBadChannel(const AliESDCaloCluster &cluster) const;
 
  private:
-
+   
+    TH2I *fBadMapM2; // Bad map
+    TH2I *fBadMapM3; // Bad map
+    TH2I *fBadMapM4; // Bad map
+    
+    // Prohibited
+    AliAnalysisEtMonteCarloPhos & operator = (const AliAnalysisEtMonteCarloPhos&) ;//cpy assignment
+    AliAnalysisEtMonteCarloPhos(const AliAnalysisEtMonteCarloPhos&) ; // cpy ctor
    ClassDef(AliAnalysisEtMonteCarloPhos, 1); 
 };
 
index d9e5f5b..3d73863 100644 (file)
 #include "AliESDpid.h"
 #include <iostream>
 #include "TH2F.h"
+#include "TH2I.h"
+#include "TH1I.h"
+#include "TFile.h"
 #include "AliAnalysisHadEtCorrections.h"
+#include "AliAnalysisEtSelector.h"
 #include "AliLog.h"
 #include "AliCentrality.h"
-
-
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSGeometry.h"
 
 
 using namespace std;
@@ -37,29 +41,37 @@ ClassImp(AliAnalysisEtReconstructed);
 
 
 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
-        AliAnalysisEt()
-        ,fCorrections(0)
-        ,fPidCut(0)
-        ,fClusterType(0)
-        ,fHistChargedPionEnergyDeposit(0)
-        ,fHistProtonEnergyDeposit(0)
-        ,fHistAntiProtonEnergyDeposit(0)
-        ,fHistChargedKaonEnergyDeposit(0)
-        ,fHistMuonEnergyDeposit(0)
-       ,fHistRemovedEnergy(0)
-        ,fGeomCorrection(1.0)
-        ,fEMinCorrection(1.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)
+    ,fGeoUtils(0)
+    ,fBadMapM2(0)
+    ,fBadMapM3(0)
+    ,fBadMapM4(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 */    
-    delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */    
-    delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */    
+    delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
+    delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
+    delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
+    delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
 
     delete fHistRemovedEnergy; // removed energy
@@ -81,155 +93,26 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         AliFatal("ERROR: ESD Event does not exist");
         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)
     {
         cent = fCentrality->GetCentralityClass10("V0M");
-       fCentClass = fCentrality->GetCentralityClass10("V0M");
+        fCentClass = fCentrality->GetCentralityClass10("V0M");
     }
 
-    Double_t protonMass = fgProtonMass;
-
-    //for PID
-    AliESDpid pID;
-    pID.MakePID(event);
-    TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event);
-
-    Int_t nGoodTracks = list->GetEntries();
-    // printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters());
-
-    for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++)
-    {
-        AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
-        if (!track)
-        {
-            AliError(Form("ERROR: Could not get track %d", iTrack));
-            continue;
-        }
-
-        fMultiplicity++;
-
-
-        Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron;
-        nSigmaPion = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kPion));
-        nSigmaProton = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kProton));
-        nSigmaKaon = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kKaon));
-        nSigmaElectron = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kElectron));
-        /*
-        bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0);
-        bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0);
-        bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0);
-        bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0);
-        */
-
-        Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
-        Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
-
-        Float_t massPart = 0;
-
-        const Double_t *pidWeights = track->PID();
-        Int_t maxpid = -1;
-        Double_t maxpidweight = 0;
-
-        if (pidWeights)
-        {
-            for (Int_t p =0; p < AliPID::kSPECIES; p++)
-            {
-                if (pidWeights[p] > maxpidweight)
-                {
-                    maxpidweight = pidWeights[p];
-                    maxpid = p;
-                }
-            }
-            if (maxpid == AliPID::kProton)
-            {
-                //by definition of ET
-                massPart = -protonMass*track->Charge();
-            }
-
-        }
-
-        Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
-       if(fMakeSparse){
-         fSparseTracks[0] = maxpid;
-         fSparseTracks[1] = track->Charge();
-         fSparseTracks[2] = track->M();
-         fSparseTracks[3] = et;
-         fSparseTracks[4] = track->Pt();
-         fSparseTracks[5] = track->Eta();
-         fSparseTracks[6] = cent;
-         fSparseHistTracks->Fill(fSparseTracks);
-       }
-        //printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout
-
-        if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
-        {
-            fTotChargedEt +=  et;
-            fChargedMultiplicity++;
-            if (maxpid != -1)
-            {
-                if (maxpid == AliPID::kProton)
-                {
-                    fProtonEt += et;
-                }
-                if (maxpid == AliPID::kPion)
-                {
-                    fPionEt += et;
-                }
-                if (maxpid == AliPID::kKaon)
-                {
-                    fChargedKaonEt += et;
-                }
-                if (maxpid == AliPID::kMuon)
-                {
-                    fMuonEt += et;
-                }
-                if (maxpid == AliPID::kElectron)
-                {
-                    fElectronEt += et;
-                }
-            }
-
-            if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
-            {
-                fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
-                if (maxpid != -1)
-                {
-                    if (maxpid == AliPID::kProton)
-                    {
-                        fProtonEtAcc += et;
-                    }
-                    if (maxpid == AliPID::kPion)
-                    {
-                        fPionEtAcc += et;
-                    }
-                    if (maxpid == AliPID::kKaon)
-                    {
-                        fChargedKaonEtAcc += et;
-                    }
-                    if (maxpid == AliPID::kMuon)
-                    {
-                        fMuonEtAcc += et;
-                    }
-                    if (maxpid == AliPID::kElectron)
-                    {
-                        fElectronEtAcc += et;
-                    }
-                }
-
-            }
-        }
-
-        if (TrackHitsCalorimeter(track, event->GetMagneticField()))
-        {
-            Double_t phi = track->Phi();
-            Double_t pt = track->Pt();
-            // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
-            if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
-            else fHistPhivsPtNeg->Fill(phi, pt);
-        }
-    }
+    //Double_t protonMass = fgProtonMass;
 
     for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
     {
@@ -239,13 +122,14 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
             AliError(Form("ERROR: Could not get cluster %d", iCluster));
             continue;
         }
-        if (cluster->GetType() != fClusterType) continue;
-
-        //if(cluster->GetTracksMatched() > 0)
-//     printf("Rec Cluster: iCluster %03d E %4.3f type %.qd NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
-
-
-        if (cluster->E() < fClusterEnergyCut) continue;
+        int x = 0;
+        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++);
 
         Float_t pos[3];
 
@@ -258,71 +142,15 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
             distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
         }
 
-       if(fMakeSparse){
-        fSparseClusters[0] = 0;
-        fSparseClusters[1] = 0;
-        fSparseClusters[2] = 0;
-        fSparseClusters[6] = 0;
-        fSparseClusters[7] = 0;
-        fSparseClusters[8] = 0;
-        fSparseClusters[9] = cent;
-        fSparseClusters[10] = 0;
-       }
-
-        if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex > -1)
-        {
-            AliVTrack *tmptrack = event->GetTrack(trackMatchedIndex);
-            if (!tmptrack)
-            {
-                AliError("Error: track does not exist");
-                return -1;
-            }
-            const Double_t *pidWeights = tmptrack->PID();
-
-            Double_t maxpidweight = 0;
-            Int_t maxpid = 0;
-            Double_t massPart = 0;
-            if (pidWeights)
-            {
-                for (Int_t p =0; p < AliPID::kSPECIES; p++)
-                {
-                    if (pidWeights[p] > maxpidweight)
-                    {
-                        maxpidweight = pidWeights[p];
-                        maxpid = p;
-                    }
-                }
-                if (maxpid == AliPID::kProton)
-                {
-                    //by definition of ET
-                    massPart = -protonMass*tmptrack->Charge();
-                }
-            }
-       if(fMakeSparse){
-            fSparseClusters[0] = maxpid;
-            fSparseClusters[1] = tmptrack->Charge();
-            fSparseClusters[2] = tmptrack->M();
-            fSparseClusters[6] = tmptrack->E() * TMath::Sin(tmptrack->Theta()) + massPart;;
-            fSparseClusters[7] = tmptrack->Pt();
-            fSparseClusters[8] = tmptrack->Eta();
-       }
-        }
-
-        
-       if(fMakeSparse){fSparseClusters[10] = distance;}
-
-        fHistTMDeltaR->Fill(distance);
-        fHistTMDxDz->Fill(cluster->GetTrackDx(), cluster->GetTrackDz());
-
-//        Float_t clusteret = cluster->E() * TMath::Sin(cp.Theta());
-
         Bool_t matched = false;
 
-        if (cluster->IsEMCAL()) matched = distance < fTrackDistanceCut;
-        else matched = (TMath::Abs(cluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(cluster->GetTrackDz()) < fTrackDzCut);
+        Double_t r = 0;
+
+        matched = !fSelector->CutTrackMatching(*cluster, r);
 
         if (matched)
         {
+
             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
             {
                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
@@ -405,58 +233,49 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
         } // distance
         else
         {
-       if(fMakeSparse){
-         fSparseClusters[0] = AliPID::kPhoton;
-         fSparseClusters[1] = 0;
-       }
+            fCutFlow->Fill(x++);
+            //std::cout << x++ << std::endl;
+            fSparseClusters[0] = AliPID::kPhoton;
+            fSparseClusters[1] = 0;
 
-            if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
-            if (cluster->E() < fClusterEnergyCut) continue;
+            //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
+            //if (cluster->E() < fClusterEnergyCut) continue;
             cluster->GetPosition(pos);
 
-            // TODO: replace with TVector3, too lazy now...
+            TVector3 p2(pos);
 
-            float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+            fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity());
 
-            float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
-            // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
-            fTotNeutralEt += cluster->E() * TMath::Sin(theta);
+            fTotNeutralEt += CalculateTransverseEnergy(cluster);
             fNeutralMultiplicity++;
         }
+        fMultiplicity++;
+    }
 
-        cluster->GetPosition(pos);
-
-        // TODO: replace with TVector3
+    fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
+    fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
+    fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
+    fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
 
-        float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
-        float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
-        float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
-       if(fMakeSparse){
-         fSparseClusters[3] = cluster->E() * TMath::Sin(theta);
-         fSparseClusters[4] = cluster->E();
-         fSparseClusters[5] = eta;
-         
-         fSparseHistClusters->Fill(fSparseClusters);
-       }
+    fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
+    fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
 
-        fMultiplicity++;
-    }
     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity);
     fHistRemovedEnergy->Fill(removedEnergy);
-//  std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
-    fTotNeutralEt = fGeomCorrection * fEMinCorrection * fTotNeutralEt - removedEnergy;
-    fTotNeutralEtAcc = fTotNeutralEt/fGeomCorrection;
+
+    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;
-       }
+    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...
     FillHistograms();
 
@@ -464,7 +283,7 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
 }
 
 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
-{ // check vertex
+{   // check vertex
 
     Float_t bxy = 999.;
     Float_t bz = 999.;
@@ -490,17 +309,25 @@ 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;
@@ -523,7 +350,7 @@ bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Doubl
 }
 
 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);
@@ -533,10 +360,15 @@ void AliAnalysisEtReconstructed::FillOutputList(TList* list)
     list->Add(fHistMuonEnergyDeposit);
 
     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;
@@ -581,6 +413,20 @@ void AliAnalysisEtReconstructed::CreateHistograms()
     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
 
+    histname = "fClusterPosition" + fHistogramNameSuffix;
+    fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13);
+    fClusterPosition->SetXTitle("Energy deposited in calorimeter");
+    fClusterPosition->SetYTitle("Energy of track");
+
+
+    histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
+    fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
+
+    histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
+    fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
+
+    histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
+    fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
 
 
 }
@@ -589,7 +435,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};
 
@@ -632,3 +478,169 @@ AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3],
     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 584290a..718487d 100644 (file)
@@ -1,6 +1,10 @@
 #ifndef ALIANALYSISETRECONSTRUCTED_H
 #define ALIANALYSISETRECONSTRUCTED_H
-//_________________________________________________________________________
+
+class TH2F;
+class TH2D;
+class AliPHOSGeometry;
+class AliPHOSGeoUtils;//_________________________________________________________________________
 //  Utility Class for transverse energy studies
 //  Base class for ESD analysis
 //  - reconstruction output
@@ -10,6 +14,7 @@
 
 #include "AliAnalysisEt.h"
 
+class TH2I;
 class AliVParticle;
 class AliESDEvent;
 class AliAnalysisHadEtCorrections;
@@ -18,7 +23,7 @@ class AliAnalysisEtReconstructed : public AliAnalysisEt
 {
 
 public:
-   
+
     AliAnalysisEtReconstructed();
     virtual ~AliAnalysisEtReconstructed();
 
@@ -28,39 +33,61 @@ public:
 
     /** Fill the objects you want to output, classes which add new histograms should overload this. */
     virtual void FillOutputList(TList *list);
-    void SetCorrections(AliAnalysisHadEtCorrections *corr){fCorrections = corr;}
+    void SetCorrections(AliAnalysisHadEtCorrections *corr) {
+        fCorrections = corr;
+    }
 
     /** Create the histograms, must be overloaded if you want to add your own */
     virtual void CreateHistograms();
     
+    void SetEMinCorrection(const Double_t factor) { fEMinCorrection = factor; }
+
 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
 
     Double_t fPidCut; // cut on the pid probability
-    
-    Char_t fClusterType; // selection on cluster type
-        
-    TH2F *fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */    
-    TH2F *fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */    
-    TH2F *fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */    
-    TH2F *fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */    
+
+    TH2F *fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
+    TH2F *fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
+    TH2F *fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
+    TH2F *fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
     TH2F *fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
 
     TH1F *fHistRemovedEnergy; // removed energy
-
+    
     Double_t fGeomCorrection; // geometry correction
     Double_t fEMinCorrection; // Emin correction
     
-  private:
+    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;
+    TH2D *fHistNeutralEnergyRemoved;
+    TH2D *fHistGammaEnergyAdded;
+    
+
+private:
+
     AliAnalysisEtReconstructed(const AliAnalysisEtReconstructed& g);
     AliAnalysisEtReconstructed & operator=(const AliAnalysisEtReconstructed&);
-    
-    Double_t CalcTrackClusterDistance(const Float_t pos[3],Int_t *trkMatchId, const AliESDEvent *event);
+
+
 
     ClassDef(AliAnalysisEtReconstructed, 1);
 };
index 8c0ec12..4814a70 100644 (file)
@@ -37,7 +37,7 @@ void AliAnalysisEtReconstructedEmcal::Init()
   fClusterEnergyCut = fCuts->GetReconstructedEmcalClusterEnergyCut();
   fSingleCellEnergyCut = fCuts->GetReconstructedEmcalSingleCellEnergyCut();
 
-  fClusterType = fCuts->GetReconstructedEmcalClusterType();
+  fClusterType = fCuts->GetEmcalClusterType();
   fTrackDistanceCut = fCuts->GetEmcalTrackDistanceCut();
   
   fDetector = fCuts->GetDetectorEmcal();
index c2b2c7b..be2bb69 100644 (file)
@@ -7,8 +7,16 @@
 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
 //_________________________________________________________________________
 #include "AliAnalysisEtReconstructedPhos.h"
+#include "AliAnalysisEtSelectorPhos.h"
 #include "AliAnalysisEtCuts.h"
 #include "AliESDtrack.h"
+#include <fstream>
+#include <iostream>
+#include <AliPHOSGeoUtils.h>
+#include "AliPHOSGeometry.h"
+#include <vector>
+#include "TH2I.h"
+
 using namespace std;
 
 ClassImp(AliAnalysisEtReconstructedPhos);
@@ -18,97 +26,109 @@ const Double_t kMEANCHARGED = 0.335;
 const Double_t kMEANNEUTRAL = 0.434;
 const Double_t kMEANGAMMA = 0.374;
 
-// Best case (pions and K0s): 
+// Best case (pions and K0s):
 const Double_t kMEANCHARGED = 0.304;
 const Double_t kMEANNEUTRAL = 0.3356;
 const Double_t kMEANGAMMA = 0.374;
 */
 // Simulated case:
-const Double_t kMEANCHARGED = 0.307;
-const Double_t kMEANNEUTRAL = 0.407;
-const Double_t kMEANGAMMA = 0.374;
+//corrEnergy =cluster->E()/(0.51 + 0.02*cluster->E());
+const Double_t kMEANCHARGED = 0.307/(0.51 + 0.02*0.307);
+const Double_t kMEANNEUTRAL = 0.407/(0.51 + 0.02*0.407);
+const Double_t kMEANGAMMA = 0.374/(0.51 + 0.02*0.374);
 
 
 AliAnalysisEtReconstructedPhos::AliAnalysisEtReconstructedPhos() :
-AliAnalysisEtReconstructed()
+        AliAnalysisEtReconstructed()
 {
-   fHistogramNameSuffix = TString("PhosRec");    
+    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;
+
 }
 
-AliAnalysisEtReconstructedPhos::~AliAnalysisEtReconstructedPhos() 
+AliAnalysisEtReconstructedPhos::~AliAnalysisEtReconstructedPhos()
 {
 }
 
 void AliAnalysisEtReconstructedPhos::Init()
 { // Init
-  AliAnalysisEtReconstructed::Init();
+    AliAnalysisEtReconstructed::Init();
+
+    fSelector = new AliAnalysisEtSelectorPhos(fCuts);
     
-  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->GetReconstructedPhosClusterType();
-  fTrackDistanceCut = fCuts->GetPhosTrackDistanceCut();
-  fTrackDxCut = fCuts->GetPhosTrackDxCut();
-  fTrackDzCut = fCuts->GetPhosTrackDzCut();
-  
-  fDetector = fCuts->GetDetectorPhos();
-  
-  fGeomCorrection = 1.0/0.036;
+    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();
   
-  fEMinCorrection = 1.0;
 
 }
 
 bool AliAnalysisEtReconstructedPhos::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
 {
-  return  AliAnalysisEtReconstructed::TrackHitsCalorimeter(track, magField);
+    return  AliAnalysisEtReconstructed::TrackHitsCalorimeter(track, magField);
 }
 
 Double_t AliAnalysisEtReconstructedPhos::GetChargedContribution(Int_t clusterMult)
 { // Charged contrib
-  if(clusterMult > 0)
-  {
-    Double_t nPart = 0.067 + 0.137*clusterMult;
-  
-    Double_t contr = nPart*kMEANCHARGED;
-  
-    return contr;
-  }
-  return 0;
-  
+    if (clusterMult > 0)
+    {
+        Double_t nPart = fChargedContributionCorrectionParameters[0] + fChargedContributionCorrectionParameters[1]*clusterMult;
+
+        Double_t contr = nPart*kMEANCHARGED;
+
+        return contr;
+    }
+    return 0;
+
 }
 
 Double_t AliAnalysisEtReconstructedPhos::GetNeutralContribution(Int_t clusterMult)
 { // Neutral contrib
-  if(clusterMult > 0)
-  {
-    Double_t nPart = 0.012 + 0.024*clusterMult - 0.00006*clusterMult*clusterMult;
-  
-    Double_t contr = nPart*kMEANNEUTRAL;
-  
-    return contr;
-  }
-  return 0;
+    if (clusterMult > 0)
+    {
+        Double_t nPart = fNeutralContributionCorrectionParameters[0] + fNeutralContributionCorrectionParameters[1]*clusterMult + fNeutralContributionCorrectionParameters[2]*clusterMult*clusterMult;
+
+        Double_t contr = nPart*kMEANNEUTRAL;
+
+        return contr;
+    }
+    return 0;
 }
 
 Double_t AliAnalysisEtReconstructedPhos::GetGammaContribution(Int_t clusterMult)
 { // Gamma contrib
-  if(clusterMult > 0)
-  {
-    Double_t nPart = -0.008 + 0.0057*clusterMult + 0.0002*clusterMult*clusterMult;
-  
-    Double_t contr = nPart*kMEANGAMMA;
-  
-    return contr;
-  }
-  return 0;
-}
+    if (clusterMult > 0)
+    {
+        Double_t nPart = fRemovedGammaContributionCorrectionParameters[0] + fRemovedGammaContributionCorrectionParameters[1]*clusterMult + fRemovedGammaContributionCorrectionParameters[2]*clusterMult*clusterMult;
 
+        Double_t contr = nPart*kMEANGAMMA;
 
+        return contr;
+    }
+    return 0;
+}
 
 
 
index 96e1cc0..0365e12 100644 (file)
@@ -16,24 +16,48 @@ class AliAnalysisEtReconstructedPhos : public AliAnalysisEtReconstructed
 {
 
 public:
-   
+
     AliAnalysisEtReconstructedPhos();
     virtual ~AliAnalysisEtReconstructedPhos();
 
     virtual void Init();
-    
+
     virtual Double_t GetChargedContribution(Int_t clusterMultiplicity);
-    
+
     virtual Double_t GetNeutralContribution(Int_t clusterMultiplicity);
-    
+
     virtual Double_t GetGammaContribution(Int_t clusterMultiplicity);
+
+    void SetChargedContributionParameters(Double_t par[2])
+    {
+        fChargedContributionCorrectionParameters[0] = par[0];
+        fChargedContributionCorrectionParameters[1] = par[1];
+    }
+    void SetNeutralContributionParameters(Double_t par[3])
+    {
+        fNeutralContributionCorrectionParameters[0] = par[0];
+        fNeutralContributionCorrectionParameters[1] = par[1];
+       fNeutralContributionCorrectionParameters[2] = par[2];
+    }
+    void SetRemovedGammaContributionParameters(Double_t par[3])
+    {
+        fRemovedGammaContributionCorrectionParameters[0] = par[0];
+        fRemovedGammaContributionCorrectionParameters[1] = par[1];
+       fRemovedGammaContributionCorrectionParameters[2] = par[2];
+    }
+
+protected:
+
+    virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
     
-   protected:
-      
-      virtual bool TrackHitsCalorimeter(AliVParticle *track, Double_t magField);
- private:
+private:
+
+    Double_t fChargedContributionCorrectionParameters[2]; // 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
+
 
-      ClassDef(AliAnalysisEtReconstructedPhos, 1);
+    ClassDef(AliAnalysisEtReconstructedPhos, 1);
 };
 
 #endif // ALIANALYSISETRECONSTRUCTEDPHOS_H
diff --git a/PWGLF/totEt/AliAnalysisEtSelector.cxx b/PWGLF/totEt/AliAnalysisEtSelector.cxx
new file mode 100644 (file)
index 0000000..d31abaa
--- /dev/null
@@ -0,0 +1,19 @@
+#include "AliAnalysisEtSelector.h"
+#include "AliAnalysisEtCuts.h"
+#include <AliESDCaloCluster.h>
+
+ClassImp(AliAnalysisEtSelector)
+
+AliAnalysisEtSelector::AliAnalysisEtSelector(AliAnalysisEtCuts *cuts) :
+fEvent(0)
+,fCuts(cuts)
+,fRunNumber(0)
+{
+
+}
+
+AliAnalysisEtSelector::~AliAnalysisEtSelector()
+{
+
+}
+
diff --git a/PWGLF/totEt/AliAnalysisEtSelector.h b/PWGLF/totEt/AliAnalysisEtSelector.h
new file mode 100644 (file)
index 0000000..b44660c
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALIANALYSISETSELECTOR_H
+#define ALIANALYSISETSELECTOR_H
+#include <Rtypes.h>
+
+class AliVEvent;
+class AliESDCaloCluster;
+class TRefArray;
+class AliAnalysisEtCuts;
+class AliAnalysisEtSelector
+{
+
+public:
+  
+    // Constructor takes cuts object
+    AliAnalysisEtSelector(AliAnalysisEtCuts *cuts);
+    
+    // Destructor 
+    virtual ~AliAnalysisEtSelector();
+    
+    // Set the current event
+    void SetEvent(const AliVEvent *event) { fEvent = event; }
+    
+    // Init 
+    virtual Int_t Init(Int_t runNumber) { fRunNumber = runNumber; return 0; }
+    
+    // Return CaloClusters for the detector
+    virtual TRefArray* GetClusters() = 0;
+    
+    // Return true if cluster has energy > cut
+    virtual Bool_t CutMinEnergy(const AliESDCaloCluster & /*cluster*/) const { return true; }
+    
+    // Cut on distance to bad channel
+    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; }
+  
+protected:
+  
+    const AliVEvent *fEvent; // Pointer to current event
+
+    AliAnalysisEtCuts *fCuts; // Pointer to the cuts object
+    
+    Int_t fRunNumber;
+    
+private:
+
+    AliAnalysisEtSelector(); // Prohibited
+    AliAnalysisEtSelector(const AliAnalysisEtSelector& other);// Prohibited
+    AliAnalysisEtSelector& operator=(const AliAnalysisEtSelector& other);// Prohibited
+    bool operator==(const AliAnalysisEtSelector& other) const;// Prohibited
+    
+    ClassDef(AliAnalysisEtSelector, 1);
+};
+
+#endif // ALIANALYSISETSELECTOR_H
diff --git a/PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx b/PWGLF/totEt/AliAnalysisEtSelectorPhos.cxx
new file mode 100644 (file)
index 0000000..085dc91
--- /dev/null
@@ -0,0 +1,234 @@
+#include "AliAnalysisEtSelectorPhos.h"
+#include "AliAnalysisEtCuts.h"
+#include "AliESDCaloCluster.h"
+#include <AliVEvent.h>
+#include "TRefArray.h"
+#include "AliPHOSGeometry.h"
+#include "TH2I.h"
+#include "TFile.h"
+#include "TMath.h"
+#include "AliLog.h"
+#include <iostream>
+ClassImp(AliAnalysisEtSelectorPhos)
+AliAnalysisEtSelectorPhos::AliAnalysisEtSelectorPhos(AliAnalysisEtCuts* cuts): AliAnalysisEtSelector(cuts)
+,fGeoUtils(0)
+,fBadMapM2(0)
+,fBadMapM3(0)
+,fBadMapM4(0)
+{
+  
+}
+
+AliAnalysisEtSelectorPhos::~AliAnalysisEtSelectorPhos()
+{
+
+}
+
+TRefArray* AliAnalysisEtSelectorPhos::GetClusters()
+{  
+  return 0;
+}
+
+Int_t AliAnalysisEtSelectorPhos::Init(Int_t runNumber)
+{
+  AliAnalysisEtSelector::Init(runNumber);
+  
+  int res = LoadGeometry();
+  if(res) return -1;
+  return LoadBadMaps();  
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutMinEnergy(const AliESDCaloCluster& cluster) const
+{
+  return cluster.E() > fCuts->GetReconstructedPhosClusterEnergyCut();
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutDistanceToBadChannel(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 kFALSE;
+                    }
+                }
+            }
+            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 kFALSE;
+                    }
+                }
+            }
+            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 kFALSE;
+                    }
+                }
+            }
+
+        }
+    }
+
+    return kTRUE;
+
+}
+
+Bool_t AliAnalysisEtSelectorPhos::CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const
+{
+
+  
+  // cluster->GetTrackDx(), cluster->GetTrackDz(), event->GetTrack(trackMatchedIndex)->Pt(), event->GetTrack(trackMatchedIndex)->Charge(), ev
+  
+  Int_t nTracksMatched = cluster.GetNTracksMatched();
+  if(nTracksMatched == 0) return kFALSE;
+  
+  Int_t trackMatchedIndex = cluster.GetTrackMatchedIndex();
+  if(trackMatchedIndex < 0) return kTRUE;
+  
+  AliVParticle *track = fEvent->GetTrack(trackMatchedIndex);
+  Double_t dx = cluster.GetTrackDx();
+  Double_t dz = cluster.GetTrackDz();
+  Double_t pt = track->Pt();
+  Int_t charge = track->Charge();
+  
+  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 = fEvent->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 ;
+  r = TMath::Sqrt(rx*rx+rz*rz);
+  if(r < fCuts->GetPhosTrackRCut()) return kFALSE;
+  
+  return kTRUE;
+}
+
+int AliAnalysisEtSelectorPhos::LoadGeometry()
+{
+
+  fGeoUtils = AliPHOSGeometry::GetInstance("IHEP");
+    // ifstream f("badchannels.txt", ios::in);
+  return 0;
+}
+
+int AliAnalysisEtSelectorPhos::LoadBadMaps()
+{
+TFile *f = TFile::Open("badchannels.root", "READ");
+
+    if(!f)
+    {
+      std::cout << "Could not open badchannels.root" << std::endl;
+      return -1;
+    }
+    
+    fBadMapM2 = (TH2I*)f->Get("bad_channels_m2");
+    if(fBadMapM2) 
+    {
+      std::cout << "Could not find bad_channels_m2 in badchannels.root" << std::endl;
+    }
+    fBadMapM3 = (TH2I*)f->Get("bad_channels_m3");
+    if(fBadMapM2) 
+    {
+      std::cout << "Could not find bad_channels_m3 in badchannels.root" << std::endl;
+    }
+    
+    fBadMapM4 = (TH2I*)f->Get("bad_channels_m4");
+    if(fBadMapM4) 
+    {
+      std::cout << "Could not find bad_channels_m4 in badchannels.root" << std::endl;
+    }
+
+    return 0;
+    
+    
+    
+}
diff --git a/PWGLF/totEt/AliAnalysisEtSelectorPhos.h b/PWGLF/totEt/AliAnalysisEtSelectorPhos.h
new file mode 100644 (file)
index 0000000..67921b3
--- /dev/null
@@ -0,0 +1,44 @@
+#ifndef ALIANALYSISETSELECTORPHOS_H
+#define ALIANALYSISETSELECTORPHOS_H
+
+#include "AliAnalysisEtSelector.h"
+
+class TH2I;
+class AliPHOSGeometry;
+
+class AliAnalysisEtSelectorPhos : public AliAnalysisEtSelector
+{
+
+public:
+
+    AliAnalysisEtSelectorPhos(AliAnalysisEtCuts *cuts);
+    virtual ~AliAnalysisEtSelectorPhos();
+    
+    virtual TRefArray* GetClusters();
+    virtual Bool_t CutMinEnergy(const AliESDCaloCluster& cluster) const;
+    virtual Bool_t CutDistanceToBadChannel(const AliESDCaloCluster& cluster) const;
+    virtual Bool_t CutTrackMatching(const AliESDCaloCluster& cluster, Double_t &r) const;
+    
+    virtual Int_t Init(int runNumber);
+  
+private:
+
+
+    int LoadGeometry();
+    int LoadBadMaps();
+    
+    AliPHOSGeometry *fGeoUtils;
+    
+    TH2I *fBadMapM2; // Bad map
+    TH2I *fBadMapM3; // Bad map
+    TH2I *fBadMapM4; // Bad map
+    
+    AliAnalysisEtSelectorPhos();
+    AliAnalysisEtSelectorPhos(const AliAnalysisEtSelectorPhos& other); // Prohibited
+    AliAnalysisEtSelectorPhos& operator=(const AliAnalysisEtSelectorPhos& other); // Prohibited
+    bool operator==(const AliAnalysisEtSelectorPhos& other) const; // Prohibited
+    
+    ClassDef(AliAnalysisEtSelectorPhos, 1);
+};
+
+#endif // ALIANALYSISETSELECTORPHOS_H
index c9dd707..4c1f4e4 100644 (file)
@@ -33,6 +33,7 @@ AliAnalysisTaskTransverseEnergy::AliAnalysisTaskTransverseEnergy(const char* nam
         ,fEsdtrackCutsTPC(0)
         ,fEsdtrackCutsITS(0)
         ,fOutputList(0)
+       ,fPhysSelTaskName("physSelTask")
         ,fCentSelTaskName("centralityTask")
        ,fIsMc(isMc)
        ,fCurrentRunNum(-1)