SPD multiplicity analysis (Domenico, Mariella)
authorenrico <enrico@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 29 Apr 2009 08:46:26 +0000 (08:46 +0000)
committerenrico <enrico@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 29 Apr 2009 08:46:26 +0000 (08:46 +0000)
PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx

index e4e1126..2d94023 100644 (file)
@@ -1,29 +1,52 @@
+/*************************************************************************
+* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+*                                                                        *
+* Author: The ALICE Off-line Project.                                    *
+* Contributors are mentioned in the code where appropriate.              *
+*                                                                        *
+* Permission to use, copy, modify and distribute this software and its   *
+* documentation strictly for non-commercial purposes is hereby granted   *
+* without fee, provided that the above copyright notice appears in all   *
+* copies and that both the copyright notice and this permission notice   *
+* appear in the supporting documentation. The authors make no claims     *
+* about the suitability of this software for any purpose. It is          *
+* provided "as is" without express or implied warranty.                  * 
+**************************************************************************/
+
+/*************************************************************************
+* Class AliAnalysisTaskSPDdNdEta                                         *
+* Analysis task for dN/dEta reconstruction with the SPD                  *
+*                                                                        *
+* Author:  M. Nicassio (INFN Bari)                                       *
+* Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it           *
+**************************************************************************/
+
 #include "TChain.h"
 #include "TTree.h"
 #include "TBranch.h"
-#include <TClonesArray.h>
+#include "TClonesArray.h"
 
 #include "TH1F.h"
 #include "TH2F.h" 
 #include "TH3F.h" 
 #include "TCanvas.h"
 
-#include <AliAnalysisTask.h>
-#include <AliAnalysisManager.h>
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
 
-#include <AliMultiplicity.h>
-#include <AliESDVertex.h>
-#include <AliESDEvent.h>
-#include <AliESDInputHandler.h>
+#include "AliMultiplicity.h"
+#include "AliESDVertex.h"
+#include "AliESDEvent.h"
+#include "AliESDInputHandler.h"
 
 #include "AliAnalysisTaskRL.h"
 #include "AliMCEventHandler.h"
 #include "AliMCEvent.h"
-#include <AliStack.h>
-#include <AliLog.h>
+#include "AliStack.h"
+#include "AliLog.h"
 #include "AliTrackReference.h"
 
-#include <AliGenEventHeader.h> 
+#include "AliGenEventHeader.h" 
 #include "AliAnalysisTaskSPDdNdEta.h"
 
 
@@ -38,13 +61,14 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
 
   fCorr(kFALSE), 
   fTrigger(0),
+  fppAna(0),
 
-  //Data to be corrected
+  // Data to be corrected
   fHistSPDRAWMultvsZ(0),
   fHistSPDRAWMultvsZTriggEvts(0),
   fHistSPDRAWEtavsZ(0),
 
-  //Clusters inner layer and tracklets
+  // Clusters inner layer and tracklets
   fHistSPDmultEtacut(0),
   fHistSPDmult(0),
   fHistSPDeta(0),               
@@ -61,7 +85,7 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
   fHistSPDphivsSPDeta(0),
   fHistSPDcl1phivsSPDcl1eta(0),
 
-  //SPD vertex 
+  // SPD vertex 
   fHistSPDvtx(0),                 
   fHistSPDvtx3D(0),
   fHistSPDvtxZ(0),
@@ -71,14 +95,14 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
   fHistNcontribSPDvtxall(0),
   fHistSPDmultvsSPDvtx(0),
 
-  //SPD fired chips  
+  // SPD fired chips  
   fHistSPDcl1multvsnFiredChipsLay1(0),
   fHistSPDmultvsnFiredChipsLay1(0),
   fHistSPDmultvsnFiredChipsLay2(0),
   fHistnFiredChipsLay2vsnFiredChipsLay1(0),
   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec(0),
 
-  //Track level correction histograms
+  // Track level correction histograms
   fHistBkgCorrNum(0),   
   fHistBkgCorrDen(0),
 
@@ -90,18 +114,18 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
   fHistTrackTrigVtxCorrNum(0),   
   fHistTrackTrigCorrDen(0),      
 
-  //Event level correction histograms
+  // Event level correction histograms
   fHistTrigVtxCorrNum(0),           
   fHistTrigVtxCorrDen(0),           
   
   fHistTrigCorrDen(0),              
 
-  // MC distribustions
+  // MC distributions
   fHistMCEtavsZTriggMCvtxEvts(0),
   fHistMCEtavsZTriggESDvtxEvts(0),
   fHistMCEtavsZ(0),                
 
-  // Check histos
+  // Additional check histos
   fHistTRradius(0),
   fHistContributorsvsDeVtx(0),
   fHistoDetectableNotr(0),
@@ -111,13 +135,17 @@ AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
   fHistoDetectedLay2(0),
   fHistoPt(0),
   fHistoDetectableTRm1(0),
-  fHistoDetectableTR0(0),
-  fHistoDetectableTR1(0),
-  fHistoDetectableTR2(0),
-  fHistoDetectableTR3(0),
-  fHistoDetectableTR4(0),
-  fHistoDetectableTR5(0),
-  fHistoDetectableTR6(0),   
+  fHistoDetectableTrITS(0),
+  fHistoDetectableTrTPC(0),
+  fHistoDetectableTrFRAME(0),
+  fHistoDetectableTrTRD(0),
+  fHistoDetectableTrTOF(0),
+  fHistoDetectableTrMUON(0),
+  fHistoDetectableTrHMPID(0),
+  fHistoDetectableTrT0(0),
+  fHistoDetectableTrEMCAL(0),  
+  fHistoDetectableTrFMD(0), 
+  fHistoDetectableTrVZERO(0), 
   fHistoRTRm1(0)
   
 {
@@ -162,7 +190,7 @@ void AliAnalysisTaskSPDdNdEta::ConnectInputData(Option_t *)
       Printf("ERROR: Could not get ESDInputHandler");
     } else fESD = esdH->GetEvent();
   }
-  // disable info messages of AliMCEvent (per event)
+  // Disable info messages of AliMCEvent (per event)
   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);  
 
 }
@@ -175,11 +203,22 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
   fOutput = new TList;
   fOutput->SetOwner();
 
-  fHistSPDRAWMultvsZ= new TH2F("fHistSPDRAWMultvsZ", "",200,0.,200.,20,-20.,20.);
+  Int_t nBinMult = 200;
+  Float_t lowBinLim = 0.;
+  Float_t highBinLim = 200.;
+
+  if (!fppAna) {
+    nBinMult = 2000;
+    lowBinLim = 0.;
+    highBinLim = 80000.;    
+  } 
+
+  // Event level correction   
+  fHistSPDRAWMultvsZ= new TH2F("fHistSPDRAWMultvsZ", "",nBinMult,lowBinLim,highBinLim,20,-20.,20.);
   fHistSPDRAWMultvsZ->Sumw2();
   fOutput->Add(fHistSPDRAWMultvsZ);
 
-  fHistSPDRAWMultvsZTriggEvts = new TH2F("fHistSPDRAWMultvsZTriggEvts", "",200,0.,200.,20,-20.,20.);
+  fHistSPDRAWMultvsZTriggEvts = new TH2F("fHistSPDRAWMultvsZTriggEvts", "",nBinMult,lowBinLim,highBinLim,20,-20.,20.);
   fHistSPDRAWMultvsZTriggEvts->Sumw2();
   fOutput->Add(fHistSPDRAWMultvsZTriggEvts);
 
@@ -188,14 +227,13 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
   fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
   fHistSPDRAWEtavsZ->Sumw2();
   fOutput->Add(fHistSPDRAWEtavsZ);
-
   
-  fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution", 200, 0., 200);
+  fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution", nBinMult,lowBinLim,highBinLim);
   fHistSPDmultEtacut->GetXaxis()->SetTitle("Reconstructed multiplicity (|#eta|<1.5)");
   fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistSPDmultEtacut);
 
-  fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", 200, 0., 200);
+  fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", nBinMult,lowBinLim,highBinLim);
   fHistSPDmult->GetXaxis()->SetTitle("Reconstructed tracklet multiplicity");
   fHistSPDmult->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistSPDmult);
@@ -208,12 +246,12 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
   fHistSPDeta->Sumw2();
   fOutput->Add(fHistSPDeta); 
 
-  fHistSPDcl1multEtacutLay1 = new TH1F("fHistSPDcl1multEtacutLay1", "Cluster multiplicity (inner layer)", 200, 0., 200);
+  fHistSPDcl1multEtacutLay1 = new TH1F("fHistSPDcl1multEtacutLay1", "Cluster multiplicity (inner layer)", nBinMult,lowBinLim,highBinLim);
   fHistSPDcl1multEtacutLay1->GetXaxis()->SetTitle("Cluster multiplicity lay1 (|#eta|<2.)");
   fHistSPDcl1multEtacutLay1->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistSPDcl1multEtacutLay1);
 
-  fHistSPDcl1mult = new TH1F("fHistSPDcl1mult", "Cluster multiplicity (inner layer)", 200, 0., 200);
+  fHistSPDcl1mult = new TH1F("fHistSPDcl1mult", "Cluster multiplicity (inner layer)", nBinMult,lowBinLim,highBinLim);
   fHistSPDcl1mult->GetXaxis()->SetTitle("Cluster multiplicity lay1");
   fHistSPDcl1mult->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistSPDcl1mult);
@@ -265,7 +303,6 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
   fHistSPDcl1phivsSPDcl1eta->GetYaxis()->SetTitle("#varphi [rad]");
   fOutput->Add(fHistSPDcl1phivsSPDcl1eta);
 
-
   fHistSPDvtx = new TH1F("fHistSPDvtx", "SPD vertex  distribution - all events",20,-20.,20.);
   fHistSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
   fHistSPDvtx->GetYaxis()->SetTitle("Entries");
@@ -277,52 +314,52 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
   fHistSPDvtxZ = new TH3F("fHistSPDvtxZ", "SPD vertex distribution",100,-5.,5.,100,-5.,5.,400,-20.,20.);
   fOutput->Add(fHistSPDvtxZ);
 
-  fHistNcontribSPDvtxvsSPDvtx= new TH2F("fHistNcontribSPDvtxvsSPDvtx", " ",100,-50.,50.,202,-2.,200.);
+  fHistNcontribSPDvtxvsSPDvtx= new TH2F("fHistNcontribSPDvtxvsSPDvtx", " ",100,-50.,50.,10002,-2.,10000.);
   fHistNcontribSPDvtxvsSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
   fHistNcontribSPDvtxvsSPDvtx->GetYaxis()->SetTitle("# contributors");
   fOutput->Add(fHistNcontribSPDvtxvsSPDvtx);
 
-  fHistNcontribSPDvtx3D= new TH1F("fHistNcontribSPDvtx_3D", "SPD vtx 3D",202,-2.,200.);
+  fHistNcontribSPDvtx3D= new TH1F("fHistNcontribSPDvtx_3D", "SPD vtx 3D",10002,-2.,10000.);
   fHistNcontribSPDvtx3D->GetXaxis()->SetTitle("# contributors");
   fHistNcontribSPDvtx3D->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistNcontribSPDvtx3D);
 
-  fHistNcontribSPDvtxZ= new TH1F("fHistNcontribSPDvtxZ", "SPD vtx Z",202,-2.,200.);
+  fHistNcontribSPDvtxZ= new TH1F("fHistNcontribSPDvtxZ", "SPD vtx Z",10002,-2.,10000.);
   fHistNcontribSPDvtxZ->GetXaxis()->SetTitle("# contributors");
   fHistNcontribSPDvtxZ->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistNcontribSPDvtxZ);
 
-  fHistNcontribSPDvtxall= new TH1F("fHistNcontribSPDvtxall", "SPD vtx - all events",202,-2.,200.);
+  fHistNcontribSPDvtxall= new TH1F("fHistNcontribSPDvtxall", "SPD vtx - all events",10002,-2.,10000.);
   fHistNcontribSPDvtxall->GetXaxis()->SetTitle("# contributors");
   fHistNcontribSPDvtxall->GetYaxis()->SetTitle("Entries");
   fOutput->Add(fHistNcontribSPDvtxall);
 
-  fHistSPDmultvsSPDvtx= new TH2F("fHistSPDmultvsSPDvtx", "",20,-20.,20.,200,0.,200.);
+  fHistSPDmultvsSPDvtx= new TH2F("fHistSPDmultvsSPDvtx", "",20,-20.,20.,nBinMult,lowBinLim,highBinLim);
   fHistSPDmultvsSPDvtx->GetXaxis()->SetTitle("z_{recvtx} [cm]");
   fHistSPDmultvsSPDvtx->GetYaxis()->SetTitle("Reconstructed multiplicity");
   fOutput->Add(fHistSPDmultvsSPDvtx);
 
-  fHistSPDcl1multvsnFiredChipsLay1 = new TH2F("fHistSPDcl1multvsnFiredChipsLay1", "",200,0.,200.,200,0.,200.);
+  fHistSPDcl1multvsnFiredChipsLay1 = new TH2F("fHistSPDcl1multvsnFiredChipsLay1", "",401,0.,401.,nBinMult,lowBinLim,highBinLim);
   fHistSPDcl1multvsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
   fHistSPDcl1multvsnFiredChipsLay1->GetYaxis()->SetTitle("Cluster lay1 multiplicity");
   fOutput->Add(fHistSPDcl1multvsnFiredChipsLay1);
 
-  fHistSPDmultvsnFiredChipsLay1 = new TH2F("fHistSPDmultvsnFiredChipsLay1","",200,0.,200.,200,0.,200.);
+  fHistSPDmultvsnFiredChipsLay1 = new TH2F("fHistSPDmultvsnFiredChipsLay1","",401,0.,401.,nBinMult,lowBinLim,highBinLim);
   fHistSPDmultvsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
   fHistSPDmultvsnFiredChipsLay1->GetYaxis()->SetTitle("Tracklet multiplicity");
   fOutput->Add(fHistSPDmultvsnFiredChipsLay1);
 
-  fHistSPDmultvsnFiredChipsLay2 = new TH2F("fHistSPDmultvsnFiredChipsLay2","",200,0.,200.,200,0.,200.);
+  fHistSPDmultvsnFiredChipsLay2 = new TH2F("fHistSPDmultvsnFiredChipsLay2","",801,0.,801.,nBinMult,lowBinLim,highBinLim);
   fHistSPDmultvsnFiredChipsLay2->GetXaxis()->SetTitle("# fired chips lay2");
   fHistSPDmultvsnFiredChipsLay2->GetYaxis()->SetTitle("Tracklet multiplicity");
   fOutput->Add(fHistSPDmultvsnFiredChipsLay2);
 
-  fHistnFiredChipsLay2vsnFiredChipsLay1 = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1","",200,0.,200.,200,0.,200.);
+  fHistnFiredChipsLay2vsnFiredChipsLay1 = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1","",401,0.,401.,801,0.,801.);
   fHistnFiredChipsLay2vsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
   fHistnFiredChipsLay2vsnFiredChipsLay1->GetYaxis()->SetTitle("# fired chip lay2");
   fOutput->Add(fHistnFiredChipsLay2vsnFiredChipsLay1);
 
-  fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec","",200,0.,200.,200,0.,200.);
+  fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec","",401,0.,401.,801,0.,801.);
   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec->GetXaxis()->SetTitle("# fired chips lay1");
   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec->GetYaxis()->SetTitle("# fired chip lay2");
   fOutput->Add(fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec);
@@ -357,16 +394,16 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
     fHistTrackTrigCorrDen->Sumw2();
     fOutput->Add(fHistTrackTrigCorrDen);
 
-    //Event level correction histograms
-    fHistTrigVtxCorrNum = new TH2F("fHistTrigVtxCorrNum","",200,0.,200.,20,-20.,20.);
+    // Event level correction histograms  
+    fHistTrigVtxCorrNum = new TH2F("fHistTrigVtxCorrNum","",nBinMult,lowBinLim,highBinLim,20,-20.,20.);
     fHistTrigVtxCorrNum->Sumw2();
     fOutput->Add(fHistTrigVtxCorrNum);
 
-    fHistTrigVtxCorrDen = new TH2F("fHistTrigVtxCorrDen","",200,0.,200.,20,-20.,20.);
+    fHistTrigVtxCorrDen = new TH2F("fHistTrigVtxCorrDen","",nBinMult,lowBinLim,highBinLim,20,-20.,20.);
     fHistTrigVtxCorrDen->Sumw2();
     fOutput->Add(fHistTrigVtxCorrDen);
 
-    fHistTrigCorrDen = new TH2F("fHistTrigCorrDen","",200,0.,200.,20,-20.,20.);
+    fHistTrigCorrDen = new TH2F("fHistTrigCorrDen","",nBinMult,lowBinLim,highBinLim,20,-20.,20.);
     fHistTrigCorrDen->Sumw2();
     fOutput->Add(fHistTrigCorrDen);
 
@@ -392,9 +429,9 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
     fHistTRradius = new TH1F("fHistTRradius","ITS track reference rad",200,0.,10.);
     fOutput->Add(fHistTRradius); 
  
-    fHistContributorsvsDeVtx = new TH2F("fHistContributorsvsDeVtx","",200,-20.,20.,22,-2.,20.);
+    fHistContributorsvsDeVtx = new TH2F("fHistContributorsvsDeVtx","",200,-20.,20.,202,-2.,200.);
     fOutput->Add(fHistContributorsvsDeVtx);
+
     fHistoDetectableNotr = new TH3F("fHistoDetectableNotr","",60,-3.00,3.00,20,-20.,20.,100,0.,10.); 
     fOutput->Add(fHistoDetectableNotr);
  
@@ -416,27 +453,39 @@ void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
     fHistoDetectableTRm1 = new TH2F("fHistoDetectableTRm1","",60,-3.00,3.00,20,-20.,20.);
     fOutput->Add(fHistoDetectableTRm1);
 
-    fHistoDetectableTR0 = new TH2F("fHistoDetectableTR0","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR0);
+    fHistoDetectableTrITS = new TH2F("fHistoDetectableTrITS","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrITS);
  
-    fHistoDetectableTR1 = new TH2F("fHistoDetectableTR1","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR1);
+    fHistoDetectableTrTPC = new TH2F("fHistoDetectableTrTPC","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrTPC);
 
-    fHistoDetectableTR2 = new TH2F("fHistoDetectableTR2","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR2);
+    fHistoDetectableTrFRAME = new TH2F("fHistoDetectableTrFRAME","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrFRAME);
  
-    fHistoDetectableTR3 = new TH2F("fHistoDetectableTR3","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR3);
+    fHistoDetectableTrTRD = new TH2F("fHistoDetectableTrTRD","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrTRD);
+
+    fHistoDetectableTrTOF = new TH2F("fHistoDetectableTrTOF","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrTOF);
+
+    fHistoDetectableTrMUON = new TH2F("fHistoDetectableTrMUON","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrMUON);
+
+    fHistoDetectableTrHMPID = new TH2F("fHistoDetectableTrHMPID","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrHMPID);
 
-    fHistoDetectableTR4 = new TH2F("fHistoDetectableTR4","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR4);
+    fHistoDetectableTrT0 = new TH2F("fHistoDetectableTrT0","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrT0);
 
-    fHistoDetectableTR5 = new TH2F("fHistoDetectableTR5","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR5);
+    fHistoDetectableTrEMCAL = new TH2F("fHistoDetectableTrEMCAL","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrEMCAL);
 
-    fHistoDetectableTR6 = new TH2F("fHistoDetectableTR6","",60,-3.00,3.00,20,-20.,20.);
-    fOutput->Add(fHistoDetectableTR6);
+    fHistoDetectableTrFMD = new TH2F("fHistoDetectableTrFMD","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrFMD);
  
+    fHistoDetectableTrVZERO = new TH2F("fHistoDetectableTrVZERO","",60,-3.00,3.00,20,-20.,20.);
+    fOutput->Add(fHistoDetectableTrVZERO);
+
     fHistoRTRm1 = new TH1F("fHistoRTRm1","",10000,0.,5000);
     fOutput->Add(fHistoRTRm1);
 
@@ -497,7 +546,8 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
   Float_t* recEtaSPDcl1 = new Float_t[multSPD+nSingleCl1];
 
 //  Printf("There are %d tracklets in this event", multSPD);
-  if (esdvtx[2]!=0.&&eventTriggered&&multSPD!=0) {
+  // Selected events: triggered with vertex
+  if (esdvtx[2]!=0.&&eventTriggered &&multSPD!=0) {
 
     fHistSPDvtx->Fill(esdvtx[2]);
     if (strcmp(vtxESD->GetTitle(),"vertexer: 3D") == 0) fHistSPDvtx3D->Fill(esdvtx[0],esdvtx[1],esdvtx[2]);
@@ -523,14 +573,14 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
       fHistSPDphivsSPDeta->Fill(recEtaSPD,phiTr);
       fHistSPDcl1phivsSPDcl1eta->Fill(recEtaSPD,phiTr);
 
-      // calculate multiplicity in etacut
+      // Calculate multiplicity in etacut
       if (TMath::Abs(recEtaSPD)<1.5) multSPDEtacut++;
       if (TMath::Abs(recEtaSPD)<2.)  multSPDcl1EtacutLay1++;
     }
 
     for (Int_t iCl1=0; iCl1<nSingleCl1; ++iCl1) {
       Float_t thetaSingleCl1 = multESD->GetThetaSingle(iCl1);
-      //calculate eta
+      // Calculate eta
       Float_t etaSingleCl1 = -TMath::Log(TMath::Tan(thetaSingleCl1/2.));
       Float_t phiSingleCl1 = multESD->GetPhiSingle(iCl1);
       recEtaSPDcl1[iCl1+multSPD] = etaSingleCl1;
@@ -557,7 +607,7 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
     fHistSPDmultvsnFiredChipsLay1->Fill(nFiredChipsLay1,multSPD);
     fHistSPDmultvsnFiredChipsLay2->Fill(nFiredChipsLay2,multSPD);
 
-  } //ev selection
+  } // End selected events
 
   if (eventTriggered) fHistSPDRAWMultvsZTriggEvts->Fill(multSPD,esdvtx[2]);
 
@@ -599,12 +649,11 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
       return;
     }
     AliGenEventHeader* genHeader = header->GenEventHeader();
-    //MC vertex
+    // MC vertex
     TArrayF vtxMC(3);
     genHeader->PrimaryVertex(vtxMC);
 
-
-    //Tracks from MC
+    // Tracks from MC
     Int_t  multMCCharged = 0;
     Int_t  nMCPart = stack->GetNprimary();
     Float_t* eta_gen = new Float_t[nMCPart];  
@@ -614,7 +663,6 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
     Int_t* detectedPrimaryPartLay1 = new Int_t[nMCPart];
     Int_t* detectedPrimaryPartLay2 = new Int_t[nMCPart];
     Int_t* detectablePrimaryPart = new Int_t[nMCPart];
-    Bool_t* stoppingPartAtLay2 = new Bool_t[nMCPart];  
 
     // Loading track references...
     Float_t rminL1 = 3.4;
@@ -622,9 +670,6 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
     Float_t rminL2 = 6.9;
     Float_t rmaxL2 = 7.9;
 
-    Int_t* nTrackRefPerPrim = new Int_t[nMCPart];
-    Float_t trackRefPerPrim[nMCPart][500][2];
-
     TTree* tRefTree = eventHandler->TreeTR();
 
     AliTrackReference *tref=0x0;
@@ -649,94 +694,77 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
       detectedPrimaryPartLay1[multMCCharged]=kFALSE;
       detectedPrimaryPartLay2[multMCCharged]=kFALSE;
       detectablePrimaryPart[multMCCharged]=kFALSE;
-      stoppingPartAtLay2[multMCCharged]=kFALSE;
 
       fHistoPt->Fill(pt_gen[multMCCharged]);
 
       AliMCParticle* mcpart = mcEvent->GetTrack(imc);
       Int_t nref = mcpart->GetNumberOfTrackReferences();
-      nTrackRefPerPrim[multMCCharged]=nref;
-
-      Bool_t nonStoppingPart=kFALSE;
-
-      // Loop over all the track refs of the track
-      for (Int_t iref=0;iref<nref;iref++) { 
-        tref = mcpart->GetTrackReference(iref);
-        if (tref) {
-          trackRefPerPrim[multMCCharged][iref][0]=tref->R();
-          trackRefPerPrim[multMCCharged][iref][1]=(Float_t)tref->DetectorId();
-          if (trackRefPerPrim[multMCCharged][iref][1]==0) fHistTRradius->Fill(trackRefPerPrim[multMCCharged][iref][0]); // ITS only 
-//        Printf("det ID: %f R: %f", trackRefPerPrim[multMCCharged][iref][1], trackRefPerPrim[multMCCharged][iref][0]);
-          if (iref==nref-1&&trackRefPerPrim[multMCCharged][iref][1]!=-1) {
-            nonStoppingPart=kTRUE;
+
+      // Detectable primaries 
+      if (nref==0) {
+        detectablePrimaryPart[multMCCharged]=kTRUE;
+        fHistoDetectableNotr->Fill(eta_gen[multMCCharged],vtxMC[2],pt_gen[multMCCharged]);
+      } else if (nref>0) {
+        tref = mcpart->GetTrackReference(nref-1);
+        if (tref->DetectorId()!=-1) {
+          detectablePrimaryPart[multMCCharged]=kTRUE;
+          fHistoNonStoppingTracks->Fill(eta_gen[multMCCharged],vtxMC[2]); 
+        } else {
+          for (Int_t iref=0;iref<nref;iref++) {
+            tref = mcpart->GetTrackReference(iref);
+            if (tref) {
+              if (tref->R()>rminL2&&tref->R()<rmaxL2) {
+                if (tref->DetectorId()==0) {
+                  detectablePrimaryPart[multMCCharged]=kTRUE;
+                  detectedPrimaryPartLay2[multMCCharged]=kTRUE; 
+                  break;
+                } 
+              } else if (tref->R()>rmaxL2) {
+                  detectablePrimaryPart[multMCCharged]=kTRUE;
+                  fHistoDetectabletr->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  if (tref->DetectorId()==-1) {
+                    fHistoDetectableTRm1->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                    fHistoRTRm1->Fill(tref->R());
+                  }
+                  else if (tref->DetectorId()==0) fHistoDetectableTrITS->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==1) fHistoDetectableTrTPC->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==2) fHistoDetectableTrFRAME->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==3) fHistoDetectableTrTRD->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==4) fHistoDetectableTrTOF->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==5) fHistoDetectableTrMUON->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==6) fHistoDetectableTrHMPID->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==7) fHistoDetectableTrT0->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==8) fHistoDetectableTrEMCAL->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==12) fHistoDetectableTrFMD->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  else if (tref->DetectorId()==14) fHistoDetectableTrVZERO->Fill(eta_gen[multMCCharged],vtxMC[2]);
+                  break;
+              }
+            }
           }
         }
-      }
-      if (nref==0||nonStoppingPart) {
-        detectablePrimaryPart[multMCCharged]=kTRUE; //0. if prim doesn't disappear it is detectable
-        if (nonStoppingPart) fHistoNonStoppingTracks->Fill(eta_gen[multMCCharged],vtxMC[2]);   
-        if (nTrackRefPerPrim[multMCCharged]==0) fHistoDetectableNotr->Fill(eta_gen[multMCCharged],vtxMC[2],pt_gen[multMCCharged]); 
-      }
+        // Find out detected prims on each layer
+        for (Int_t iref=0; iref<nref; iref++) {
+          tref = mcpart->GetTrackReference(iref);
+          if (tref->R()>rminL1&&tref->R()<rmaxL1&&tref->DetectorId()==0) {
+            detectedPrimaryPartLay1[multMCCharged] = kTRUE;
+            fHistoDetectedLay1->Fill(eta_gen[multMCCharged],vtxMC[2]);
+          }
+          if (tref->R()>rminL2&&tref->R()<rmaxL2&&tref->DetectorId()==0) {
+            detectedPrimaryPartLay2[multMCCharged] = kTRUE;
+            fHistoDetectedLay2->Fill(eta_gen[multMCCharged],vtxMC[2]);
+          }
+        }
+      } 
 
       multMCCharged++;
-    } //end of MC particle loop
+    } // End of MC particle loop
 
-    // Check histos
-    if (multSPD==0 && eventTriggered && esdvtx[2]!=0.) {
-      fHistContributorsvsDeVtx->Fill(esdvtx[2]-vtxMC[2],nContrib);
+    if (esdvtx[2]==0.) {
+      fHistContributorsvsDeVtx->Fill(vtxMC[2],nContrib);
     }
-
-    //Event selection
+    // Event selection
     if (eventTriggered && esdvtx[2]!=0. && multSPD!=0) {
 
-      for (Int_t imc=0; imc<multMCCharged; imc++) {
-
-        // Find out detected prims on each layer
-        for (Int_t iref=0; iref<nTrackRefPerPrim[imc]; iref++) {
-          if (trackRefPerPrim[imc][iref][0]>rminL1&&trackRefPerPrim[imc][iref][0]<rmaxL1&&trackRefPerPrim[imc][iref][1]==0) {
-            detectedPrimaryPartLay1[imc] = kTRUE;
-            fHistoDetectedLay1->Fill(eta_gen[imc],vtxMC[2]);
-            break; 
-          }
-        }
-        for (Int_t iref=0; iref<nTrackRefPerPrim[imc]; iref++) {
-          if (trackRefPerPrim[imc][iref][0]>rminL2&&trackRefPerPrim[imc][iref][0]<rmaxL2) {
-            if (trackRefPerPrim[imc][iref][1]==0) {
-              detectedPrimaryPartLay2[imc] = kTRUE;
-              fHistoDetectedLay2->Fill(eta_gen[imc],vtxMC[2]);
-              detectablePrimaryPart[imc] = kTRUE;  //1. if prim is detected it is also detectable
-              break;
-            } 
-            if (trackRefPerPrim[imc][iref][1]==-1) {
-              stoppingPartAtLay2[imc]=kTRUE;
-            }
-          }
-        }
-   
-        // Find out SPD detectable prims 
-        if (!detectablePrimaryPart[imc]) {
-          for (Int_t iref=0; iref<nTrackRefPerPrim[imc]; iref++) {
-            if (trackRefPerPrim[imc][iref][0]>rmaxL2) { 
-              detectablePrimaryPart[imc] = kTRUE; //2. if it has crossed layer 2 it is detectable
-              fHistoDetectabletr->Fill(eta_gen[imc],vtxMC[2]);
-              if (trackRefPerPrim[imc][iref][1]==-1) {
-                fHistoDetectableTRm1->Fill(eta_gen[imc],vtxMC[2]);
-                fHistoRTRm1->Fill(trackRefPerPrim[imc][iref][0]);
-              }
-              else if (trackRefPerPrim[imc][iref][1]==0) fHistoDetectableTR0->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==1) fHistoDetectableTR1->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==2) fHistoDetectableTR2->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==3) fHistoDetectableTR3->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==4) fHistoDetectableTR4->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==5) fHistoDetectableTR5->Fill(eta_gen[imc],vtxMC[2]);
-              else if (trackRefPerPrim[imc][iref][1]==6) fHistoDetectableTR6->Fill(eta_gen[imc],vtxMC[2]);
-              break;
-            }
-          }
-        }
-//        Printf("Primary nr %d", stackIndexOfPrimaryParts[imc]);  
-      }
       for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 
         Int_t labL1 = multESD->GetLabel(itracklet,0);
@@ -770,9 +798,8 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
         fHistMCEtavsZTriggMCvtxEvts->Fill(eta_gen[imc],vtxMC[2]);
       }
 
-      fHistTrigVtxCorrDen->Fill(multSPD,vtxMC[2]);
-
-    } //ev selection
+      if (esdvtx[2]>=-10. && esdvtx[2]<10.) fHistTrigVtxCorrDen->Fill(multSPD,vtxMC[2]);
+    } // End of selected events
  
     if (eventTriggered) {
       fHistTrigCorrDen->Fill(multSPD,vtxMC[2]);
@@ -790,13 +817,11 @@ void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
 
     delete[] eta_gen;
     delete[] pt_gen;
-    delete[] nTrackRefPerPrim;
     delete[] stackIndexOfPrimaryParts;
     delete[] reconstructedPrimaryPart;
     delete[] detectedPrimaryPartLay1;
     delete[] detectedPrimaryPartLay2;
     delete[] detectablePrimaryPart;
-    delete[] stoppingPartAtLay2;
   }
 }      
 
@@ -845,6 +870,7 @@ void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec= dynamic_cast<TH2F*> (fOutput->FindObject("fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec"));
 
   if (fCorr) {
+
     fHistBkgCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistBkgCorrNum"));
     fHistBkgCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistBkgCorrDen"));
 
@@ -857,18 +883,15 @@ void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
 
     fHistTrackTrigCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrackTrigCorrDen")); 
 
-    //Event level
     fHistTrigVtxCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigVtxCorrNum"));
     fHistTrigVtxCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigVtxCorrDen"));
 
     fHistTrigCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigCorrDen"));
 
-    //MC distrib
     fHistMCEtavsZTriggMCvtxEvts = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZTriggMCvtxEvts"));
     fHistMCEtavsZTriggESDvtxEvts = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZTriggESDvtxEvts"));
     fHistMCEtavsZ = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZ"));
 
-    //Check histos  
     fHistTRradius = dynamic_cast<TH1F*> (fOutput->FindObject("fHistTRradius"));
 
     fHistContributorsvsDeVtx = dynamic_cast<TH2F*> (fOutput->FindObject("fHistContributorsvsDeVtx"));
@@ -884,13 +907,18 @@ void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
     fHistoPt = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoPt"));
 
     fHistoDetectableTRm1 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTRm1"));
-    fHistoDetectableTR0 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR0"));
-    fHistoDetectableTR1 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR1"));
-    fHistoDetectableTR2 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR2"));
-    fHistoDetectableTR3 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR3"));
-    fHistoDetectableTR4 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR4"));
-    fHistoDetectableTR5 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR5"));
-    fHistoDetectableTR6 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTR6"));
+    fHistoDetectableTrITS = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrITS"));
+    fHistoDetectableTrTPC = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrTPC"));
+    fHistoDetectableTrFRAME = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrFRAME"));
+    fHistoDetectableTrTRD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrTRD"));
+    fHistoDetectableTrTOF = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrTOF"));
+    fHistoDetectableTrMUON = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrMUON"));
+    fHistoDetectableTrHMPID = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrHMPID"));
+    fHistoDetectableTrT0 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrT0"));
+    fHistoDetectableTrEMCAL = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrEMCAL"));
+    fHistoDetectableTrFMD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrFMD"));
+    fHistoDetectableTrVZERO = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectableTrVZERO"));
+    
     fHistoRTRm1 = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoRTRm1"));
   }
 }