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