- Adding sparse histograms to analysis
authorodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Apr 2011 18:08:33 +0000 (18:08 +0000)
committerodjuvsla <odjuvsla@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Apr 2011 18:08:33 +0000 (18:08 +0000)
- Getting reconstructed clusters based on detector cut
- Adding new histograms to check hadronic deposit

PWG4/totEt/AliAnalysisEtMonteCarlo.cxx
PWG4/totEt/AliAnalysisEtMonteCarlo.h

index 8b2fe68..50e57df 100644 (file)
@@ -5,7 +5,7 @@
 //  implementation file
 //
 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
-              //_________________________________________________________________________
+//_________________________________________________________________________
 
 #include "AliAnalysisEtMonteCarlo.h"
 #include "AliAnalysisEtCuts.h"
 #include "AliMCEvent.h"
 #include "AliESDEvent.h"
 #include "TH2F.h"
+#include "TH3F.h"
 #include "TParticle.h"
 #include "AliGenHijingEventHeader.h"
 #include "AliGenPythiaEventHeader.h"
 #include "TList.h"
 #include "AliESDCaloCluster.h"
 #include "AliLog.h"
+#include <iostream>
+#include <AliCentrality.h>
 
-  using namespace std;
+using namespace std;
 
 ClassImp(AliAnalysisEtMonteCarlo);
 
 
 // ctor
 AliAnalysisEtMonteCarlo::AliAnalysisEtMonteCarlo():AliAnalysisEt()
-                                                 ,fImpactParameter(0)
-                                                 ,fNcoll(0)
-                                                 ,fNpart(0)
-                                                 ,fHistPrimElectronEtaEET(0)  
-                                                 ,fHistSecElectronEtaEET(0) 
-                                                 ,fHistConvElectronEtaEET(0)  
-                                                 ,fHistPrimGammaEtaEET(0)  
-                                                 ,fHistPion0GammaEtaEET(0)  
-                                                 ,fHistEtaGammaEtaEET(0)  
-                                                 ,fHistOmega0GammaEtaEET(0)  
-                                                 ,fHistSecGammaEtaEET(0)  
-
-                                                 ,fHistPrimElectronEtaE(0) 
-                                                 ,fHistSecElectronEtaE(0)
-                                                 ,fHistConvElectronEtaE(0)
-                                                 ,fHistPrimGammaEtaE(0) 
-                                                 ,fHistPion0GammaEtaE(0)
-                                                 ,fHistEtaGammaEtaE(0)
-                                                 ,fHistOmega0GammaEtaE(0)
-                                                 ,fHistSecGammaEtaE(0) 
-
-                                                 ,fHistPrimElectronEtaERec(0) 
-                                                 ,fHistSecElectronEtaERec(0)
-                                                 ,fHistConvElectronEtaERec(0)
-                                                 ,fHistPrimGammaEtaERec(0) 
-                                                 ,fHistSecGammaEtaERec(0) 
-                                                 ,fHistPion0GammaEtaERec(0)
-                                                 ,fHistEtaGammaEtaERec(0)
-                                                 ,fHistOmega0GammaEtaERec(0)
-
-                                                 ,fHistAllERecEMC(0)    
-                                                 ,fHistGammaERecEMC(0) 
-                                                 ,fHistElectronERecEMC(0)      
-
-                                                 ,fHistElectronFirstMother(0)  
-                                                 ,fHistElectronLastMother(0)  
-                                                 ,fHistElectronFirstMotherEtaAcc(0)  
-                                                 ,fHistElectronFirstMotherNPP(0)  
-                                                 ,fHistElectronFirstMotherNPPAcc(0)  
-
-                                                 ,fHistGammaFirstMother(0)  
-                                                 ,fHistGammaLastMother(0)
-                                                 ,fHistGammaFirstMotherEtaAcc(0)
-                                                 ,fHistGammaFirstMotherNPP(0)
-                                                 ,fHistGammaFirstMotherNPPAcc(0)
+        ,fImpactParameter(0)
+        ,fNcoll(0)
+        ,fNpart(0)
+        ,fHistPrimElectronEtaEET(0)
+        ,fHistSecElectronEtaEET(0)
+        ,fHistConvElectronEtaEET(0)
+        ,fHistPrimGammaEtaEET(0)
+        ,fHistPion0GammaEtaEET(0)
+        ,fHistEtaGammaEtaEET(0)
+        ,fHistOmega0GammaEtaEET(0)
+        ,fHistSecGammaEtaEET(0)
+
+        ,fHistPrimElectronEtaE(0)
+        ,fHistSecElectronEtaE(0)
+        ,fHistConvElectronEtaE(0)
+        ,fHistPrimGammaEtaE(0)
+        ,fHistPion0GammaEtaE(0)
+        ,fHistEtaGammaEtaE(0)
+        ,fHistOmega0GammaEtaE(0)
+        ,fHistSecGammaEtaE(0)
+
+        ,fHistPrimElectronEtaERec(0)
+        ,fHistSecElectronEtaERec(0)
+        ,fHistConvElectronEtaERec(0)
+        ,fHistPrimGammaEtaERec(0)
+        ,fHistSecGammaEtaERec(0)
+        ,fHistPion0GammaEtaERec(0)
+        ,fHistEtaGammaEtaERec(0)
+        ,fHistOmega0GammaEtaERec(0)
+
+        ,fHistAllERecEMC(0)
+        ,fHistGammaERecEMC(0)
+        ,fHistElectronERecEMC(0)
+
+        ,fHistElectronFirstMother(0)
+        ,fHistElectronLastMother(0)
+        ,fHistElectronFirstMotherEtaAcc(0)
+        ,fHistElectronFirstMotherNPP(0)
+        ,fHistElectronFirstMotherNPPAcc(0)
+
+        ,fHistGammaFirstMother(0)
+        ,fHistGammaLastMother(0)
+        ,fHistGammaFirstMotherEtaAcc(0)
+        ,fHistGammaFirstMotherNPP(0)
+        ,fHistGammaFirstMotherNPPAcc(0)
+
+        ,fHistDecayVertexNonRemovedCharged(0)
+        ,fHistDecayVertexRemovedCharged(0)
+        ,fHistDecayVertexNonRemovedNeutral(0)
+        ,fHistDecayVertexRemovedNeutral(0)
+
+        ,fHistRemovedOrNot(0)
+        ,fHistDxDzRemovedCharged(0)
+        ,fHistDxDzNonRemovedCharged(0)
+        ,fHistDxDzRemovedNeutral(0)
+        ,fHistDxDzNonRemovedNeutral(0)
+        ,fNeutralRemoved(0)
+        ,fChargedRemoved(0)
+        ,fChargedNotRemoved(0)
+        ,fNeutralNotRemoved(0)
+        ,fEnergyNeutralRemoved(0)
+        ,fEnergyChargedRemoved(0)
+        ,fEnergyChargedNotRemoved(0)
+        ,fEnergyNeutralNotRemoved(0)
+
 {
+    fTrackDistanceCut = 7.0;
 }
 
 // dtor
-AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo() 
+AliAnalysisEtMonteCarlo::~AliAnalysisEtMonteCarlo()
 {
 }
 
 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev)
 { // analyse MC event
-  ResetEventValues();
-       
-  // Get us an mc event
-  if(!ev){  
-    AliFatal("ERROR: Event does not exist");
-    return 0;
-  }
-  AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
-  if(!event){  
-    AliFatal("ERROR: MC Event does not exist");
-    return 0;
-  }
-       
-  Double_t protonMass =fgProtonMass;
-       
-  // Hijing header
-  AliGenEventHeader* genHeader = event->GenEventHeader();
-  if(!genHeader){
-    Printf("ERROR: Event generation header does not exist");   
-    return 0;
-  }
-  AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
-  if (hijingGenHeader) {
-    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());
+    ResetEventValues();
+
+    fPiPlusMult = 0;
+    fPiMinusMult = 0;
+    fPiZeroMult = 0;
+    if (fCentrality)
+    {
+        fCentClass = fCentrality->GetCentralityClass10(fCentralityMethod);
+
+    }
+    fSparseTracks[0] = 0;
+    fSparseTracks[1] = 0;
+    fSparseTracks[2] = 0;
+    fSparseTracks[3] = 0;
+    fSparseTracks[4] = 0;
+    fSparseTracks[5] = 0;
+    fSparseTracks[6] = fCentClass;
+
+
+    // Get us an mc event
+    if (!ev) {
+        AliFatal("ERROR: Event does not exist");
+        return 0;
+    }
+    AliMCEvent *event = dynamic_cast<AliMCEvent*>(ev);
+    if (!event) {
+        AliFatal("ERROR: MC Event does not exist");
+        return 0;
+    }
+
+    Double_t protonMass =fgProtonMass;
+
+    // Hijing header
+    AliGenEventHeader* genHeader = event->GenEventHeader();
+    if (!genHeader) {
+        Printf("ERROR: Event generation header does not exist");
+        return 0;
+    }
+    AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
+    if (hijingGenHeader) {
+        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);
+       if (pythiaGenHeader) { // not Hijing; try with Pythia
+       printf("Pythia: ProcessType %d  GetPtHard %g \n",
+       pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
+       }
     */
-  }
-       
-  /* // placeholder if we want to get some Pythia info later
-     AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
-     if (pythiaGenHeader) { // not Hijing; try with Pythia      
-     printf("Pythia: ProcessType %d  GetPtHard %g \n",
-     pythiaGenHeader->ProcessType(), pythiaGenHeader->GetPtHard());
-     }
-  */
-       
-  // Let's play with the stack!
-  AliStack *stack = event->Stack();
-       
-  Int_t nPrim = stack->GetNtrack();
-       
-  for (Int_t iPart = 0; iPart < nPrim; iPart++)
+
+    // Let's play with the stack!
+    AliStack *stack = event->Stack();
+
+    Int_t nPrim = stack->GetNtrack();
+
+    Int_t partCount = 0;
+    for (Int_t iPart = 0; iPart < nPrim; iPart++)
     {
-               
-      TParticle *part = stack->Particle(iPart);
-      TParticle *partMom = 0;
-      TParticle *partMomLast = 0;
-               
-      if (!part)
+
+        TParticle *part = stack->Particle(iPart);
+        TParticle *partMom = 0;
+        TParticle *partMomLast = 0;
+
+        if (!part)
+        {
+            Printf("ERROR: Could not get particle %d", iPart);
+            continue;
+        }
+
+        Int_t iPartMom = part->GetMother(0);
+        Int_t iPartLastMom = part->GetMother(1);
+
+        TParticlePDG *pdg = part->GetPDG(0);
+        TParticlePDG *pdgMom = 0;
+        TParticlePDG *pdgMomLast = 0;
+
+        if (!pdg)
         {
-         Printf("ERROR: Could not get particle %d", iPart);
-         continue;
+            //Printf("ERROR: Could not get particle PDG %d", iPart);
+            continue;
         }
-               
-      Int_t iPartMom = part->GetMother(0);
-      Int_t iPartLastMom = part->GetMother(1);
 
-      TParticlePDG *pdg = part->GetPDG(0);
-      TParticlePDG *pdgMom = 0;
-      TParticlePDG *pdgMomLast = 0;
+        if (iPartMom>0)
+        {
+            partMom = stack->Particle(iPartMom);
+            pdgMom = partMom->GetPDG(0);
+        }
 
-      if (!pdg)
+        if (iPartLastMom>0)
         {
-         //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;
-               
-      // it goes to the next particle in case it is not a physical primary or not a electron or gamma (want to check the contribution from secondary)
-      if (!stack->IsPhysicalPrimary(iPart))
-       {                                       
-         // this part is used only to check the origin of secondary electrons and gammas 
-         if (iPartMom>0)
-           {
-             //if (!stack->IsPhysicalPrimary(iPartMom)) continue;
-                               
-             if (pdgMom)
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   {
-                     fHistElectronFirstMotherNPP->Fill(pdgMom->PdgCode());
-                     // inside EMCal acceptance
-                     if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
-                       fHistElectronFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
-                   }
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   {           
-                     fHistGammaFirstMotherNPP->Fill(pdgMom->PdgCode());
-                     if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
-                       fHistGammaFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
-                   }
-               }
-             else 
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   fHistElectronFirstMotherNPP->Fill(-598);
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   fHistGammaFirstMotherNPP->Fill(-598);
-                                       
-                 continue;
-               }
-           }
-         else
-           {
-             if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-               fHistElectronFirstMotherNPP->Fill(-599);
-             else if (pdg->PdgCode() == fgGammaCode)
-               fHistGammaFirstMotherNPP->Fill(-599);
-                               
-             continue;
-           }
-                       
-         if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
-           continue;
-       }
-      // this part is used only to check the origin of physical primary electrons and gammas 
-      else
-       {                                       
-         if (iPartMom>0)
-           {
-             if (pdgMom)
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   fHistElectronFirstMother->Fill(pdgMom->PdgCode());
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   fHistGammaFirstMother->Fill(pdgMom->PdgCode());
-                                       
-                 if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
-                   {
-                     if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                       fHistElectronFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
-                     else if (pdg->PdgCode() == fgGammaCode)
-                       fHistGammaFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
-                   }
-               }
-             else 
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   fHistElectronFirstMother->Fill(-598);
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   fHistGammaFirstMother->Fill(-598);
-               }
-           }
-         else
-           {
-             if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-               fHistElectronFirstMother->Fill(-599);
-             else if (pdg->PdgCode() == fgGammaCode)
-               fHistGammaFirstMother->Fill(-599);
-           }
-                       
-         if (iPartLastMom>0)
-           {
-             if (pdgMomLast)
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   fHistElectronLastMother->Fill(pdgMomLast->PdgCode());
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   fHistGammaLastMother->Fill(pdgMomLast->PdgCode());
-               }
-             else
-               {
-                 if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-                   fHistElectronLastMother->Fill(-598);
-                 else if (pdg->PdgCode() == fgGammaCode)
-                   fHistGammaLastMother->Fill(-598);
-               }                                               
-           }
-         else
-           {
-             if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
-               fHistElectronLastMother->Fill(-599);
-             else if (pdg->PdgCode() == fgGammaCode)
-               fHistGammaLastMother->Fill(-599);
-           }
-       }               
-
-      //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
-               
-      // Check for reasonable (for now neutral and singly charged) charge on the particle
-      //TODO:Maybe not only singly charged?
-      if (TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 && TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3) continue;
-               
-      fMultiplicity++;
-               
-      if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
+            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;
+
+        // it goes to the next particle in case it is not a physical primary or not a electron or gamma (want to check the contribution from secondary)
+        if (!stack->IsPhysicalPrimary(iPart))
         {
-         // 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
-             )
+            // this part is used only to check the origin of secondary electrons and gammas
+            if (iPartMom>0)
             {
-             if (pdg->PdgCode() > 0) { particleMassPart = - protonMass;}
-             if (pdg->PdgCode() < 0) { particleMassPart = protonMass;}
-           }
-         Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
-                       
-         // Fill up total E_T counters for each particle species
-         if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
-           {
-             fProtonEt += et;
-           }
-         if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
-           {
-             fPionEt += et;
-           }
-         if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
-           {
-             fChargedKaonEt += et;
-           }
-         if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
-           {
-             fMuonEt += et;
-           }
-         if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
-           {
-             fElectronEt += et;
-           }
-                       
-         // some neutrals also
-         if(pdg->PdgCode() == fgNeutronCode)
-           {
-             fNeutronEt += et;
-           }
-         if(pdg->PdgCode() == fgAntiNeutronCode)
-           {
-             fAntiNeutronEt += et;
-           }
-         if(pdg->PdgCode() == fgGammaCode)
-           {
-             fGammaEt += et;
-           }
-                       
-         // Neutral particles
-         if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+                //if (!stack->IsPhysicalPrimary(iPartMom)) continue;
+
+                if (pdgMom)
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                    {
+                        fHistElectronFirstMotherNPP->Fill(pdgMom->PdgCode());
+                        // inside EMCal acceptance
+                        if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
+                            fHistElectronFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
+                    }
+                    else if (pdg->PdgCode() == fgGammaCode)
+                    {
+                        fHistGammaFirstMotherNPP->Fill(pdgMom->PdgCode());
+                        if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
+                            fHistGammaFirstMotherNPPAcc->Fill(pdgMom->PdgCode());
+                    }
+                }
+                else
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                        fHistElectronFirstMotherNPP->Fill(-598);
+                    else if (pdg->PdgCode() == fgGammaCode)
+                        fHistGammaFirstMotherNPP->Fill(-598);
+
+                    continue;
+                }
+            }
+            else
             {
-             fNeutralMultiplicity++;
-             fTotNeutralEt += et;
-                               
-             // inside EMCal acceptance
-             if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
-                {
-                 fTotNeutralEtAcc += et;
-                 fTotEtAcc += et;                              
-
-                 if(pdg->PdgCode() == fgGammaCode)
-                   {
-                     if (!stack->IsPhysicalPrimary(iPart))
-                       {
-                         fHistSecGammaEtaE->Fill(part->Energy(),part->Eta());
-                         fHistSecGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                       }
-                     else 
-                       {
-                         if (pdgMom)
-                           {
-                             if (pdgMom->PdgCode() == fgPi0Code) 
-                               {
-                                 fHistPion0GammaEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistPion0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                             else if (pdgMom->PdgCode() == fgOmega0Code) 
-                               {
-                                 fHistOmega0GammaEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistOmega0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                             else if (pdgMom->PdgCode() == fgEtaCode)
-                               {
-                                 fHistEtaGammaEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistEtaGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                             else
-                               {
-                                 fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                           }
-                         else
-                           {
-                             fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
-                             fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
-                           }
-                       }
-                   }
-               }
+                if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                    fHistElectronFirstMotherNPP->Fill(-599);
+                else if (pdg->PdgCode() == fgGammaCode)
+                    fHistGammaFirstMotherNPP->Fill(-599);
+
+                continue;
             }
-         //Charged particles
-         else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
+
+            if (!((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode) || (pdg->PdgCode() == fgGammaCode)))
+                continue;
+        }
+        // this part is used only to check the origin of physical primary electrons and gammas
+        else
+        {
+            if (iPartMom>0)
+            {
+                if (pdgMom)
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                        fHistElectronFirstMother->Fill(pdgMom->PdgCode());
+                    else if (pdg->PdgCode() == fgGammaCode)
+                        fHistGammaFirstMother->Fill(pdgMom->PdgCode());
+
+                    if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
+                    {
+                        if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                            fHistElectronFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
+                        else if (pdg->PdgCode() == fgGammaCode)
+                            fHistGammaFirstMotherEtaAcc->Fill(pdgMom->PdgCode());
+                    }
+                }
+                else
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                        fHistElectronFirstMother->Fill(-598);
+                    else if (pdg->PdgCode() == fgGammaCode)
+                        fHistGammaFirstMother->Fill(-598);
+                }
+            }
+            else
+            {
+                if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                    fHistElectronFirstMother->Fill(-599);
+                else if (pdg->PdgCode() == fgGammaCode)
+                    fHistGammaFirstMother->Fill(-599);
+            }
+
+            if (iPartLastMom>0)
+            {
+                if (pdgMomLast)
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                        fHistElectronLastMother->Fill(pdgMomLast->PdgCode());
+                    else if (pdg->PdgCode() == fgGammaCode)
+                        fHistGammaLastMother->Fill(pdgMomLast->PdgCode());
+                }
+                else
+                {
+                    if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                        fHistElectronLastMother->Fill(-598);
+                    else if (pdg->PdgCode() == fgGammaCode)
+                        fHistGammaLastMother->Fill(-598);
+                }
+            }
+            else
+            {
+                if ((pdg->PdgCode() == fgEPlusCode) || (pdg->PdgCode() == fgEMinusCode))
+                    fHistElectronLastMother->Fill(-599);
+                else if (pdg->PdgCode() == fgGammaCode)
+                    fHistGammaLastMother->Fill(-599);
+            }
+        }
+
+        //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
+
+        // Check for reasonable (for now neutral and singly charged) charge on the particle
+        //TODO:Maybe not only singly charged?
+        if ((stack->IsPhysicalPrimary(iPart))&&( TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloSingleChargedParticle())<1e-3 || TMath::Abs(TMath::Abs(pdg->Charge()) - fCuts->GetMonteCarloNeutralParticle())<1e-3))
+        {
+
+            fMultiplicity++;
+
+            if (TMath::Abs(part->Eta()) < fCuts->GetCommonEtaCut())
             {
-             fChargedMultiplicity++;
-             fTotChargedEt += et;
-                               
-             // inside EMCal acceptance
-             if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
-                {
-                 fTotChargedEtAcc += et;
-                 fTotEtAcc += et;
-                                       
-                 if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
-                   {
-                     fProtonEtAcc += et;
-                   }
-                 if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
-                   {
-                     fPionEtAcc += et;
-                   }
-                 if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
-                   {
-                     fChargedKaonEtAcc += et;
-                   }
-                 if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
-                   {
-                     fMuonEtAcc += et;
-                   }
-                 /*
-                   if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
-                   {
-                   fElectronEtAcc += et;
-                   }
-                 */
-                 if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
-                   {
-                     if (!stack->IsPhysicalPrimary(iPart))
-                       {
-                         if (pdgMom)
-                           {
-                             if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
-                               {
-                                 fHistConvElectronEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistConvElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                             else
-                               {
-                                 fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
-                                 fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
-                               }
-                           }
-                         else
-                           {
-                             fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
-                             fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
-                           }                                                   
-                       }
-                     else 
-                       {
-                         fElectronEtAcc += et;
-                         fHistPrimElectronEtaE->Fill(part->Energy(),part->Eta());
-                         fHistPrimElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
-                       }
-                   } // electron
-               } // inside EMCal acceptance
-                               
-             //          if (TrackHitsCalorimeter(part, event->GetMagneticField()))
-             if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
-               {
-                 if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
-                 else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
-               }
-           }
-       }
+                // 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
+                )
+                {
+                    if (pdg->PdgCode() > 0) {
+                        particleMassPart = - protonMass;
+                    }
+                    if (pdg->PdgCode() < 0) {
+                        particleMassPart = protonMass;
+                    }
+                }
+                Double_t et = part->Energy() * TMath::Sin(part->Theta()) + particleMassPart;
+
+                fSparseTracks[0] = pdg->PdgCode();
+                fSparseTracks[1] = pdg->Charge()*3;
+                fSparseTracks[2] = pdg->Mass();
+                fSparseTracks[3] = et;
+                fSparseTracks[4] = part->Pt();
+                fSparseTracks[5] = part->Eta();
+                fSparseHistTracks->Fill(fSparseTracks);
+
+                // Fill up total E_T counters for each particle species
+                if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+                {
+                    fProtonEt += et;
+                }
+                if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+                {
+                    fPionEt += et;
+                    if (pdg->PdgCode() == fgPiPlusCode)
+                    {
+                        fPiPlusMult++;
+                    }
+                    else
+                    {
+                        fPiMinusMult++;
+                    }
+                }
+                if (pdg->PdgCode() == fgGammaCode)
+                {
+                    fPiZeroMult++;
+                }
+                if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+                {
+                    fChargedKaonEt += et;
+                }
+                if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+                {
+                    fMuonEt += et;
+                }
+                if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+                {
+                    fElectronEt += et;
+                }
+
+                // some neutrals also
+                if (pdg->PdgCode() == fgNeutronCode)
+                {
+                    fNeutronEt += et;
+                }
+                if (pdg->PdgCode() == fgAntiNeutronCode)
+                {
+                    fAntiNeutronEt += et;
+                }
+                if (pdg->PdgCode() == fgGammaCode)
+                {
+                    fGammaEt += et;
+                }
+
+                // Neutral particles
+                if (TMath::Abs(pdg->Charge() - fCuts->GetMonteCarloNeutralParticle()) <1e-3 )
+                {
+
+                    if (et > fCuts->GetCommonClusterEnergyCut()) fTotNeutralEt += et;
+
+                    // inside EMCal acceptance
+                    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;
+
+                        //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;
+                        if (part->Energy() > 0.05) partCount++;
+                        if (pdg->PdgCode() == fgGammaCode)
+                        {
+                            if (!stack->IsPhysicalPrimary(iPart))
+                            {
+                                fHistSecGammaEtaE->Fill(part->Energy(),part->Eta());
+                                fHistSecGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                            }
+                            else
+                            {
+                                if (pdgMom)
+                                {
+                                    if (pdgMom->PdgCode() == fgPi0Code)
+                                    {
+                                        fHistPion0GammaEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistPion0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                    else if (pdgMom->PdgCode() == fgOmega0Code)
+                                    {
+                                        fHistOmega0GammaEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistOmega0GammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                    else if (pdgMom->PdgCode() == fgEtaCode)
+                                    {
+                                        fHistEtaGammaEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistEtaGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                    else
+                                    {
+                                        fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                }
+                                else
+                                {
+                                    fHistPrimGammaEtaE->Fill(part->Energy(),part->Eta());
+                                    fHistPrimGammaEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                }
+                            }
+                        }
+                    }
+                }
+                //Charged particles
+                else if (TMath::Abs( pdg->Charge() - fCuts->GetMonteCarloNeutralParticle())>1e-3 )
+                {
+                    fChargedMultiplicity++;
+                    fTotChargedEt += et;
+
+                    // inside EMCal acceptance
+                    if (TMath::Abs(part->Eta()) < fEtaCutAcc && part->Phi() < fPhiCutAccMax && part->Phi() > fPhiCutAccMin)
+                    {
+                        fTotChargedEtAcc += et;
+                        fTotEtAcc += et;
+
+                        if (pdg->PdgCode() == fgProtonCode || pdg->PdgCode() == fgAntiProtonCode)
+                        {
+                            fProtonEtAcc += et;
+                        }
+                        if (pdg->PdgCode() == fgPiPlusCode || pdg->PdgCode() == fgPiMinusCode)
+                        {
+                            fPionEtAcc += et;
+                        }
+                        if (pdg->PdgCode() == fgKPlusCode || pdg->PdgCode() == fgKMinusCode)
+                        {
+                            fChargedKaonEtAcc += et;
+                        }
+                        if (pdg->PdgCode() == fgMuPlusCode || pdg->PdgCode() == fgMuMinusCode)
+                        {
+                            fMuonEtAcc += et;
+                        }
+                        /*
+                          if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+                          {
+                          fElectronEtAcc += et;
+                          }
+                        */
+                        if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+                        {
+                            if (!stack->IsPhysicalPrimary(iPart))
+                            {
+                                if (pdgMom)
+                                {
+                                    if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
+                                    {
+                                        fHistConvElectronEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistConvElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                    else
+                                    {
+                                        fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
+                                        fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                    }
+                                }
+                                else
+                                {
+                                    fHistSecElectronEtaE->Fill(part->Energy(),part->Eta());
+                                    fHistSecElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
+                                }
+                            }
+                            else
+                            {
+                                fElectronEtAcc += et;
+                                fHistPrimElectronEtaE->Fill(part->Energy(),part->Eta());
+                                fHistPrimElectronEtaEET->Fill(part->Energy(),part->Eta(),et);
+                            }
+                        } // electron
+                    } // inside EMCal acceptance
+
+                    //   if (TrackHitsCalorimeter(part, event->GetMagneticField()))
+                    if (TrackHitsCalorimeter(part)) // magnetic field info not filled?
+                    {
+                        if (pdg->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt());
+                        else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt());
+                    }
+                }
+            }
+
+        }
     }
-    
-  fTotEt = fTotChargedEt + fTotNeutralEt;
-  fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
-    
-  FillHistograms();
-       
-  return 0;    
+
+    fTotEt = fTotChargedEt + fTotNeutralEt;
+    fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
+
+
+    //std::cout << "Event done! # of particles: " << partCount << std::endl;
+    return 0;
 }
 //Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliMCEvent* mcEvent,AliESDEvent* realEvent)
 Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev,AliVEvent* ev2)
 { // analyse MC and real event info
-  //if(!mcEvent || !realEvent){
-  if(!ev || !ev2){
-    AliFatal("ERROR: Event does not exist");   
-    return 0;
-  }
-  //AnalyseEvent(mcEvent);
-  AnalyseEvent(ev);
-  AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
-  AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
-  if(!mcEvent || !realEvent){  
-    AliFatal("ERROR: mcEvent or realEvent does not exist");
-    return 0;
-  }
-
-  AliStack *stack = mcEvent->Stack();
-       
-  // get all emcal clusters
-  TRefArray* caloClusters = new TRefArray();
-  realEvent->GetEMCALClusters( caloClusters );
-       
-  Int_t nCluster = caloClusters->GetEntries();
-       
-  // loop the clusters
-  for (int iCluster = 0; iCluster < nCluster; iCluster++ ) 
+    //if(!mcEvent || !realEvent){
+    if (!ev || !ev2) {
+        AliFatal("ERROR: Event does not exist");
+        return 0;
+    }
+
+    fSparseClusters[0] = 0;
+    fSparseClusters[1] = 0;
+    fSparseClusters[2] = 0;
+    fSparseClusters[3] = 0;
+    fSparseClusters[4] = 0;
+    fSparseClusters[5] = 0;
+    fSparseClusters[6] = 0;
+    fSparseClusters[7] = 0;
+    fSparseClusters[8] = 0;
+    fSparseClusters[9] = fCentClass;
+    fSparseClusters[10] = 0;
+
+    //AnalyseEvent(mcEvent);
+    AnalyseEvent(ev);
+    AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(ev);
+    AliESDEvent *realEvent = dynamic_cast<AliESDEvent*>(ev2);
+    if (!mcEvent || !realEvent) {
+        AliFatal("ERROR: mcEvent or realEvent does not exist");
+        return 0;
+    }
+
+    AliStack *stack = mcEvent->Stack();
+
+    // get all detector clusters
+    TRefArray* caloClusters = new TRefArray();
+    if (fDetector == fCuts->GetDetectorEmcal()) realEvent->GetEMCALClusters( caloClusters );
+    else if (fDetector == fCuts->GetDetectorPhos()) realEvent->GetPHOSClusters( caloClusters );
+    else {
+        AliFatal("Detector ID has not been specified");
+        return -1;
+    }
+
+    // Track loop to fill a pT spectrum
+    for (Int_t iTracks = 0; iTracks < realEvent->GetNumberOfTracks(); iTracks++)
     {
-      AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
-      Float_t caloE = caloCluster->E();
-
-      UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
-      TParticle *part  = stack->Particle(iPart);
-      TParticle *partMom = 0;
-               
-      if (!part)
+        AliESDtrack* track = realEvent->GetTrack(iTracks);
+        if (!track)
         {
-         Printf("No MC particle %d", iCluster);
-         continue;
+            Printf("ERROR: Could not receive track %d", iTracks);
+            continue;
         }
-               
-      // compare MC and Rec energies for all particles
-      fHistAllERecEMC->Fill(part->Energy(),caloE);
-               
-      TParticlePDG *pdg = part->GetPDG(0);
-      TParticlePDG *pdgMom = 0;
-                               
-      if (!pdg)
+
+        if ( track->Phi() < fCuts->GetGeometryPhosPhiAccMaxCut() * TMath::Pi()/180. && track->Phi() > fCuts->GetGeometryPhosPhiAccMinCut()  * TMath::Pi()/180. && TMath::Abs(track->Eta()) < fCuts->GetGeometryPhosEtaAccCut())
+        {
+            fTrackMultInAcc++;
+        }
+    }
+
+    Int_t nCluster = caloClusters->GetEntries();
+
+    // loop the clusters
+    for (int iCluster = 0; iCluster < nCluster; iCluster++ )
+    {
+        AliESDCaloCluster* caloCluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
+        Float_t caloE = caloCluster->E();
+
+        UInt_t iPart = (UInt_t)TMath::Abs(caloCluster->GetLabel());
+        TParticle *part  = stack->Particle(iPart);
+        TParticle *partMom = 0;
+
+        if (!part)
         {
-         Printf("ERROR: Could not get particle PDG %d", iPart);
-         continue;
-        }              
-
-      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;
-      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;
-       } // 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;
-               
-      if(pdg->PdgCode() == fgGammaCode)
-       {
-         // compare MC and Rec energies for gammas
-         fHistGammaERecEMC->Fill(part->Energy(),caloE);
-                       
-         if (!stack->IsPhysicalPrimary(iPart))
-           {
-             fHistSecGammaEtaERec->Fill(part->Energy(),part->Eta());
-           }
-         else
-           {
-             if (pdgMom)
-               {
-                 if (pdgMom->PdgCode() == fgPi0Code) 
-                   {
-                     fHistPion0GammaEtaERec->Fill(part->Energy(),part->Eta());
-                   }
-                 else if (partMom->GetPDG(0)->PdgCode() == fgOmega0Code) 
-                   {
-                     fHistOmega0GammaEtaERec->Fill(part->Energy(),part->Eta());
-                   }
-                 else if (partMom->GetPDG(0)->PdgCode() == fgEtaCode)
-                   {
-                     fHistEtaGammaEtaERec->Fill(part->Energy(),part->Eta());
-                   }
-                 else
-                   {
-                     fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
-                   }
-               }
-             else
-               {
-                 fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
-               }
-           }
-       } // gamma
-
-      if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
-       {
-         // compare MC and Rec energies for electrons
-         fHistElectronERecEMC->Fill(part->Energy(),caloE);
-
-         if (!stack->IsPhysicalPrimary(iPart))
-           {
-             if (pdgMom)
-               {
-                 if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
-                   {
-                     fHistConvElectronEtaERec->Fill(part->Energy(),part->Eta());
-                   }
-                 else
-                   {
-                     fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());                                        
-                   }                                   
-               }
-             else
-               {
-                 fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());                                    
-               }
-           }
-         else 
-           {
-             fHistPrimElectronEtaERec->Fill(part->Energy(),part->Eta());
-           }
-       }
-    } // end of loop over clusters     
-  return 0;
-}      
+            Printf("No MC particle %d", iCluster);
+            continue;
+        }
+
+        // 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;
+        }
+
+        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;
+        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;
+        } // 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;
+
+        if (pdg->PdgCode() == fgGammaCode)
+        {
+            // compare MC and Rec energies for gammas
+            fHistGammaERecEMC->Fill(part->Energy(),caloE);
+
+            if (!stack->IsPhysicalPrimary(iPart))
+            {
+                fHistSecGammaEtaERec->Fill(part->Energy(),part->Eta());
+            }
+            else
+            {
+                if (pdgMom)
+                {
+                    if (pdgMom->PdgCode() == fgPi0Code)
+                    {
+                        fHistPion0GammaEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                    else if (partMom->GetPDG(0)->PdgCode() == fgOmega0Code)
+                    {
+                        fHistOmega0GammaEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                    else if (partMom->GetPDG(0)->PdgCode() == fgEtaCode)
+                    {
+                        fHistEtaGammaEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                    else
+                    {
+                        fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                }
+                else
+                {
+                    fHistPrimGammaEtaERec->Fill(part->Energy(),part->Eta());
+                }
+            }
+        } // gamma
+
+        if (pdg->PdgCode() == fgEPlusCode || pdg->PdgCode() == fgEMinusCode)
+        {
+            // compare MC and Rec energies for electrons
+            fHistElectronERecEMC->Fill(part->Energy(),caloE);
+
+            if (!stack->IsPhysicalPrimary(iPart))
+            {
+                if (pdgMom)
+                {
+                    if ((pdgMom->PdgCode() == fgGammaCode) && (stack->IsPhysicalPrimary(iPartMom)))
+                    {
+                        fHistConvElectronEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                    else
+                    {
+                        fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());
+                    }
+                }
+                else
+                {
+                    fHistSecElectronEtaERec->Fill(part->Energy(),part->Eta());
+                }
+            }
+            else
+            {
+                fHistPrimElectronEtaERec->Fill(part->Energy(),part->Eta());
+            }
+        }
+
+
+        Double_t clEt = CalculateTransverseEnergy(caloCluster);
+        if (caloCluster->E() < fCuts->GetCommonClusterEnergyCut()) continue;
+        Float_t pos[3];
+
+        caloCluster->GetPosition(pos);
+        TVector3 cp(pos);
+
+        fSparseClusters[0] = pdg->PdgCode();
+        fSparseClusters[1] = pdg->Charge()/3;
+        fSparseClusters[2] = pdg->Mass();
+        fSparseClusters[3] = clEt;
+        fSparseClusters[4] = caloCluster->E();
+        fSparseClusters[5] = cp.Eta();
+        fSparseClusters[6] = part->Energy() * TMath::Sin(part->Theta());
+        fSparseClusters[7] = part->Pt();
+        fSparseClusters[8] = part->Eta();
+        fSparseClusters[9] = fCentClass;
+        fSparseClusters[10] = 0;
+        fSparseHistClusters->Fill(fSparseClusters);
+
+        //if(caloCluster->GetEmcCpvDistance() < fTrackDistanceCut)
+
+        if (TMath::Abs(caloCluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(caloCluster->GetTrackDz()) < fTrackDzCut)
+        {
+
+            if (pdg->Charge() != 0)
+            {
+                //std::cout << "Removing charged: " << pdg->PdgCode() << ", with energy: " << caloCluster->E() << ", dist matched: " << caloCluster->GetTrackDx() << " " << caloCluster->GetTrackDz() << std::endl;
+                fChargedRemoved++;
+                fEnergyChargedRemoved += caloCluster->E();
+                fHistRemovedOrNot->Fill(0.0, fCentClass);
+                //std::cout << part->Vx() << " " << part->Vy() << " " << part->Vz() << " " << std::endl;
+                //fHistDecayVertexRemovedCharged->Fill(part->Vy(), part->Vx(), part->Vz());
+                fHistDxDzRemovedCharged->Fill(caloCluster->GetTrackDx(), caloCluster->GetTrackDz());
+
+
+            }
+            else
+            {
+                //std::cout << "Removing neutral: " << pdg->PdgCode() << ", 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();
+                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)
+                {
+                    fEtRemovedGammas += clEt;
+                    fMultRemovedGammas++;
+                }
+                else if (pdg->PdgCode() == fgNeutronCode)
+                {
+                    fEtRemovedNeutrons += clEt;
+                    fMultRemovedNeutrons++;
+                }
+                else if (pdg->PdgCode() == fgAntiNeutronCode)
+                {
+                    fEtRemovedAntiNeutrons += clEt;
+                    fMultRemovedAntiNeutrons++;
+                }
+                //else std::cout << "Hmm, what is this (neutral: " << pdg->PdgCode() << std::endl;
+            }
+        }
+        else
+        {
+
+            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 << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << std::endl;
+                fChargedNotRemoved++;
+                fEnergyChargedNotRemoved += caloCluster->E();
+                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)
+                {
+                    //std::cout << clEt << std::endl;
+                    fEtNonRemovedProtons += clEt;
+                    fMultNonRemovedProtons++;
+                }
+                else if (pdg->PdgCode() == fgAntiProtonCode)
+                {
+                    //std::cout << clEt << std::endl;
+                    fEtNonRemovedAntiProtons += clEt;
+                    fMultNonRemovedAntiProtons++;
+                }
+                else if (pdg->PdgCode() == fgPiPlusCode)
+                {
+                    //std::cout << "PI+" <<  clEt << std::endl;
+                    fEtNonRemovedPiPlus += clEt;
+                    fMultNonRemovedPiPlus++;
+                }
+                else if (pdg->PdgCode() == fgPiMinusCode)
+                {
+                    // std::cout << "PI-"  << clEt << std::endl;
+                    fEtNonRemovedPiMinus += clEt;
+                    fMultNonRemovedPiMinus++;
+                }
+                else if (pdg->PdgCode() == fgKPlusCode)
+                {
+                    //std::cout << clEt << std::endl;
+                    fEtNonRemovedKaonPlus += clEt;
+                    fMultNonRemovedKaonPlus++;
+                }
+                else if (pdg->PdgCode() == fgKMinusCode)
+                {
+                    //std::cout << clEt << std::endl;
+                    fEtNonRemovedKaonMinus += clEt;
+                    fMultNonRemovedKaonMinus++;
+                }
+                else if (pdg->PdgCode() == fgEPlusCode)
+                {
+                    //std::cout << clEt << std::endl;
+                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
+                    {
+                        fEtNonRemovedPositrons += clEt;
+                        fMultNonRemovedPositrons++;
+                    }
+                }
+                else if (pdg->PdgCode() == fgEMinusCode)
+                {
+                    //std::cout << clEt << std::endl;
+                    if (TMath::Sqrt(part->Vx()*part->Vx() + part->Vy()*part->Vy()) < 430)
+                    {
+                        fEtNonRemovedElectrons += clEt;
+                        fMultNonRemovedElectrons++;
+                    }
+                }
+                else if (pdg->PdgCode() == fgMuPlusCode)
+                {
+                    std::cout << clEt << std::endl;
+                    fEtNonRemovedMuPlus += clEt;
+                    fMultNonRemovedMuPlus++;
+                }
+                else if (pdg->PdgCode() == fgMuMinusCode)
+                {
+                    std::cout << clEt << std::endl;
+                    fEtNonRemovedMuMinus += clEt;
+                    fMultNonRemovedMuMinus++;
+                }
+
+            }
+            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 << "Mother is: " << stack->Particle(part->GetMother(0))->GetPdgCode() << stack->Particle(part->GetMother(0))->GetPDG()->GetName() << ", E: " << part->Energy() << std::endl;
+                fNeutralNotRemoved++;
+                fEnergyNeutralNotRemoved += caloCluster->E();
+                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)
+                {
+                    fEtNonRemovedGammas += clEt;
+                    fMultNonRemovedGammas++;
+                    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;
+                        }
+                    }
+                }
+                else if (pdg->PdgCode() == fgNeutronCode)
+                {
+                    fEtNonRemovedNeutrons += clEt;
+                    fMultNonRemovedNeutrons++;
+                }
+                else if (pdg->PdgCode() == fgAntiNeutronCode)
+                {
+                    fEtNonRemovedAntiNeutrons += clEt;
+                    fMultNonRemovedAntiNeutrons++;
+                }
+                else if (pdg->PdgCode() == fgK0LCode || pdg->PdgCode() == fgK0SCode)
+                {
+                    fEtNonRemovedK0s += clEt;
+                    fMultNonRemovedK0s++;
+
+                }
+                else if (TMath::Abs(pdg->PdgCode()) == fgLambdaCode)
+                {
+                    fEtNonRemovedLambdas += clEt;
+                    fMultNonRemovedLambdas++;
+                }
+                else std::cout << "Hmm, what is this (neutral not removed): " << pdg->PdgCode() << " " << pdg->GetName() << ", ET: " << clEt << std::endl;
+            }
+        }
+    } // end of loop over clusters
+
+    delete caloClusters;
+    //std::cout << "Distance cut: " << fTrackDistanceCut << std::endl;
+    //std::cout << "Number of removed neutrals: " << fNeutralRemoved << std::endl;
+    //std::cout << "Number of removed charged: " << fChargedRemoved << std::endl;
+    //std::cout << "Number of non-removed charged: " << fChargedNotRemoved << std::endl;
+    //std::cout << "Number of non-removed neutral: " << fNeutralNotRemoved << std::endl;
+
+//  std::cout << "Energy of removed neutrals: " << fEnergyNeutralRemoved << std::endl;
+//  std::cout << "Energy of removed charged: " << fEnergyChargedRemoved << std::endl;
+//  std::cout << "Energy of non-removed charged: " << fEnergyChargedNotRemoved << std::endl;
+//  std::cout << "Energy of non-removed neutral: " << fEnergyNeutralNotRemoved << std::endl;
+
+    FillHistograms();
+    return 0;
+}
 
 void AliAnalysisEtMonteCarlo::Init()
 { // init
-  AliAnalysisEt::Init();
+    AliAnalysisEt::Init();
 }
 
 void AliAnalysisEtMonteCarlo::ResetEventValues()
 { // reset event values
-  AliAnalysisEt::ResetEventValues();
-       
-  // collision geometry defaults for p+p:
-  fImpactParameter = 0;
-  fNcoll = 1;
-  fNpart = 2;  
+    AliAnalysisEt::ResetEventValues();
+
+    // collision geometry defaults for p+p:
+    fImpactParameter = 0;
+    fNcoll = 1;
+    fNpart = 2;
+
+    fEtNonRemovedProtons = 0;
+    fEtNonRemovedAntiProtons = 0;
+    fEtNonRemovedPiPlus = 0;
+    fEtNonRemovedPiMinus = 0;
+    fEtNonRemovedKaonPlus = 0;
+    fEtNonRemovedKaonMinus = 0;
+    fEtNonRemovedK0s = 0;
+    fEtNonRemovedLambdas = 0;
+    fEtNonRemovedElectrons = 0;
+    fEtNonRemovedPositrons = 0;
+    fEtNonRemovedMuPlus = 0;
+    fEtNonRemovedMuMinus = 0;
+    fEtNonRemovedNeutrons = 0;
+    fEtNonRemovedAntiNeutrons = 0;
+    fEtNonRemovedGammas = 0;
+    fEtNonRemovedGammasFromPi0 = 0;
+
+    fEtRemovedGammas = 0;
+    fEtRemovedNeutrons = 0;
+    fEtRemovedAntiNeutrons = 0;
+
+    fMultNonRemovedProtons = 0;
+    fMultNonRemovedAntiProtons = 0;
+    fMultNonRemovedPiPlus = 0;
+    fMultNonRemovedPiMinus = 0;
+    fMultNonRemovedKaonPlus = 0;
+    fMultNonRemovedKaonMinus = 0;
+    fMultNonRemovedK0s = 0;
+    fMultNonRemovedLambdas = 0;
+    fMultNonRemovedElectrons = 0;
+    fMultNonRemovedPositrons = 0;
+    fMultNonRemovedMuPlus = 0;
+    fMultNonRemovedMuMinus = 0;
+    fMultNonRemovedNeutrons = 0;
+    fMultNonRemovedAntiNeutrons = 0;
+    fMultNonRemovedGammas = 0;
+
+    fMultRemovedGammas = 0;
+    fMultRemovedNeutrons = 0;
+    fMultRemovedAntiNeutrons = 0;
+
+    fTrackMultInAcc = 0;
+
 }
 
 void AliAnalysisEtMonteCarlo::CreateHistograms()
 { // histogram related additions
-  AliAnalysisEt::CreateHistograms();
-  if (fTree) {
-    fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
-    fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
-    fTree->Branch("fNpart",&fNpart,"fNpart/I");
-  }
-  fHistPrimElectronEtaEET = CreateEtaEHisto2D("fHistPrimElectronEtaEET","MC E_{T}, primary electrons","E_{T}(GeV)");
-  fHistSecElectronEtaEET = CreateEtaEHisto2D("fHistSecElectronEtaEET","MC E_{T}, secondary (no conversion) electrons","E_{T}(GeV)");
-  fHistConvElectronEtaEET = CreateEtaEHisto2D("fHistConvElectronEtaEET","MC E_{T}, electrons from conversion","E_{T}(GeV)");
-  fHistPrimGammaEtaEET = CreateEtaEHisto2D("fHistPrimGammaEtaEET","MC E_{T}, primary gammas","E_{T}(GeV)"); 
-  fHistPion0GammaEtaEET = CreateEtaEHisto2D("fHistPion0GammaEtaEET","MC E_{T}, #pi^{0}","E_{T}(GeV)");
-  fHistOmega0GammaEtaEET = CreateEtaEHisto2D("fHistOmega0GammaEtaEET","MC E_{T}, #omega^{0}","E_{T}(GeV)");
-  fHistEtaGammaEtaEET = CreateEtaEHisto2D("fHistEtaGammaEtaEET","MC E_{T}, #eta","E_{T}(GeV)");
-  fHistSecGammaEtaEET = CreateEtaEHisto2D("fHistSecGammaEtaEET","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","E_{T}(GeV)"); 
-
-  fHistPrimElectronEtaE = CreateEtaEHisto2D("fHistPrimElectronEtaE","MC E_{T}, primary electrons","#");
-  fHistSecElectronEtaE = CreateEtaEHisto2D("fHistSecElectronEtaE","MC E_{T}, secondary (no conversion) electrons","#");
-  fHistConvElectronEtaE = CreateEtaEHisto2D("fHistConvElectronEtaE","MC E_{T}, electrons from conversion","#");
-  fHistPrimGammaEtaE = CreateEtaEHisto2D("fHistPrimGammaEtaE","MC E_{T}, primary gammas","#"); 
-  fHistPion0GammaEtaE = CreateEtaEHisto2D("fHistPion0GammaEtaE","MC E_{T}, #pi^{0}","#");
-  fHistOmega0GammaEtaE = CreateEtaEHisto2D("fHistOmega0GammaEtaE","MC E_{T}, #omega^{0}","#");
-  fHistEtaGammaEtaE = CreateEtaEHisto2D("fHistEtaGammaEtaE","MC E_{T}, #eta","#");
-  fHistSecGammaEtaE = CreateEtaEHisto2D("fHistSecGammaEtaE","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#"); 
-
-  fHistPrimElectronEtaERec = CreateEtaEHisto2D("fHistPrimElectronEtaERec","MC E_{T}, primary electrons","#");
-  fHistSecElectronEtaERec = CreateEtaEHisto2D("fHistSecElectronEtaERec","MC E_{T}, secondary (no conversion) electrons","#");
-  fHistConvElectronEtaERec = CreateEtaEHisto2D("fHistConvElectronEtaERec","MC E_{T}, electrons from conversion","#");
-  fHistPrimGammaEtaERec = CreateEtaEHisto2D("fHistPrimGammaEtaERec","MC E_{T}, primary gammas","#"); 
-  fHistPion0GammaEtaERec = CreateEtaEHisto2D("fHistPion0GammaEtaERec","MC E_{T}, #pi^{0}","#");
-  fHistOmega0GammaEtaERec = CreateEtaEHisto2D("fHistOmega0GammaEtaERec","MC E_{T}, #omega","#");
-  fHistEtaGammaEtaERec = CreateEtaEHisto2D("fHistEtaGammaEtaERec","MC E_{T}, #eta","#");
-  fHistSecGammaEtaERec = CreateEtaEHisto2D("fHistSecGammaEtaERec","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#"); 
-
-  fHistAllERecEMC = new TH2F("fHistAllERecEMC","E cluster Rec vs MC, all particles",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
-  fHistAllERecEMC->SetXTitle("E_{MC}(GeV)");
-  fHistAllERecEMC->SetYTitle("E_{Rec}(GeV)");
-
-  fHistElectronERecEMC = new TH2F("fHistElectronERecEMC","E cluster Rec vs MC, Electrons",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
-  fHistElectronERecEMC->SetXTitle("E_{MC}(GeV)");
-  fHistElectronERecEMC->SetYTitle("E_{Rec}(GeV)");
-       
-  fHistGammaERecEMC = new TH2F("fHistGammaERecEMC","E cluster Rec vs MC, Gammas",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
-  fHistGammaERecEMC->SetXTitle("E_{MC}(GeV)");
-  fHistGammaERecEMC->SetYTitle("E_{Rec}(GeV)");
-       
-  fHistElectronFirstMother = new TH1F("fHistElectronFirstMother","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistElectronLastMother = new TH1F("fHistElectronLastMother","Electron Last Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistElectronFirstMotherEtaAcc = new TH1F("fHistElectronFirstMotherEtaAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistElectronFirstMotherNPP = new TH1F("fHistElectronFirstMotherNPP","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistElectronFirstMotherNPPAcc = new TH1F("fHistElectronFirstMotherNPPAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
-
-  fHistGammaFirstMother = new TH1F("fHistGammaFirstMother","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistGammaLastMother = new TH1F("fHistGammaLastMother","Gamma Last Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistGammaFirstMotherEtaAcc = new TH1F("fHistGammaFirstMotherEtaAcc","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);  
-  fHistGammaFirstMotherNPP = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
-  fHistGammaFirstMotherNPPAcc = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);     
+    AliAnalysisEt::CreateHistograms();
+    if (fTree) {
+        fTree->Branch("fImpactParameter",&fImpactParameter,"fImpactParameter/D");
+        fTree->Branch("fNcoll",&fNcoll,"fNcoll/I");
+        fTree->Branch("fNpart",&fNpart,"fNpart/I");
+    }
+    fHistPrimElectronEtaEET = CreateEtaEHisto2D("fHistPrimElectronEtaEET","MC E_{T}, primary electrons","E_{T}(GeV)");
+    fHistSecElectronEtaEET = CreateEtaEHisto2D("fHistSecElectronEtaEET","MC E_{T}, secondary (no conversion) electrons","E_{T}(GeV)");
+    fHistConvElectronEtaEET = CreateEtaEHisto2D("fHistConvElectronEtaEET","MC E_{T}, electrons from conversion","E_{T}(GeV)");
+    fHistPrimGammaEtaEET = CreateEtaEHisto2D("fHistPrimGammaEtaEET","MC E_{T}, primary gammas","E_{T}(GeV)");
+    fHistPion0GammaEtaEET = CreateEtaEHisto2D("fHistPion0GammaEtaEET","MC E_{T}, #pi^{0}","E_{T}(GeV)");
+    fHistOmega0GammaEtaEET = CreateEtaEHisto2D("fHistOmega0GammaEtaEET","MC E_{T}, #omega^{0}","E_{T}(GeV)");
+    fHistEtaGammaEtaEET = CreateEtaEHisto2D("fHistEtaGammaEtaEET","MC E_{T}, #eta","E_{T}(GeV)");
+    fHistSecGammaEtaEET = CreateEtaEHisto2D("fHistSecGammaEtaEET","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","E_{T}(GeV)");
+
+    fHistPrimElectronEtaE = CreateEtaEHisto2D("fHistPrimElectronEtaE","MC E_{T}, primary electrons","#");
+    fHistSecElectronEtaE = CreateEtaEHisto2D("fHistSecElectronEtaE","MC E_{T}, secondary (no conversion) electrons","#");
+    fHistConvElectronEtaE = CreateEtaEHisto2D("fHistConvElectronEtaE","MC E_{T}, electrons from conversion","#");
+    fHistPrimGammaEtaE = CreateEtaEHisto2D("fHistPrimGammaEtaE","MC E_{T}, primary gammas","#");
+    fHistPion0GammaEtaE = CreateEtaEHisto2D("fHistPion0GammaEtaE","MC E_{T}, #pi^{0}","#");
+    fHistOmega0GammaEtaE = CreateEtaEHisto2D("fHistOmega0GammaEtaE","MC E_{T}, #omega^{0}","#");
+    fHistEtaGammaEtaE = CreateEtaEHisto2D("fHistEtaGammaEtaE","MC E_{T}, #eta","#");
+    fHistSecGammaEtaE = CreateEtaEHisto2D("fHistSecGammaEtaE","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#");
+
+    fHistPrimElectronEtaERec = CreateEtaEHisto2D("fHistPrimElectronEtaERec","MC E_{T}, primary electrons","#");
+    fHistSecElectronEtaERec = CreateEtaEHisto2D("fHistSecElectronEtaERec","MC E_{T}, secondary (no conversion) electrons","#");
+    fHistConvElectronEtaERec = CreateEtaEHisto2D("fHistConvElectronEtaERec","MC E_{T}, electrons from conversion","#");
+    fHistPrimGammaEtaERec = CreateEtaEHisto2D("fHistPrimGammaEtaERec","MC E_{T}, primary gammas","#");
+    fHistPion0GammaEtaERec = CreateEtaEHisto2D("fHistPion0GammaEtaERec","MC E_{T}, #pi^{0}","#");
+    fHistOmega0GammaEtaERec = CreateEtaEHisto2D("fHistOmega0GammaEtaERec","MC E_{T}, #omega","#");
+    fHistEtaGammaEtaERec = CreateEtaEHisto2D("fHistEtaGammaEtaERec","MC E_{T}, #eta","#");
+    fHistSecGammaEtaERec = CreateEtaEHisto2D("fHistSecGammaEtaERec","MC E_{T}, secondary (no #pi^{0}, #eta or #omega) gammas","#");
+
+    fHistAllERecEMC = new TH2F("fHistAllERecEMC","E cluster Rec vs MC, all particles",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
+    fHistAllERecEMC->SetXTitle("E_{MC}(GeV)");
+    fHistAllERecEMC->SetYTitle("E_{Rec}(GeV)");
+
+    fHistElectronERecEMC = new TH2F("fHistElectronERecEMC","E cluster Rec vs MC, Electrons",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
+    fHistElectronERecEMC->SetXTitle("E_{MC}(GeV)");
+    fHistElectronERecEMC->SetYTitle("E_{Rec}(GeV)");
+
+    fHistGammaERecEMC = new TH2F("fHistGammaERecEMC","E cluster Rec vs MC, Gammas",fgNumOfEBins, fgEAxis,fgNumOfEBins, fgEAxis);
+    fHistGammaERecEMC->SetXTitle("E_{MC}(GeV)");
+    fHistGammaERecEMC->SetYTitle("E_{Rec}(GeV)");
+
+    fHistElectronFirstMother = new TH1F("fHistElectronFirstMother","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistElectronLastMother = new TH1F("fHistElectronLastMother","Electron Last Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistElectronFirstMotherEtaAcc = new TH1F("fHistElectronFirstMotherEtaAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistElectronFirstMotherNPP = new TH1F("fHistElectronFirstMotherNPP","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistElectronFirstMotherNPPAcc = new TH1F("fHistElectronFirstMotherNPPAcc","Electron First Mother PDG Code Distribution",1201,-600.5,600.5);
+
+    fHistGammaFirstMother = new TH1F("fHistGammaFirstMother","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistGammaLastMother = new TH1F("fHistGammaLastMother","Gamma Last Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistGammaFirstMotherEtaAcc = new TH1F("fHistGammaFirstMotherEtaAcc","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistGammaFirstMotherNPP = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
+    fHistGammaFirstMotherNPPAcc = new TH1F("fHistGammaFirstMotherNPP","Gamma First Mother PDG Code Distribution",1201,-600.5,600.5);
+
+    //fHistDecayVertexNonRemovedCharged = new TH3F("fHistDecayVertexNonRemovedCharged","fHistDecayVertexNonRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
+    //fHistDecayVertexRemovedCharged = new TH3F("fHistDecayVertexRemovedCharged","fHistDecayVertexRemovedCharged", 500, -470, 30, 500, -300, 300, 40, -20, 20);
+    //fHistDecayVertexNonRemovedNeutral = new TH3F("fHistDecayVertexNonRemovedNeutral","fHistDecayVertexNonRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
+    //fHistDecayVertexRemovedNeutral = new TH3F("fHistDecayVertexRemovedNeutral","fHistDecayVertexRemovedNeutral", 500, -470, 30, 500, -300, 300, 40, -20, 20);
+
+    fHistRemovedOrNot = new TH2F("fHistRemovedOrNot", "fHistRemovedOrNot", 4, -0.5, 3.5, 10, -0.5, 9.5);
+
+    fHistEtNonRemovedProtons = new TH2F("fHistEtNonRemovedProtons", "fHistEtNonRemovedProtons", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedAntiProtons = new TH2F("fHistEtNonRemovedAntiProtons", "fHistEtNonRemovedAntiProtons", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedPiPlus = new TH2F("fHistEtNonRemovedPiPlus", "fHistEtNonRemovedPiPlus", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedPiMinus = new TH2F("fHistEtNonRemovedPiMinus", "fHistEtNonRemovedPiMinus", 1500, 0, 30, 10, -0.5, 9.5);
+    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);
+    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);
+    fHistEtNonRemovedMuPlus = new TH2F("fHistEtNonRemovedMuPlus", "fHistEtNonRemovedMuPlus", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedMuMinus = new TH2F("fHistEtNonRemovedMuMinus", "fHistEtNonRemovedMuMinus", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedNeutrons = new TH2F("fHistEtNonRemovedNeutrons", "fHistEtNonRemovedNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedAntiNeutrons = new TH2F("fHistEtNonRemovedAntiNeutrons", "fHistEtNonRemovedAntiNeutrons", 1500, 0, 30, 10, -0.5, 9.5);
+
+    fHistEtNonRemovedGammas = new  TH2F("fHistEtNonRemovedGammas", "fHistEtNonRemovedGammas", 1500, 0, 30, 10, -0.5, 9.5);
+    fHistEtNonRemovedGammasFromPi0 = new  TH2F("fHistEtNonRemovedGammasFromPi0", "fHistEtNonRemovedGammasFromPi0", 1500, 0, 30, 10, -0.5, 9.5);
+
+    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);
+
+    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);
+    fHistMultNonRemovedPiMinus = new TH2F("fHistMultNonRemovedPiMinus", "fHistMultNonRemovedPiMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    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);
+    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);
+    fHistMultNonRemovedMuPlus = new TH2F("fHistMultNonRemovedMuPlus", "fHistMultNonRemovedMuPlus", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    fHistMultNonRemovedMuMinus = new TH2F("fHistMultNonRemovedMuMinus", "fHistMultNonRemovedMuMinus", 100, -0.5, 99.5, 10, -0.5, 9.5);
+    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);
+
+    fHistMultRemovedGammas = new  TH2F("fHistMultRemovedGammas", "fHistMultRemovedGammas", 100, -0.5, 99.5, 10, -0.5, 9.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);
+
+    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);
+    fHistTrackMultvsRemovedGamma = new TH2F("fHistTrackMultvsRemovedGamma", "fHistTrackMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
+
+    fHistClusterMultvsNonRemovedCharged = new TH2F("fHistClusterMultvsNonRemovedCharged", "fHistClusterMultvsNonRemovedCharged", 1000, -0.5, 999.5, 100, -0.5, 99.5);
+    fHistClusterMultvsNonRemovedNeutral = new TH2F("fHistClusterMultvsNonRemovedNeutral", "fHistClusterMultvsNonRemovedNeutral", 1000, -0.5, 999.5, 100, -0.5, 99.5);
+    fHistClusterMultvsRemovedGamma = new TH2F("fHistClusterMultvsRemovedGamma", "fHistClusterMultvsRemovedGamma", 1000, -0.5, 999.5, 100, -0.5, 99.5);
+
+    fHistDxDzNonRemovedCharged = new TH2F("fHistDxDzNonRemovedCharged", "fHistDxDzNonRemovedCharged", 800, -200, 200, 800, -200, 200);
+    fHistDxDzRemovedCharged = new TH2F("fHistDxDzRemovedCharged", "fHistDxDzRemovedCharged", 800, -200, 200, 800, -200, 200);
+    fHistDxDzNonRemovedNeutral = new TH2F("fHistDxDzNonRemovedNeutral", "fHistDxDzNonRemovedNeutral", 800, -200, 200, 800, -200, 200);
+    fHistDxDzRemovedNeutral = new TH2F("fHistDxDzRemovedNeutral", "fHistDxDzRemovedNeutral", 800, -200, 200, 800, -200, 200);
+
+    fHistPiPlusMult = new TH1F("fHistPiPlusMult", "fHistPiPlusMult", 2000, -0.5, 1999.5);
+    fHistPiMinusMult = new TH1F("fHistPiMinusMult", "fHistPiMinusMult", 2000, -0.5, 1999.5);
+    fHistPiZeroMult = new TH1F("fHistPiZeroMult", "fHistPiZeroMult", 2000, -0.5, 1999.5);
+
+    fHistPiPlusMultAcc = new TH1F("fHistPiPlusMultAcc", "fHistPiPlusMultAcc", 2000, -0.5, 1999.5);
+    fHistPiMinusMultAcc = new TH1F("fHistPiMinusMultAcc", "fHistPiMinusMultAcc", 2000, -0.5, 1999.5);
+    fHistPiZeroMultAcc = new TH1F("fHistPiZeroMultAcc", "fHistPiZeroMultAcc", 2000, -0.5, 1999.5);
+
 }
 
 void AliAnalysisEtMonteCarlo::FillOutputList(TList *list)
 {//fill the output list
-  AliAnalysisEt::FillOutputList(list);
-
-  list->Add(fHistPrimElectronEtaEET);
-  list->Add(fHistSecElectronEtaEET);
-  list->Add(fHistConvElectronEtaEET);
-  list->Add(fHistPrimGammaEtaEET); 
-  list->Add(fHistPion0GammaEtaEET);
-  list->Add(fHistOmega0GammaEtaEET);
-  list->Add(fHistEtaGammaEtaEET);
-  list->Add(fHistSecGammaEtaEET); 
-       
-  list->Add(fHistPrimElectronEtaE);
-  list->Add(fHistSecElectronEtaE);
-  list->Add(fHistConvElectronEtaE);
-  list->Add(fHistPrimGammaEtaE); 
-  list->Add(fHistPion0GammaEtaE);
-  list->Add(fHistOmega0GammaEtaE);
-  list->Add(fHistEtaGammaEtaE);
-  list->Add(fHistSecGammaEtaE); 
-       
-  list->Add(fHistPrimElectronEtaERec);
-  list->Add(fHistSecElectronEtaERec);
-  list->Add(fHistConvElectronEtaERec);
-  list->Add(fHistPrimGammaEtaERec); 
-  list->Add(fHistPion0GammaEtaERec);
-  list->Add(fHistOmega0GammaEtaERec);
-  list->Add(fHistEtaGammaEtaERec);
-  list->Add(fHistSecGammaEtaERec); 
-       
-  list->Add(fHistAllERecEMC);
-  list->Add(fHistElectronERecEMC);
-  list->Add(fHistGammaERecEMC);
-       
-  list->Add(fHistElectronFirstMother);
-  list->Add(fHistElectronLastMother);
-  list->Add(fHistElectronFirstMotherEtaAcc);
-  list->Add(fHistElectronFirstMotherNPP);
-  list->Add(fHistElectronFirstMotherNPPAcc);
-
-  list->Add(fHistGammaFirstMother);
-  list->Add(fHistGammaLastMother);
-  list->Add(fHistGammaFirstMotherEtaAcc);
-  list->Add(fHistGammaFirstMotherNPP);
-  list->Add(fHistGammaFirstMotherNPPAcc);
+    AliAnalysisEt::FillOutputList(list);
+
+    list->Add(fHistPrimElectronEtaEET);
+    list->Add(fHistSecElectronEtaEET);
+    list->Add(fHistConvElectronEtaEET);
+    list->Add(fHistPrimGammaEtaEET);
+    list->Add(fHistPion0GammaEtaEET);
+    list->Add(fHistOmega0GammaEtaEET);
+    list->Add(fHistEtaGammaEtaEET);
+    list->Add(fHistSecGammaEtaEET);
+
+    list->Add(fHistPrimElectronEtaE);
+    list->Add(fHistSecElectronEtaE);
+    list->Add(fHistConvElectronEtaE);
+    list->Add(fHistPrimGammaEtaE);
+    list->Add(fHistPion0GammaEtaE);
+    list->Add(fHistOmega0GammaEtaE);
+    list->Add(fHistEtaGammaEtaE);
+    list->Add(fHistSecGammaEtaE);
+
+    list->Add(fHistPrimElectronEtaERec);
+    list->Add(fHistSecElectronEtaERec);
+    list->Add(fHistConvElectronEtaERec);
+    list->Add(fHistPrimGammaEtaERec);
+    list->Add(fHistPion0GammaEtaERec);
+    list->Add(fHistOmega0GammaEtaERec);
+    list->Add(fHistEtaGammaEtaERec);
+    list->Add(fHistSecGammaEtaERec);
+
+    list->Add(fHistAllERecEMC);
+    list->Add(fHistElectronERecEMC);
+    list->Add(fHistGammaERecEMC);
+
+    list->Add(fHistElectronFirstMother);
+    list->Add(fHistElectronLastMother);
+    list->Add(fHistElectronFirstMotherEtaAcc);
+    list->Add(fHistElectronFirstMotherNPP);
+    list->Add(fHistElectronFirstMotherNPPAcc);
+
+    list->Add(fHistGammaFirstMother);
+    list->Add(fHistGammaLastMother);
+    list->Add(fHistGammaFirstMotherEtaAcc);
+    list->Add(fHistGammaFirstMotherNPP);
+    list->Add(fHistGammaFirstMotherNPPAcc);
+
+    list->Add(fHistRemovedOrNot);
+
+    list->Add(fHistEtNonRemovedProtons);
+    list->Add(fHistEtNonRemovedAntiProtons);
+    list->Add(fHistEtNonRemovedPiPlus);
+    list->Add(fHistEtNonRemovedPiMinus);
+    list->Add(fHistEtNonRemovedKaonPlus);
+    list->Add(fHistEtNonRemovedKaonMinus);
+    list->Add(fHistEtNonRemovedK0s);
+    list->Add(fHistEtNonRemovedLambdas);
+    list->Add(fHistEtNonRemovedElectrons);
+    list->Add(fHistEtNonRemovedPositrons);
+    list->Add(fHistEtNonRemovedMuPlus);
+    list->Add(fHistEtNonRemovedMuMinus);
+    list->Add(fHistEtNonRemovedNeutrons);
+    list->Add(fHistEtNonRemovedAntiNeutrons);
+    list->Add(fHistEtNonRemovedGammas);
+    list->Add(fHistEtNonRemovedGammasFromPi0);
+
+    list->Add(fHistEtRemovedGammas);
+    list->Add(fHistEtRemovedNeutrons);
+    list->Add(fHistEtRemovedAntiNeutrons);
+
+
+    list->Add(fHistMultNonRemovedProtons);
+    list->Add(fHistMultNonRemovedAntiProtons);
+    list->Add(fHistMultNonRemovedPiPlus);
+    list->Add(fHistMultNonRemovedPiMinus);
+    list->Add(fHistMultNonRemovedKaonPlus);
+    list->Add(fHistMultNonRemovedKaonMinus);
+    list->Add(fHistMultNonRemovedK0s);
+    list->Add(fHistMultNonRemovedLambdas);
+    list->Add(fHistMultNonRemovedElectrons);
+    list->Add(fHistMultNonRemovedPositrons);
+    list->Add(fHistMultNonRemovedMuPlus);
+    list->Add(fHistMultNonRemovedMuMinus);
+    list->Add(fHistMultNonRemovedNeutrons);
+    list->Add(fHistMultNonRemovedAntiNeutrons);
+    list->Add(fHistMultNonRemovedGammas);
+
+    list->Add(fHistMultRemovedGammas);
+    list->Add(fHistMultRemovedNeutrons);
+    list->Add(fHistMultRemovedAntiNeutrons);
+
+    list->Add(fHistTrackMultvsNonRemovedCharged);
+    list->Add(fHistTrackMultvsNonRemovedNeutral);
+    list->Add(fHistTrackMultvsRemovedGamma);
+
+    list->Add(fHistClusterMultvsNonRemovedCharged);
+    list->Add(fHistClusterMultvsNonRemovedNeutral);
+    list->Add(fHistClusterMultvsRemovedGamma);
+
+    //list->Add(fHistDecayVertexNonRemovedCharged);
+    //list->Add(fHistDecayVertexNonRemovedNeutral);
+    //list->Add(fHistDecayVertexRemovedCharged);
+    //list->Add(fHistDecayVertexRemovedNeutral);
+
+    list->Add(fHistDxDzNonRemovedCharged);
+    list->Add(fHistDxDzRemovedCharged);
+    list->Add(fHistDxDzNonRemovedNeutral);
+    list->Add(fHistDxDzRemovedNeutral);
+
+    list->Add(fHistPiPlusMult);
+    list->Add(fHistPiMinusMult);
+    list->Add(fHistPiZeroMult);
+    list->Add(fHistPiPlusMultAcc);
+    list->Add(fHistPiMinusMultAcc);
+    list->Add(fHistPiZeroMultAcc);
+
 }
 
 
 bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField)
 {
-  //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
-  AliESDtrack *esdTrack = new AliESDtrack(part);
-  // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
-       
-  Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
-       
-  // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
-       
-  bool status = prop && 
-    TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && 
-    esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && 
-    esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
-  delete esdTrack;
-       
-  return status;
+    //  printf(" TrackHitsCalorimeter - magField %f\n", magField);
+    AliESDtrack *esdTrack = new AliESDtrack(part);
+    // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
+
+    Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
+
+    // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
+
+    bool status = prop &&
+                  TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
+                  esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
+                  esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
+    delete esdTrack;
+
+    return status;
+}
+
+void AliAnalysisEtMonteCarlo::FillHistograms()
+{
+    AliAnalysisEt::FillHistograms();
+    //std::cout << fEtNonRemovedPiPlus << " " << fCentClass << std::endl;
+
+    fHistEtNonRemovedProtons->Fill(fEtNonRemovedProtons, fCentClass);
+    fHistEtNonRemovedAntiProtons->Fill(fEtNonRemovedAntiProtons, fCentClass);
+    fHistEtNonRemovedKaonPlus->Fill(fEtNonRemovedKaonPlus, fCentClass);
+    fHistEtNonRemovedKaonMinus->Fill(fEtNonRemovedKaonMinus, fCentClass);
+    fHistEtNonRemovedK0s->Fill(fEtNonRemovedK0s, fCentClass);
+    fHistEtNonRemovedLambdas->Fill(fEtNonRemovedLambdas, fCentClass);
+    fHistEtNonRemovedPiPlus->Fill(fEtNonRemovedPiPlus, fCentClass);
+    fHistEtNonRemovedPiMinus->Fill(fEtNonRemovedPiMinus, fCentClass);
+    fHistEtNonRemovedElectrons->Fill(fEtNonRemovedElectrons, fCentClass);
+    fHistEtNonRemovedPositrons->Fill(fEtNonRemovedPositrons, fCentClass);
+    fHistEtNonRemovedMuPlus->Fill(fEtNonRemovedMuPlus, fCentClass);
+    fHistEtNonRemovedMuMinus->Fill(fEtNonRemovedMuMinus, fCentClass);
+    fHistEtNonRemovedNeutrons->Fill(fEtNonRemovedNeutrons, fCentClass);
+    fHistEtNonRemovedAntiNeutrons->Fill(fEtNonRemovedAntiNeutrons, fCentClass);
+    fHistEtNonRemovedGammas->Fill(fEtNonRemovedGammas, fCentClass);
+    fHistEtNonRemovedGammasFromPi0->Fill(fEtNonRemovedGammasFromPi0, fCentClass);
+
+    fHistEtRemovedGammas->Fill(fEtRemovedGammas, fCentClass);
+    fHistEtRemovedNeutrons->Fill(fEtRemovedNeutrons, fCentClass);
+    fHistEtRemovedAntiNeutrons->Fill(fEtRemovedAntiNeutrons, fCentClass);
+
+    fHistMultNonRemovedProtons->Fill(fMultNonRemovedProtons, fCentClass);
+    fHistMultNonRemovedAntiProtons->Fill(fMultNonRemovedAntiProtons, fCentClass);
+    fHistMultNonRemovedKaonPlus->Fill(fMultNonRemovedKaonPlus, fCentClass);
+    fHistMultNonRemovedKaonMinus->Fill(fMultNonRemovedKaonMinus, fCentClass);
+    fHistMultNonRemovedK0s->Fill(fMultNonRemovedK0s, fCentClass);
+    fHistMultNonRemovedLambdas->Fill(fMultNonRemovedLambdas, fCentClass);
+    fHistMultNonRemovedPiPlus->Fill(fMultNonRemovedPiPlus, fCentClass);
+    fHistMultNonRemovedPiMinus->Fill(fMultNonRemovedPiMinus, fCentClass);
+    fHistMultNonRemovedElectrons->Fill(fMultNonRemovedElectrons, fCentClass);
+    fHistMultNonRemovedPositrons->Fill(fMultNonRemovedPositrons, fCentClass);
+    fHistMultNonRemovedMuPlus->Fill(fMultNonRemovedMuPlus, fCentClass);
+    fHistMultNonRemovedMuMinus->Fill(fMultNonRemovedMuMinus, fCentClass);
+    fHistMultNonRemovedNeutrons->Fill(fMultNonRemovedNeutrons, fCentClass);
+    fHistMultNonRemovedAntiNeutrons->Fill(fMultNonRemovedAntiNeutrons, fCentClass);
+    fHistMultNonRemovedGammas->Fill(fMultNonRemovedGammas, fCentClass);
+
+    fHistMultRemovedGammas->Fill(fMultRemovedGammas, fCentClass);
+    fHistMultRemovedNeutrons->Fill(fMultRemovedNeutrons, fCentClass);
+    fHistMultRemovedAntiNeutrons->Fill(fMultRemovedAntiNeutrons, fCentClass);
+
+    fHistTrackMultvsNonRemovedCharged->Fill(fTrackMultInAcc,
+                                            fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus+fMultNonRemovedKaonPlus
+                                            +fMultNonRemovedMuMinus+fMultNonRemovedMuPlus+fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons
+                                            +fMultNonRemovedProtons);
+
+    fHistTrackMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
+                                            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+
+    fHistTrackMultvsRemovedGamma->Fill(fTrackMultInAcc,
+                                       fMultRemovedGammas);
+
+    fHistClusterMultvsNonRemovedCharged->Fill(fNeutralMultiplicity,
+            fMultNonRemovedAntiProtons+fMultNonRemovedElectrons+fMultNonRemovedKaonMinus
+            +fMultNonRemovedKaonPlus+fMultNonRemovedMuMinus+fMultNonRemovedMuPlus
+            +fMultNonRemovedPiMinus+fMultNonRemovedPiPlus+fMultNonRemovedPositrons+fMultNonRemovedProtons);
+
+    fHistClusterMultvsNonRemovedNeutral->Fill(fTrackMultInAcc,
+            fMultNonRemovedNeutrons+fMultNonRemovedAntiNeutrons+fMultNonRemovedK0s+fMultNonRemovedLambdas);
+
+    fHistClusterMultvsRemovedGamma->Fill(fTrackMultInAcc,
+                                         fMultRemovedGammas);
+
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
index bb67918..671ab57 100644 (file)
@@ -10,6 +10,7 @@
 
 #include "AliAnalysisEt.h"
 class TParticle;
+class TH3F;
 //class AliMCEvent;
 //class AliESDEvent;
 
@@ -30,6 +31,8 @@ public:
     virtual void CreateHistograms();
     virtual void FillOutputList(TList* list);
 
+    virtual void FillHistograms();
+
 protected:
 
     virtual bool TrackHitsCalorimeter(TParticle *part, Double_t magField=0.5);
@@ -83,6 +86,144 @@ protected:
        TH1F *fHistGammaFirstMotherNPP;//enter comment here
        TH1F *fHistGammaFirstMotherNPPAcc;//enter comment here
        
+       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; 
+       TH2F *fHistEtNonRemovedAntiProtons; 
+       TH2F *fHistEtNonRemovedPiPlus; 
+       TH2F *fHistEtNonRemovedPiMinus; 
+       TH2F *fHistEtNonRemovedKaonPlus; 
+       TH2F *fHistEtNonRemovedKaonMinus; 
+       TH2F *fHistEtNonRemovedK0s; 
+       TH2F *fHistEtNonRemovedLambdas; 
+       TH2F *fHistEtNonRemovedElectrons; 
+       TH2F *fHistEtNonRemovedPositrons; 
+       TH2F *fHistEtNonRemovedMuPlus;
+       TH2F *fHistEtNonRemovedMuMinus;
+       TH2F *fHistEtNonRemovedNeutrons;
+       TH2F *fHistEtNonRemovedAntiNeutrons;
+       TH2F *fHistEtNonRemovedGammas;
+       TH2F *fHistEtNonRemovedGammasFromPi0;
+       
+       TH2F *fHistEtRemovedGammas;
+       TH2F *fHistEtRemovedNeutrons;
+       TH2F *fHistEtRemovedAntiNeutrons;
+       
+       
+       TH2F *fHistMultNonRemovedProtons; 
+       TH2F *fHistMultNonRemovedAntiProtons; 
+       TH2F *fHistMultNonRemovedPiPlus; 
+       TH2F *fHistMultNonRemovedPiMinus; 
+       TH2F *fHistMultNonRemovedKaonPlus; 
+       TH2F *fHistMultNonRemovedKaonMinus; 
+       TH2F *fHistMultNonRemovedK0s; 
+       TH2F *fHistMultNonRemovedLambdas; 
+       TH2F *fHistMultNonRemovedElectrons; 
+       TH2F *fHistMultNonRemovedPositrons; 
+       TH2F *fHistMultNonRemovedMuPlus;
+       TH2F *fHistMultNonRemovedMuMinus;
+       TH2F *fHistMultNonRemovedNeutrons;
+       TH2F *fHistMultNonRemovedAntiNeutrons;
+       TH2F *fHistMultNonRemovedGammas;
+       
+       TH2F *fHistMultRemovedGammas;
+       TH2F *fHistMultRemovedNeutrons;
+       TH2F *fHistMultRemovedAntiNeutrons;
+       
+       TH2F *fHistTrackMultvsNonRemovedCharged;
+       TH2F *fHistTrackMultvsNonRemovedNeutral;
+       TH2F *fHistTrackMultvsRemovedGamma;
+       
+       TH2F *fHistClusterMultvsNonRemovedCharged;
+       TH2F *fHistClusterMultvsNonRemovedNeutral;
+       TH2F *fHistClusterMultvsRemovedGamma;
+       
+       TH2F *fHistMultvsNonRemovedChargedE;
+       TH2F *fHistMultvsNonRemovedNeutralE;
+       TH2F *fHistMultvsRemovedGammaE;
+       
+       Float_t fEtNonRemovedProtons; 
+       Float_t fEtNonRemovedAntiProtons; 
+       Float_t fEtNonRemovedPiPlus; 
+       Float_t fEtNonRemovedPiMinus; 
+       Float_t fEtNonRemovedKaonPlus; 
+       Float_t fEtNonRemovedKaonMinus; 
+       Float_t fEtNonRemovedK0s; 
+       Float_t fEtNonRemovedLambdas; 
+       Float_t fEtNonRemovedElectrons; 
+       Float_t fEtNonRemovedPositrons; 
+       Float_t fEtNonRemovedMuMinus;
+       Float_t fEtNonRemovedMuPlus;
+       Float_t fEtNonRemovedGammas;
+       Float_t fEtNonRemovedGammasFromPi0;
+       Float_t fEtNonRemovedNeutrons;
+       Float_t fEtNonRemovedAntiNeutrons;
+       
+       Float_t fEtRemovedGammas;
+       Float_t fEtRemovedNeutrons;
+       Float_t fEtRemovedAntiNeutrons;
+               
+       Int_t fMultNonRemovedProtons; 
+       Int_t fMultNonRemovedAntiProtons; 
+       Int_t fMultNonRemovedPiPlus; 
+       Int_t fMultNonRemovedPiMinus; 
+       Int_t fMultNonRemovedKaonPlus; 
+       Int_t fMultNonRemovedKaonMinus; 
+       Int_t fMultNonRemovedK0s; 
+       Int_t fMultNonRemovedLambdas; 
+       Int_t fMultNonRemovedElectrons; 
+       Int_t fMultNonRemovedPositrons; 
+       Int_t fMultNonRemovedMuMinus;
+       Int_t fMultNonRemovedMuPlus;
+       Int_t fMultNonRemovedGammas;
+       Int_t fMultNonRemovedNeutrons;
+       Int_t fMultNonRemovedAntiNeutrons;
+       
+       Int_t fMultRemovedGammas;
+       Int_t fMultRemovedNeutrons;
+       Int_t fMultRemovedAntiNeutrons;
+       
+       Int_t fTrackMultInAcc;
+       
+       
+       TH2F *fHistDxDzNonRemovedCharged;
+       TH2F *fHistDxDzRemovedCharged;
+       TH2F *fHistDxDzNonRemovedNeutral;
+       TH2F *fHistDxDzRemovedNeutral;
+       
+       TH1F *fHistPiPlusMult;
+       TH1F *fHistPiMinusMult;
+       TH1F *fHistPiZeroMult;
+
+       TH1F *fHistPiPlusMultAcc;
+       TH1F *fHistPiMinusMultAcc;
+       TH1F *fHistPiZeroMultAcc;
+       
+       Int_t fPiPlusMult;
+       Int_t fPiMinusMult;
+       Int_t fPiZeroMult;
+
+       Int_t fPiPlusMultAcc;
+       Int_t fPiMinusMultAcc;
+       Int_t fPiZeroMultAcc;
+       
+       
+       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:
 
     //Declare it private to avoid compilation warning