PWG2/EVCHAR/AliAnalysisTaskSPDdNdEta.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 //////////////////////////////////////////////////////////////////////////
  17 // Class AliAnalysisTaskSPDdNdEta                                       //
  18 // Analysis task for dN/dEta reconstruction with the SPD                //
  19 //                                                                      //
  20 // Author:  M. Nicassio (INFN Bari)                                     //
  21 // Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it         //
  22 //////////////////////////////////////////////////////////////////////////
  23
  24 #include "TChain.h"
  25 #include "TTree.h"
  26
  27 #include "TH1F.h"
  28 #include "TH2F.h"
  29 #include "TH3F.h"
  30 #include "TList.h"
  31 #include "TGeoGlobalMagField.h"
  32 #include "AliMagF.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliMultiplicity.h"
  36 #include "AliESDEvent.h"
  37 #include "AliESDInputHandler.h"
  38
  39 #include "AliESDInputHandlerRP.h"
  40 #include "AliESDVZERO.h"
  41 #include "../ITS/AliITSRecPoint.h"
  42 #include "AliCDBPath.h"
  43 #include "AliCDBManager.h"
  44 #include "AliCDBEntry.h"
  45 #include "AliCDBStorage.h"
  46 #include "AliGeomManager.h"
  47
  48 #include "AliMCEventHandler.h"
  49 #include "AliMCEvent.h"
  50 #include "AliStack.h"
  51 #include "AliTrackReference.h"
  52
  53 #include "AliGenHijingEventHeader.h"
  54 #include "AliGenDPMjetEventHeader.h"
  55
  56 #include "AliLog.h"
  57
  58 #include "AliTriggerAnalysis.h"
  59 #include "AliPhysicsSelection.h"
  60 #include "AliESDCentrality.h"
  61 #include "AliTrackletAlg.h"
  62 #include "AliAnalysisTaskSPDdNdEta.h"
  63
  64
  65 ClassImp(AliAnalysisTaskSPDdNdEta)
  66 //________________________________________________________________________
  67 AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
  68   : AliAnalysisTaskSE(name),
  69
  70   fmyESD(0),
  71   fOutput(0),
  72
  73   fUseMC(kFALSE),
  74   fTrigger(AliTriggerAnalysis::kAcceptAll),
  75   fTR(kFALSE),
  76   fRecoTracklets(kFALSE),
  77
  78   fMCCentralityBin(AliAnalysisTaskSPDdNdEta::kall),
  79   fCentrLowLim(0),
  80   fCentrUpLim(0),
  81   fCentrEst(kFALSE),
  82   fMinClMultLay2(0),
  83   fMaxClMultLay2(0),
  84   fMinMultV0(0),
  85   fVtxLim(0),
  86   fPartSpecies(0),
  87
  88   fPhiWindow(0),
  89   fThetaWindow(0),
  90   fPhiShift(0),
  91   fRemoveClustersFromOverlaps(0),
  92   fPhiOverlapCut(0),
  93   fZetaOverlapCut(0),
  94   fPhiRotationAngle(0),
  95   fPhiWindowAna(0),
  96
  97   fMultReco(0),
  98
  99   fV0Ampl(0),
 100   fHistSPDRAWMultvsZ(0),
 101   fHistSPDRAWMultvsZTriggCentrEvts(0),
 102   fHistSPDRAWMultvsZCentrEvts(0),
 103   fHistSPDRAWEtavsZ(0),
 104
 105   fHistSPDmultEtacut(0),
 106   fHistSPDmult(0),
 107   fHistSPDmultcl1(0),
 108   fHistSPDmultcl2(0),
 109   fHistSPDmultcl1vscl2(0),
 110   fHistSPDmultvscl1(0),
 111   fHistSPDmultvscl2(0),
 112
 113   fHistSPDeta(0),
 114   fHistSPDphi(0),
 115   fHistSPDtheta(0),
 116   fHistSPDdePhi(0),
 117   fHistSPDphivsSPDeta(0),
 118   fHistSPDdeTheta(0),
 119   fHistSPDvtx(0),
 120   fHistSPDvtxAnalysis(0),
 121   fHistSPDdePhideTheta(0),
 122
 123   fHistSPDphicl1(0),
 124   fHistSPDphicl2(0),
 125
 126   fHistBkgCorrDen(0),
 127   fHistReconstructedProtons(0),
 128   fHistReconstructedKaons(0),
 129   fHistReconstructedPions(0),
 130   fHistReconstructedOthers(0),
 131   fHistReconstructedSec(0),
 132   fHistBkgCorrDenPrimGen(0),
 133   fHistBkgCombLabels(0),
 134   fHistBkgCorrNum(0),
 135   fHistAlgEffNum(0),
 136   fHistNonDetectableCorrDen(0),
 137
 138   fHistNonDetectableCorrNum(0),
 139   fHistProtons(0),
 140   fHistKaons(0),
 141   fHistPions(0),
 142   fHistOthers(0),
 143   fHistAllPrimaries(0),
 144   fHistTrackCentrEvts(0),
 145   fHistTrackTrigCentrEvts(0),
 146
 147   fHistAllEvts(0),
 148   fHistCentrEvts(0),
 149   fHistTrigCentrEvts(0),
 150   fHistSelEvts(0),
 151
 152   fHistMCmultEtacut(0),
 153   fHistMCmultEtacutvsSPDmultEtacut(0),
 154   fHistMCmultEtacutvsSPDmultcl1(0),
 155   fHistMCmultEtacutvsSPDmultcl2(0),
 156
 157   fHistMCvtxx(0),
 158   fHistMCvtxy(0),
 159   fHistMCvtxz(0),
 160
 161   fHistRecvsGenImpactPar(0),
 162   fHistMCNpart(0),
 163
 164   fHistdPhidThetaPP(0),
 165   fHistdPhidThetaSS(0),
 166   fHistdPhidThetaComb(0),
 167
 168   fHistDeVtx(0)
 169
 170 {
 171
 172   // Constructor
 173
 174   // Define input and output slots here
 175   // Input slot #0 works with a TChain
 176 //  DefineInput(0, TChain::Class());
 177   // Output slot #0 id reserved by the base class for AOD
 178   DefineOutput(1, TList::Class());
 179
 180 }
 181 //________________________________________________________________________
 182 AliAnalysisTaskSPDdNdEta::~AliAnalysisTaskSPDdNdEta()
 183 {
 184
 185   // Destructor
 186
 187   // histograms are in the output list and deleted when the output
 188   // list is deleted by the TSelector dtor
 189   delete fMultReco;
 190
 191   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 192     delete fOutput;
 193     fOutput = 0;
 194   }
 195
 196 }
 197 //________________________________________________________________________
 198 void AliAnalysisTaskSPDdNdEta::UserCreateOutputObjects()
 199 {
 200   if (fRecoTracklets) {
 201     AliCDBManager *man = AliCDBManager::Instance();
 202     man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
 203     if (fUseMC) man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
 204     else man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB");
 205     man->SetRun(137161);
 206     AliCDBEntry* obj = man->Get(AliCDBPath("GRP", "Geometry", "Data"));
 207     AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 208     AliGeomManager::GetNalignable("ITS");
 209     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 210   }
 211
 212   fMultReco = new AliTrackletAlg();
 213
 214   // Create histograms
 215   fOutput = new TList();
 216   fOutput->SetOwner();
 217
 218   Int_t nBinVtx = 40;
 219   Double_t MaxVtx = 20.;
 220
 221   Int_t nBinMultCorr = 200;
 222   Float_t nMaxMult = 20000.;
 223   Double_t binLimMultCorr[nBinMultCorr+1];
 224   binLimMultCorr[0] = 0.;
 225   binLimMultCorr[1] = 1.;
 226   for (Int_t i = 2; i<=nBinMultCorr;++i) {
 227     binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
 228   }
 229
 230   Int_t nBinEtaCorr = 60;
 231   Float_t etaMax = 3.;
 232   Float_t etaMin = -3.;
 233   Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
 234   for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
 235     binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
 236   }
 237
 238   // Data to be corrected
 239   // ...event level
 240
 241   fV0Ampl = new TH1F("fV0Ampl","",500,0.,30000);
 242   fOutput->Add(fV0Ampl);
 243
 244   fHistSPDRAWMultvsZ = new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 245   fHistSPDRAWMultvsZ->GetXaxis()->SetTitle("Tracklet multiplicity");
 246   fHistSPDRAWMultvsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 247   fOutput->Add(fHistSPDRAWMultvsZ);
 248
 249   fHistSPDRAWMultvsZTriggCentrEvts = new TH2F("fHistSPDRAWMultvsZTriggCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 250   fHistSPDRAWMultvsZTriggCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
 251   fHistSPDRAWMultvsZTriggCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 252   fOutput->Add(fHistSPDRAWMultvsZTriggCentrEvts);
 253
 254   fHistSPDRAWMultvsZCentrEvts = new TH2F("fHistSPDRAWMultvsZCentrEvts", "",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 255   fHistSPDRAWMultvsZCentrEvts->GetXaxis()->SetTitle("Tracklet multiplicity");
 256   fHistSPDRAWMultvsZCentrEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 257   fOutput->Add(fHistSPDRAWMultvsZCentrEvts);
 258
 259   // ...track level
 260   fHistSPDRAWEtavsZ = new TH2F("fHistSPDRAWEtavsZ", "Tracklet pseudorapidity distribution", nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 261   fHistSPDRAWEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
 262   fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 263   fOutput->Add(fHistSPDRAWEtavsZ);
 264
 265   fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution",nBinMultCorr,binLimMultCorr);
 266   fHistSPDmultEtacut->GetXaxis()->SetTitle("Tracklet multiplicity (|#eta|<1.4)");
 267   fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
 268   fOutput->Add(fHistSPDmultEtacut);
 269
 270   fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", nBinMultCorr,binLimMultCorr);
 271   fHistSPDmult->GetXaxis()->SetTitle("Tracklet multiplicity");
 272   fHistSPDmult->GetYaxis()->SetTitle("Entries");
 273   fOutput->Add(fHistSPDmult);
 274
 275   fHistSPDmultcl1 = new TH1F("fHistSPDmultcl1", "Inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
 276   fHistSPDmultcl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
 277   fHistSPDmultcl1->GetYaxis()->SetTitle("Entries");
 278   fOutput->Add(fHistSPDmultcl1);
 279
 280   fHistSPDmultcl2 = new TH1F("fHistSPDmultcl2", "Outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr);
 281   fHistSPDmultcl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
 282   fHistSPDmultcl2->GetYaxis()->SetTitle("Entries");
 283   fOutput->Add(fHistSPDmultcl2);
 284
 285   fHistSPDmultcl1vscl2 = new TH2F("fHistSPDmultcl1vscl2", "Inner layer cluster vs outer multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
 286   fHistSPDmultcl1vscl2->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
 287   fHistSPDmultcl1vscl2->GetYaxis()->SetTitle("Outer layer cluster multiplicity");
 288   fOutput->Add(fHistSPDmultcl1vscl2);
 289
 290   fHistSPDmultvscl1 = new TH2F("fHistSPDmultvscl1", "Tracklet vs inner layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
 291   fHistSPDmultvscl1->GetXaxis()->SetTitle("Inner layer cluster multiplicity");
 292   fHistSPDmultvscl1->GetYaxis()->SetTitle("Tracklet multiplicity");
 293   fOutput->Add(fHistSPDmultvscl1);
 294
 295   fHistSPDmultvscl2 = new TH2F("fHistSPDmultvscl2", "Tracklet vs outer layer cluster multiplicity", nBinMultCorr,binLimMultCorr, nBinMultCorr,binLimMultCorr);
 296   fHistSPDmultvscl2->GetXaxis()->SetTitle("Outer layer cluster multiplicity");
 297   fHistSPDmultvscl2->GetYaxis()->SetTitle("Tracklet multiplicity");
 298   fOutput->Add(fHistSPDmultvscl2);
 299
 300   fHistSPDphi = new TH1F("fHistSPDphi", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
 301   fHistSPDphi->GetXaxis()->SetTitle("#varphi [rad]");
 302   fHistSPDphi->GetYaxis()->SetTitle("Entries");
 303   fOutput->Add(fHistSPDphi);
 304
 305   fHistSPDtheta = new TH1F("fHistSPDtheta", "Tracklet #theta  distribution", 180, 0.,TMath::Pi());
 306   fHistSPDtheta->GetXaxis()->SetTitle("#theta [rad]");
 307   fHistSPDtheta->GetYaxis()->SetTitle("Entries");
 308   fOutput->Add(fHistSPDtheta);
 309
 310   fHistSPDdePhi= new TH1F("fHistSPDdePhi", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
 311   fHistSPDdePhi->GetXaxis()->SetTitle("#Delta#varphi [rad]");
 312   fHistSPDdePhi->GetYaxis()->SetTitle("Entries");
 313   fOutput->Add(fHistSPDdePhi);
 314
 315   fHistSPDphivsSPDeta= new TH2F("fHistSPDphivsSPDeta", "Tracklets - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());
 316   fHistSPDphivsSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 317   fHistSPDphivsSPDeta->GetYaxis()->SetTitle("#varphi [rad]");
 318   fOutput->Add(fHistSPDphivsSPDeta);
 319
 320   fHistSPDdeTheta= new TH1F("fHistSPDdeTheta", "Tracklet #Delta#theta distribution",100,-0.05,.05);
 321   fHistSPDdeTheta->GetXaxis()->SetTitle("#Delta#theta [rad]");
 322   fHistSPDdeTheta->GetYaxis()->SetTitle("Entries");
 323   fOutput->Add(fHistSPDdeTheta);
 324
 325   fHistSPDvtx = new TH1F("fHistSPDvtx", "SPD vertex distribution",nBinVtx,-MaxVtx,MaxVtx);
 326   fHistSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 327   fHistSPDvtx->GetYaxis()->SetTitle("Entries");
 328   fOutput->Add(fHistSPDvtx);
 329
 330   fHistSPDvtxAnalysis = new TH1F("fHistSPDvtxAnalysis", "SPD vertex distribution: events selected for the analysis",nBinVtx,-MaxVtx,MaxVtx);
 331   fHistSPDvtxAnalysis->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 332   fHistSPDvtxAnalysis->GetYaxis()->SetTitle("Entries");
 333   fOutput->Add(fHistSPDvtxAnalysis);
 334
 335   fHistSPDdePhideTheta= new TH2F("fHistSPDdePhideTheta", "Tracklet #Delta#varphi  distribution",2000,-1.,1.,1000,-0.25,.25);
 336   fHistSPDdePhideTheta->GetXaxis()->SetTitle("#Delta#varphi [rad]");
 337   fHistSPDdePhideTheta->GetYaxis()->SetTitle("#Delta#theta [rad]");
 338   fOutput->Add(fHistSPDdePhideTheta);
 339
 340   fHistSPDphicl1 = new TH1F("fHistSPDphicl1", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
 341   fHistSPDphicl1->GetXaxis()->SetTitle("#varphi [rad]");
 342   fHistSPDphicl1->GetYaxis()->SetTitle("Entries");
 343   fOutput->Add(fHistSPDphicl1);
 344
 345   fHistSPDphicl2 = new TH1F("fHistSPDphicl2", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
 346   fHistSPDphicl2->GetXaxis()->SetTitle("#varphi [rad]");
 347   fHistSPDphicl2->GetYaxis()->SetTitle("Entries");
 348   fOutput->Add(fHistSPDphicl2);
 349
 350
 351   if (fUseMC) {
 352
 353     // Track level correction histograms
 354     fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 355     fOutput->Add(fHistBkgCorrDen);
 356
 357     fHistReconstructedProtons = new TH2F("fHistReconstructedProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 358     fOutput->Add(fHistReconstructedProtons);
 359     fHistReconstructedKaons = new TH2F("fHistReconstructedKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 360     fOutput->Add(fHistReconstructedKaons);
 361     fHistReconstructedPions = new TH2F("fHistReconstructedPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 362     fOutput->Add(fHistReconstructedPions);
 363     fHistReconstructedOthers = new TH2F("fHistReconstructedOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 364     fOutput->Add(fHistReconstructedOthers);
 365
 366     fHistBkgCorrDenPrimGen = new TH2F("fHistBkgCorrDenPrimGen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 367     fOutput->Add(fHistBkgCorrDenPrimGen);
 368
 369     fHistBkgCombLabels = new TH2F("fHistBkgCombLabels","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 370     fOutput->Add(fHistBkgCombLabels);
 371
 372     if (fTR) {
 373       fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 374       fOutput->Add(fHistBkgCorrNum);
 375
 376       fHistAlgEffNum = new TH2F("fHistAlgEffNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 377       fOutput->Add(fHistAlgEffNum);
 378
 379       fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 380       fOutput->Add(fHistNonDetectableCorrDen);
 381
 382     }
 383
 384     fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 385     fOutput->Add(fHistNonDetectableCorrNum);
 386
 387     fHistProtons = new TH2F("fHistProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 388     fOutput->Add(fHistProtons);
 389     fHistKaons = new TH2F("fHistKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 390     fOutput->Add(fHistKaons);
 391     fHistPions = new TH2F("fHistPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 392     fOutput->Add(fHistPions);
 393     fHistOthers = new TH2F("fHistOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 394     fOutput->Add(fHistOthers);
 395     fHistReconstructedSec = new TH2F("fHistReconstructedSec","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 396     fOutput->Add(fHistReconstructedSec);
 397
 398
 399     fHistAllPrimaries = new TH2F("fHistAllPrimaries","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 400     fOutput->Add(fHistAllPrimaries);
 401
 402     fHistTrackCentrEvts = new TH2F("fHistTrackCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 403     fOutput->Add(fHistTrackCentrEvts);
 404
 405     fHistTrackTrigCentrEvts = new TH2F("fHistTrackTrigCentrEvts","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 406     fOutput->Add(fHistTrackTrigCentrEvts);
 407
 408
 409     // Event level correction histograms
 410     fHistAllEvts = new TH2F("fHistAllEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 411     fOutput->Add(fHistAllEvts);
 412
 413     fHistCentrEvts = new TH2F("fHistCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 414     fOutput->Add(fHistCentrEvts);
 415
 416     fHistTrigCentrEvts = new TH2F("fHistTrigCentrEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 417     fOutput->Add(fHistTrigCentrEvts);
 418
 419     fHistSelEvts = new TH2F("fHistSelEvts","",nBinMultCorr,binLimMultCorr,nBinVtx,-MaxVtx,MaxVtx);
 420     fOutput->Add(fHistSelEvts);
 421
 422
 423     // MC distributions
 424     fHistMCmultEtacut = new TH1F("fHistMCmultEtacut","Generated multiplicity",nBinMultCorr,binLimMultCorr);
 425     fHistMCmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
 426     fHistMCmultEtacut->GetYaxis()->SetTitle("Entries");
 427     fOutput->Add(fHistMCmultEtacut);
 428
 429     fHistMCmultEtacutvsSPDmultEtacut = new TH2F("fHistMCmultEtacutvsSPDmultEtacut","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
 430     fHistMCmultEtacutvsSPDmultEtacut->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
 431     fHistMCmultEtacutvsSPDmultEtacut->GetYaxis()->SetTitle("Tracklet multiplicity |#eta|<1.4");
 432     fOutput->Add(fHistMCmultEtacutvsSPDmultEtacut);
 433
 434     fHistMCmultEtacutvsSPDmultcl1 = new TH2F("fHistMCmultEtacutvsSPDmultcl1","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
 435     fHistMCmultEtacutvsSPDmultcl1->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
 436     fHistMCmultEtacutvsSPDmultcl1->GetYaxis()->SetTitle("Cluster inner layer multiplicity");
 437     fOutput->Add(fHistMCmultEtacutvsSPDmultcl1);
 438
 439     fHistMCmultEtacutvsSPDmultcl2 = new TH2F("fHistMCmultEtacutvsSPDmultcl2","",nBinMultCorr,binLimMultCorr,nBinMultCorr,binLimMultCorr);
 440     fHistMCmultEtacutvsSPDmultcl2->GetXaxis()->SetTitle("Generated multiplicity |#eta|<1.4");
 441     fHistMCmultEtacutvsSPDmultcl2->GetYaxis()->SetTitle("Cluster outer layer multiplicity");
 442     fOutput->Add(fHistMCmultEtacutvsSPDmultcl2);
 443
 444
 445     fHistMCvtxx = new TH1F("fHistMCvtxx", "MC vertex distribution - x",100,-.5,.5);
 446     fOutput->Add(fHistMCvtxx);
 447     fHistMCvtxy = new TH1F("fHistMCvtxy", "MC vertex distribution - y",100,-.5,.5);
 448     fOutput->Add(fHistMCvtxy);
 449     fHistMCvtxz = new TH1F("fHistMCvtxz", "MC vertex distribution - z",500,-50.,50.);
 450     fOutput->Add(fHistMCvtxz);
 451
 452     fHistRecvsGenImpactPar= new TH2F("fHistRecvsGenImpactPar","",20,0.,20.,20,0.,20.);
 453     fHistRecvsGenImpactPar->GetXaxis()->SetTitle("b_{gen} [fm]");
 454     fHistRecvsGenImpactPar->GetYaxis()->SetTitle("b_{rec} [fm]");
 455     fOutput->Add(fHistRecvsGenImpactPar);
 456
 457     fHistMCNpart = new TH1F("fHistMCNpart","Number of participants",450,0,450);
 458     fHistMCNpart->GetXaxis()->SetTitle("N_{part}");
 459     fHistMCNpart->GetYaxis()->SetTitle("Entries");
 460     fOutput->Add(fHistMCNpart);
 461
 462     fHistdPhidThetaPP = new TH2F("fHistdPhidThetaPP","",2000,-1.,1.,1000,-0.25,.25);
 463     fOutput->Add(fHistdPhidThetaPP);
 464     fHistdPhidThetaSS = new TH2F("fHistdPhidThetaSS","",2000,-1.,1.,1000,-0.25,.25);
 465     fOutput->Add(fHistdPhidThetaSS);
 466     fHistdPhidThetaComb = new TH2F("fHistdPhidThetaComb","",2000,-1.,1.,1000,-0.25,.25);
 467     fOutput->Add(fHistdPhidThetaComb);
 468
 469     fHistDeVtx = new TH2F("fHistDeVtx","",80,-20.,20.,5000,-0.5,0.5);
 470     fHistDeVtx->GetXaxis()->SetTitle("z_{MCvtx} [cm]");
 471     fHistDeVtx->GetYaxis()->SetTitle("z_{MCvtx}-z_{SPDvtx} [cm]");
 472     fOutput->Add(fHistDeVtx);
 473
 474   }
 475   PostData(1, fOutput);
 476 //  Printf("Output objects created...");
 477 }
 478
 479 //________________________________________________________________________
 480 void AliAnalysisTaskSPDdNdEta::UserExec(Option_t *)
 481 {
 482   // Main loop
 483
 484   // Called for each event
 485 //  Printf("User exec..........");
 486
 487   AliESDInputHandlerRP *hand = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 488
 489   fmyESD = dynamic_cast<AliESDEvent*>(InputEvent());
 490   if (!fmyESD) {
 491     Printf("ERROR: fESD not available");
 492     return;
 493   }
 494
 495
 496   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 497   if (!field && !fmyESD->InitMagneticField()) {Printf("Failed to initialize the B field\n");return;}
 498
 499
 500   // Trigger selection
 501   Bool_t eventTriggered = kTRUE;
 502   static AliTriggerAnalysis* triggerAnalysis = 0;
 503   if (eventTriggered) // applying an offline trigger
 504     eventTriggered = triggerAnalysis->IsTriggerFired(fmyESD, fTrigger);
 505   if (!eventTriggered)
 506     Printf("No trigger");
 507
 508
 509   // Centrality selection
 510   Bool_t eventInCentralityBin = kFALSE;
 511 /*  AliESDCentrality *centrality = fmyESD->GetCentrality();
 512   if (fCentrEst=="") eventInCentralityBin = kTRUE;
 513   else {
 514     if(!centrality) {
 515       AliError("Centrality object not available");
 516     }  else {
 517       if (centrality->IsEventInCentralityClass(fCentrLowLim,fCentrUpLim,fCentrEst.Data())) eventInCentralityBin = kTRUE;
 518     }
 519   }
 520 */
 521   if (fCentrEst) {
 522     AliESDVZERO* esdV0 = fmyESD->GetVZEROData();
 523     Float_t multV0A=esdV0->GetMTotV0A();
 524     Float_t multV0C=esdV0->GetMTotV0C();
 525     fV0Ampl->Fill(multV0A+multV0C);
 526     if (multV0A+multV0C>=fMinMultV0) eventInCentralityBin = kTRUE;
 527   } else if (!fCentrEst) {
 528     eventInCentralityBin = kTRUE;
 529   }
 530
 531   const AliMultiplicity* multESD = fmyESD->GetMultiplicity();
 532
 533   // ESD vertex
 534   Bool_t eventWithVertex = kFALSE;
 535   const AliESDVertex* vtxESD = fmyESD->GetVertex();
 536   const AliESDVertex* vtxTPC = fmyESD->GetPrimaryVertexTPC();
 537   Double_t esdvtx[3];
 538   vtxESD->GetXYZ(esdvtx);
 539   Int_t nContrib = vtxESD->GetNContributors();
 540   Int_t nContribTPC = vtxTPC->GetNContributors();
 541   if (nContrib>0&&nContribTPC>0) {
 542     if (vtxESD->GetDispersion()<0.04)
 543       if (vtxESD->GetZRes()<0.25)
 544         if (nContribTPC>(-10.+0.25*multESD->GetNumberOfITSClusters(0))) {
 545           fHistSPDvtx->Fill(esdvtx[2]);
 546           if (TMath::Abs(esdvtx[2])<fVtxLim)
 547             eventWithVertex = kTRUE;
 548         }
 549   }
 550
 551   // Reconstructing or loading tracklets...
 552   Int_t multSPD = multESD->GetNumberOfTracklets();
 553   Int_t nSingleCl1 = multESD->GetNumberOfSingleClusters();
 554   Int_t multSPDcl1 = nSingleCl1 + multSPD;
 555   Int_t multSPDEtacut = 0;
 556   Int_t multSPDAna = 0;
 557
 558   Int_t multSPDcl2 = multESD->GetNumberOfITSClusters(1);
 559
 560   Float_t trackletCoord[multSPDcl1][5];
 561   Float_t trackletLab[multSPDcl1][2];
 562
 563   Bool_t eventSelected = kFALSE;
 564   // Event selection: in centrality bin, triggered with vertex
 565   if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
 566     eventSelected = kTRUE;
 567     fHistSPDmultcl1->Fill(multSPDcl1);
 568     fHistSPDmultcl2->Fill(multSPDcl2);
 569     fHistSPDmultcl1vscl2->Fill(multSPDcl1,multSPDcl2);
 570     fHistSPDmultvscl1->Fill(multSPDcl1,multSPD);
 571     fHistSPDmultvscl2->Fill(multSPDcl2,multSPD);
 572     if (fRecoTracklets) {
 573       // Load ITS rec points tree
 574       TTree *itsClusterTree = hand->GetTreeR("ITS");
 575       if (!itsClusterTree) {
 576         AliError(" Invalid ITS cluster tree !\n");
 577         return;
 578       }
 579       float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()};
 580
 581       fMultReco->SetPhiWindow(fPhiWindow);
 582       fMultReco->SetThetaWindow(fThetaWindow);
 583       fMultReco->SetPhiShift(fPhiShift);
 584       fMultReco->SetRemoveClustersFromOverlaps(fRemoveClustersFromOverlaps);
 585       fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
 586       fMultReco->SetZetaOverlapCut(fZetaOverlapCut);
 587       fMultReco->SetPhiRotationAngle(fPhiRotationAngle);
 588       fMultReco->SetHistOn(kFALSE);
 589
 590       fMultReco->Reconstruct(itsClusterTree,vtxf,vtxf);
 591 //      Printf("cl 1 from alg %d",fMultReco->GetNClustersLayer1());
 592 //      Printf("cl 2 from alg %d",fMultReco->GetNClustersLayer2());
 593       multSPD = fMultReco->GetNTracklets();
 594 //      Printf("tracklets found...%d",multSPD);
 595       nSingleCl1 = fMultReco->GetNSingleClusters();
 596 //      Printf("singl found...%d",nSingleCl1);
 597       for (Int_t itr = 0; itr<multSPD;++itr) {
 598         trackletCoord[itr][3] = fMultReco->GetTracklet(itr)[0]; //theta
 599         trackletCoord[itr][2] = fMultReco->GetTracklet(itr)[1]; //phi
 600         trackletCoord[itr][1] = fMultReco->GetTracklet(itr)[3]; //deTheta
 601         trackletCoord[itr][0] = fMultReco->GetTracklet(itr)[2]; //dephi
 602         trackletCoord[itr][4] = -TMath::Log(TMath::Tan(trackletCoord[itr][3]/2.)); //eta
 603         trackletLab[itr][0] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[4]);  //label lay1
 604         trackletLab[itr][1] = static_cast<Int_t>(fMultReco->GetTracklet(itr)[5]);  //label lay2
 605       }
 606       for (Int_t icl1 = 0; icl1<multSPDcl1;++icl1) {
 607         fHistSPDphicl1->Fill(fMultReco->GetClusterLayer1(icl1)[1]);
 608       }
 609       for (Int_t icl2 = 0; icl2<multSPDcl2; ++icl2) {
 610         fHistSPDphicl2->Fill(fMultReco->GetClusterLayer2(icl2)[1]);
 611       }
 612
 613     } else {
 614       // set variables from ESD
 615       for (Int_t itr = 0; itr<multSPD;++itr) {
 616         trackletCoord[itr][3] = multESD->GetTheta(itr);       //theta
 617         trackletCoord[itr][2] = multESD->GetPhi(itr);         //phi
 618         trackletCoord[itr][1] = multESD->GetDeltaTheta(itr);  //deTheta
 619         trackletCoord[itr][0] = multESD->GetDeltaPhi(itr);    //dePhi
 620         trackletCoord[itr][4] = multESD->GetEta(itr);         //eta
 621         trackletLab[itr][0] = multESD->GetLabel(itr,0);       //label lay1
 622         trackletLab[itr][1] = multESD->GetLabel(itr,1);       //label lay2
 623       }
 624     }
 625 //    Printf("tracklets in ESD...%d",multESD->GetNumberOfTracklets());
 626
 627     for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 628       if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) { // select tracklets with tighter cuts
 629         fHistSPDRAWEtavsZ->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 630         fHistSPDphi->Fill(trackletCoord[itracklet][2]);
 631         fHistSPDtheta->Fill(trackletCoord[itracklet][3]);
 632         fHistSPDdePhi->Fill(trackletCoord[itracklet][0]);
 633         fHistSPDdeTheta->Fill(trackletCoord[itracklet][1]);
 634
 635         fHistSPDphivsSPDeta->Fill(trackletCoord[itracklet][4],trackletCoord[itracklet][2]);
 636         multSPDAna++;
 637         // Calculate multiplicity in etacut
 638         if (TMath::Abs(trackletCoord[itracklet][4])<1.4) multSPDEtacut++;
 639       }
 640       // for combinatorial background normalization
 641       fHistSPDdePhideTheta->Fill(trackletCoord[itracklet][0], trackletCoord[itracklet][1]);
 642
 643     }
 644     if (multSPDAna!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
 645     fHistSPDmultEtacut->Fill(multSPDEtacut);
 646     fHistSPDmult->Fill(multSPDAna);
 647
 648     fHistSPDRAWMultvsZ->Fill(multSPDAna,esdvtx[2]);
 649
 650   } // End selected events
 651
 652   if (eventInCentralityBin) {
 653     fHistSPDRAWMultvsZCentrEvts->Fill(multSPDAna,esdvtx[2]);
 654     if (eventTriggered) fHistSPDRAWMultvsZTriggCentrEvts->Fill(multSPDAna,esdvtx[2]);
 655   }
 656
 657   if (fUseMC) {
 658
 659     if (!fMCEvent) {
 660       AliError("No MC info found");
 661       return;
 662     }
 663     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 664     if (!eventHandler) {
 665       Printf("ERROR: Could not retrieve MC event handler");
 666       return;
 667     }
 668
 669     AliStack* stack = fMCEvent->Stack();
 670     if (!stack) {
 671       AliDebug(AliLog::kError, "Stack not available");
 672       return;
 673     }
 674
 675     Float_t impactParameter = 0.;
 676     Int_t npart = 0;
 677
 678     AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(fMCEvent->Header()->GenEventHeader());
 679     AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(fMCEvent->Header()->GenEventHeader());
 680
 681     if (hijingHeader) impactParameter = hijingHeader->ImpactParameter();
 682     else if (dpmHeader) impactParameter = dpmHeader->ImpactParameter();
 683
 684     Bool_t IsEventInMCCentralityBin = kFALSE;
 685     switch (fMCCentralityBin) {
 686
 687       case 1: if (impactParameter<3) IsEventInMCCentralityBin = kTRUE;
 688       break;
 689       case 2: if (impactParameter>3&&impactParameter<6) IsEventInMCCentralityBin = kTRUE;
 690       break;
 691       case 3: if (impactParameter>6&&impactParameter<9) IsEventInMCCentralityBin = kTRUE;
 692       break;
 693       case 4: if (impactParameter>9&&impactParameter<12) IsEventInMCCentralityBin = kTRUE;
 694       break;
 695       case 5: if (impactParameter>12&&impactParameter<15) IsEventInMCCentralityBin = kTRUE;
 696       break;
 697       case 6: if (impactParameter>15) IsEventInMCCentralityBin = kTRUE;
 698       break;
 699       case 7:  IsEventInMCCentralityBin = kTRUE;
 700       break;
 701
 702     }
 703
 704     if (IsEventInMCCentralityBin) {
 705       // MC vertex
 706       TArrayF vtxMC(3);
 707       fMCEvent->GenEventHeader()->PrimaryVertex(vtxMC);
 708       fHistMCvtxx->Fill(vtxMC[0]);
 709       fHistMCvtxy->Fill(vtxMC[1]);
 710       fHistMCvtxz->Fill(vtxMC[2]);
 711
 712 //      Printf("Impact parameter gen: %f", impactParameter);
 713        if (hijingHeader) npart = hijingHeader->TargetParticipants()+hijingHeader->ProjectileParticipants();
 714        else if (dpmHeader)npart = dpmHeader->TargetParticipants()+dpmHeader->ProjectileParticipants();
 715
 716       //Rec centrality vs gen centrality
 717       AliESDZDC *zdcRec = fmyESD->GetESDZDC();
 718       Double_t impactParameterZDC = zdcRec->GetImpactParameter();
 719 //      Printf("Impact parameter rec: %f", impactParameterZDC);
 720
 721       fHistRecvsGenImpactPar->Fill(impactParameter,impactParameterZDC);
 722       fHistMCNpart->Fill(npart);
 723
 724       // Tracks from MC
 725       Int_t  multMCCharged = 0;
 726       Int_t  multMCChargedEtacut = 0;
 727 //      Int_t  nMCPart = stack->GetNprimary();
 728       Int_t  nMCPart = stack->GetNtrack();  // decay products of D and B mesons are also primaries and produced in HIJING during transport
 729       Float_t* etagen = new Float_t[nMCPart];
 730       Int_t* stackIndexOfPrimaryParts = new Int_t[nMCPart];
 731       Bool_t* reconstructedPrimaryPart = new Bool_t[nMCPart];
 732       Bool_t* detectedPrimaryPartLay1 = new Bool_t[nMCPart];
 733       Bool_t* detectedPrimaryPartLay2 = new Bool_t[nMCPart];
 734       Bool_t* detectablePrimaryPart = new Bool_t[nMCPart];
 735
 736       TTree* tRefTree;
 737       if (fTR) {
 738         tRefTree = eventHandler->TreeTR();
 739         fMCEvent->ConnectTreeTR(tRefTree);
 740       }
 741
 742       // Loop over MC particles
 743       for (Int_t imc=0; imc<nMCPart; imc++) {
 744         AliMCParticle *mcpart  = (AliMCParticle*)fMCEvent->GetTrack(imc);
 745
 746         Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
 747         if (!isPrimary)            continue;
 748         if (mcpart->Charge() == 0) continue;
 749         Float_t theta = mcpart->Theta();
 750         if (theta==0 || theta==TMath::Pi())    continue;
 751         Float_t eta = mcpart->Eta();
 752 //        Float_t pt = mcpart->Pt();
 753         etagen[multMCCharged] = eta;
 754         stackIndexOfPrimaryParts[multMCCharged] = imc;
 755
 756         reconstructedPrimaryPart[multMCCharged]=kFALSE;
 757         detectedPrimaryPartLay1[multMCCharged]=kFALSE;
 758         detectedPrimaryPartLay2[multMCCharged]=kFALSE;
 759         detectablePrimaryPart[multMCCharged]=kFALSE;
 760
 761         if (fTR) {
 762           Int_t nref = mcpart->GetNumberOfTrackReferences();
 763           if (nref==0) {
 764             detectablePrimaryPart[multMCCharged]= kTRUE;
 765           } else if (nref>0) {
 766             // Check if the primary is detectable
 767             detectablePrimaryPart[multMCCharged] = IsDetectablePrimary(nref,mcpart);
 768             // Check if the primary is detected (if SPD 100% efficient)
 769             if (detectablePrimaryPart[multMCCharged]) {
 770               detectedPrimaryPartLay1[multMCCharged] = IsDetectedPrimary(nref,mcpart,0);
 771               detectedPrimaryPartLay2[multMCCharged] = IsDetectedPrimary(nref,mcpart,1);
 772             }
 773           }
 774         }
 775         if (eventSelected&&fPartSpecies) {
 776           if (TMath::Abs(mcpart->PdgCode())==2212) fHistProtons->Fill(etagen[multMCCharged],vtxMC[2]);
 777           else if (TMath::Abs(mcpart->PdgCode())==321) fHistKaons->Fill(etagen[multMCCharged],vtxMC[2]);
 778           else if (TMath::Abs(mcpart->PdgCode())==211) fHistPions->Fill(etagen[multMCCharged],vtxMC[2]);
 779           else fHistOthers->Fill(etagen[multMCCharged],vtxMC[2]);  //includes leptons pdg->11,13
 780         }
 781         multMCCharged++;
 782         if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
 783       } // end of MC particle loop
 784
 785       fHistMCmultEtacut->Fill(multMCChargedEtacut);
 786
 787       // Event selection: in centrality bin, triggered with vertex
 788       if (eventTriggered&&eventWithVertex&&eventInCentralityBin&&multSPDcl2>fMinClMultLay2&&multSPDcl2<fMaxClMultLay2) {
 789
 790         fHistDeVtx->Fill(vtxMC[2],vtxMC[2]-esdvtx[2]);
 791
 792         for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 793           if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna)
 794             fHistBkgCorrDen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 795
 796           if (trackletLab[itracklet][0]==trackletLab[itracklet][1]) {
 797             Bool_t trakletByPrim = kFALSE;
 798             for (Int_t imc=0; imc<multMCCharged; imc++) {
 799               if (trackletLab[itracklet][0]==stackIndexOfPrimaryParts[imc]) {
 800                 fHistdPhidThetaPP->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
 801                 if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
 802                   fHistBkgCorrDenPrimGen->Fill(etagen[imc],vtxMC[2]);
 803                   if (fPartSpecies) {
 804                     if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==2212)
 805                       fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 806                     else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==321)
 807                       fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 808                     else if (TMath::Abs(((AliMCParticle*)fMCEvent->GetTrack(stackIndexOfPrimaryParts[imc]))->PdgCode())==211)
 809                       fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 810                     else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 811                   }
 812                 }
 813                 trakletByPrim = kTRUE;
 814                 if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
 815                   reconstructedPrimaryPart[imc]=kTRUE; // tracklet by prim and tr (=rp) on both layers
 816                 }
 817                 break;
 818               }
 819             }
 820             if (!trakletByPrim) {
 821               fHistdPhidThetaSS->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
 822               if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
 823
 824                 fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 825                 if (fPartSpecies) {
 826                   Int_t motherlab = ((AliMCParticle*)fMCEvent->GetTrack((Int_t)trackletLab[itracklet][0]))->GetMother();
 827                   if (motherlab>-1) {
 828                     if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==2212) fHistReconstructedProtons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 829                     else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==321) fHistReconstructedKaons->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 830                     else if (TMath::Abs(fMCEvent->GetTrack(motherlab)->PdgCode())==211) fHistReconstructedPions->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 831                     else fHistReconstructedOthers->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 832                   } else fHistReconstructedSec->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 833                 }
 834               }
 835             }
 836           } else {
 837             fHistdPhidThetaComb->Fill(trackletCoord[itracklet][0],trackletCoord[itracklet][1]);
 838             if (TMath::Abs(trackletCoord[itracklet][0])<fPhiWindowAna) {
 839               fHistBkgCorrDenPrimGen->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 840               fHistBkgCombLabels->Fill(trackletCoord[itracklet][4],esdvtx[2]);
 841             }
 842           }
 843         } // end loop tracklets
 844
 845         for (Int_t imc=0; imc<multMCCharged; imc++) {
 846           if (fTR) {
 847             if (reconstructedPrimaryPart[imc]) fHistBkgCorrNum->Fill(etagen[imc],vtxMC[2]); // only for separate corrections
 848             if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) fHistAlgEffNum->Fill(etagen[imc],vtxMC[2]);
 849             if (detectablePrimaryPart[imc]) fHistNonDetectableCorrDen->Fill(etagen[imc],vtxMC[2]);
 850           }
 851           fHistNonDetectableCorrNum->Fill(etagen[imc],vtxMC[2]);
 852         }
 853
 854         fHistSelEvts->Fill(multSPDAna,vtxMC[2]);
 855       } // end of selected events
 856
 857       fHistMCmultEtacutvsSPDmultEtacut->Fill(multMCChargedEtacut,multSPDEtacut);
 858       fHistMCmultEtacutvsSPDmultcl1->Fill(multMCChargedEtacut,multSPDcl1);
 859       fHistMCmultEtacutvsSPDmultcl2->Fill(multMCChargedEtacut,multSPDcl2);
 860       fHistAllEvts->Fill(multSPDAna,vtxMC[2]);
 861
 862       if (eventInCentralityBin) {
 863         fHistCentrEvts->Fill(multSPDAna,vtxMC[2]);
 864         if (eventTriggered) {
 865           fHistTrigCentrEvts->Fill(multSPDAna,vtxMC[2]);
 866         }
 867       }
 868
 869       for (Int_t imc=0; imc<multMCCharged; imc++) {
 870         fHistAllPrimaries->Fill(etagen[imc],vtxMC[2]);
 871         if (eventInCentralityBin) {
 872           fHistTrackCentrEvts->Fill(etagen[imc],vtxMC[2]);
 873           if (eventTriggered) fHistTrackTrigCentrEvts->Fill(etagen[imc],vtxMC[2]);
 874         }
 875       }
 876
 877       delete[] etagen;
 878       delete[] stackIndexOfPrimaryParts;
 879       delete[] reconstructedPrimaryPart;
 880       delete[] detectedPrimaryPartLay1;
 881       delete[] detectedPrimaryPartLay2;
 882       delete[] detectablePrimaryPart;
 883     }
 884   }
 885   PostData(1, fOutput);
 886 }
 887 //________________________________________________________________________
 888 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectedPrimary(Int_t nref, AliMCParticle* mcpart, Int_t Layer) {
 889   Double_t rMinLay1 = 3.4; //min rad for track ref first SPD layer
 890   Double_t rMaxLay1 = 4.4; //max rad for track ref first SPD layer
 891   Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer
 892   Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer
 893   Double_t rmin = (Layer > 0 ? rMinLay2 : rMinLay1);
 894   Double_t rmax = (Layer > 0 ? rMaxLay2 : rMaxLay1);
 895   AliTrackReference *tref=0x0;
 896
 897   for (Int_t iref=0;iref<nref;iref++) {
 898      tref = mcpart->GetTrackReference(iref);
 899    if (tref) {
 900        if (tref->R()>rmin&&tref->R()<rmax) {
 901          if (tref->DetectorId()==0) {
 902            return kTRUE;
 903          }
 904        }
 905      }
 906   }
 907   return kFALSE;
 908 }
 909 //________________________________________________________________________
 910 Bool_t AliAnalysisTaskSPDdNdEta::IsDetectablePrimary(Int_t nref, AliMCParticle* mcpart) {
 911
 912   Double_t rMinLay2 = 6.9; //min rad for track ref second SPD layer
 913   Double_t rMaxLay2 = 7.9; //max rad for track ref second SPD layer
 914
 915   AliTrackReference *tref= mcpart->GetTrackReference(nref-1);
 916   if (tref->DetectorId()!=-1) {
 917     return kTRUE;
 918   } else { //last is -1 -> particle disappeared. Where?
 919     if (tref->R()>rMaxLay2) {
 920       return kTRUE;
 921     } else if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2) { // this last tr is in lay 2
 922       for (Int_t iref=0;iref<nref;iref++) {  // look for other tr in lay 2
 923         tref = mcpart->GetTrackReference(iref);
 924         if (tref)
 925           if (tref->R()>=rMinLay2&&tref->R()<=rMaxLay2)
 926             if (tref->DetectorId()==0) return kTRUE;
 927       }
 928     } // else particle disappeared before lay2
 929   }
 930   return kFALSE;
 931 }
 932 //________________________________________________________________________
 933 void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
 934 {
 935   Printf("Terminating...");
 936 //  AliAnalysisTaskSE::Terminate();
 937
 938 }