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
  31 #include "AliAnalysisManager.h"
  32
  33 #include "AliMultiplicity.h"
  34 #include "AliESDInputHandler.h"
  35
  36 #include "AliMCEventHandler.h"
  37 #include "AliMCEvent.h"
  38 #include "AliStack.h"
  39 #include "AliLog.h"
  40 #include "AliTrackReference.h"
  41
  42 #include "AliGenEventHeader.h"
  43 #include "AliGenPythiaEventHeader.h"
  44 #include "AliGenDPMjetEventHeader.h"
  45 #include "AliAnalysisTaskSPDdNdEta.h"
  46
  47
  48 ClassImp(AliAnalysisTaskSPDdNdEta)
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskSPDdNdEta::AliAnalysisTaskSPDdNdEta(const char *name)
  52   : AliAnalysisTask(name, "SPDdNdEtaTask"),
  53
  54   fESD(0),
  55   fOutput(0),
  56
  57   fpythia(kTRUE),
  58   fCorr(kFALSE),
  59   fTrigger(0),
  60
  61   // Data to be corrected
  62   fHistSPDRAWMultvsZ(0),
  63   fHistSPDRAWMultvsZTriggEvts(0),
  64   fHistSPDRAWEtavsZ(0),
  65
  66   // Clusters inner layer and tracklets
  67   fHistSPDmultEtacut(0),
  68   fHistSPDmult(0),
  69   fHistSPDeta(0),
  70   fHistSPDcl1multEtacutLay1(0),
  71   fHistSPDcl1mult(0),
  72   fHistSPDcl1eta(0),
  73   fHistSPDphi(0),
  74   fHistSPDcl1phi(0),
  75   fHistSPDtheta(0),
  76   fHistSPDcl1theta(0),
  77   fHistSPDdePhi(0),
  78   fHistSPDdePhiZ(0),
  79   fHistSPDdePhi3D(0),
  80   fHistSPDphivsSPDeta(0),
  81   fHistSPDcl1phivsSPDcl1eta(0),
  82   fHistSPDdeTheta(0),
  83
  84   // SPD vertex
  85   fHistSPDvtxAnalysis(0),
  86   fHistSPDvtx3D(0),
  87   fHistSPDvtxZ(0),
  88   fHistNcontribSPDvtxvsSPDvtx(0),
  89   fHistNcontribSPDvtx3D(0),
  90   fHistNcontribSPDvtxZ(0),
  91   fHistNcontribSPDvtxall(0),
  92
  93   // SPD fired chips
  94   fHistSPDcl1multvsnFiredChipsLay1(0),
  95   fHistSPDmultvsnFiredChipsLay1(0),
  96   fHistSPDmultvsnFiredChipsLay2(0),
  97   fHistnFiredChipsLay2vsnFiredChipsLay1(0),
  98   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec(0),
  99   fHistSPDvtxRec(0),
 100
 101   // Track level correction histograms
 102   fHistBkgCorrNum(0),
 103   fHistBkgCorrDen(0),
 104
 105   fHistAlgEffNum(0),
 106
 107   fHistNonDetectableCorrNum(0),
 108   fHistNonDetectableCorrDen(0),
 109
 110   fHistTrackTrigVtxCorrNum(0),
 111   fHistTrackTrigCorrDen(0),
 112
 113   fHistTrackTrigVtxCorrNumNSD(0),
 114   fHistTrackTrigNSD(0),
 115
 116   // Event level correction histograms
 117   fHistTrigVtxCorrNum(0),
 118   fHistTrigVtxCorrDen(0),
 119
 120   fHistTrigCorrDen(0),
 121   fHistTrigVtxCorrNumNSD(0),
 122   fHistEvTrigNSD(0),
 123
 124   // MC distributions
 125   fHistMCEtavsZTriggMCvtxEvts(0),
 126   fHistMCEtavsZTriggESDvtxEvts(0),
 127   fHistMCEtavsZ(0),
 128
 129   // MC dN/dEta for each process type
 130   fHistMCEtaInel(0),
 131   fHistMCEtaNonDiffractive(0),
 132   fHistMCEtaNonSingleDiffractive(0),
 133   fHistoProcessType(0),
 134   fHistoProcessTypeTriggered(0),
 135
 136   // Additional check histos
 137   fHistContributorsvsMCVtx(0),
 138   fHistoDetectableNotr(0),
 139   fHistoDetectabletr(0),
 140   fHistoNonStoppingTracks(0),
 141   fHistoDetectedLay1(0),
 142   fHistoDetectedLay2(0),
 143   fHistoPt(0),
 144   fHistoRTRm1(0),
 145   fHistMCvtx(0),
 146
 147   // Trigger efficiency vs MC multiplicity
 148   fHistMultAllNonDiff(0),
 149
 150   fHistMultAllSingleDiff(0),
 151   fHistMultAllDoubleDiff(0),
 152   fHistMultTrVtxNonDiff(0),
 153   fHistMultTrVtxSingleDiff(0),
 154   fHistMultTrVtxDoubleDiff(0),
 155
 156   fHistMCEtaNonSingleDiffractiveLargeBin(0)
 157
 158 {
 159
 160   // Constructor
 161
 162   // Define input and output slots here
 163
 164   DefineInput(0, TChain::Class());
 165   DefineOutput(0, TList::Class());
 166
 167 }
 168 //________________________________________________________________________
 169 AliAnalysisTaskSPDdNdEta::~AliAnalysisTaskSPDdNdEta()
 170 {
 171
 172   // Destructor
 173
 174   // histograms are in the output list and deleted when the output
 175   // list is deleted by the TSelector dtor
 176
 177   if (fOutput) {
 178     delete fOutput;
 179     fOutput = 0;
 180   }
 181
 182 }
 183 //________________________________________________________________________
 184 void AliAnalysisTaskSPDdNdEta::ConnectInputData(Option_t *)
 185 {
 186
 187   // Connect ESD
 188   // Called once
 189
 190   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 191
 192   if (!tree) {
 193     Printf("ERROR: Could not read chain from input slot 0");
 194   } else {
 195     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 196     if (!esdH) {
 197       Printf("ERROR: Could not get ESDInputHandler");
 198     } else fESD = esdH->GetEvent();
 199   }
 200   // Disable info messages of AliMCEvent (per event)
 201   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 202
 203 }
 204
 205 //________________________________________________________________________
 206 void AliAnalysisTaskSPDdNdEta::CreateOutputObjects()
 207 {
 208
 209   // Create histograms
 210   fOutput = new TList;
 211   fOutput->SetOwner();
 212
 213   const Int_t nBinMultCorr = 28;
 214   Double_t binLimMultCorr[nBinMultCorr+1];
 215   for (Int_t i = 0; i<nBinMultCorr+1; ++i) {
 216     if (i<21) binLimMultCorr[i] = (Double_t) i;
 217     else if (i<27) binLimMultCorr[i] = (Double_t) 20 +5*(i-20);
 218     else if (i==27) binLimMultCorr[i] = 100;
 219     else if (i==28) binLimMultCorr[i] = 200;
 220   }
 221
 222   // Event level correction
 223   fHistSPDRAWMultvsZ= new TH2F("fHistSPDRAWMultvsZ", "",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 224   fHistSPDRAWMultvsZ->Sumw2();
 225   fOutput->Add(fHistSPDRAWMultvsZ);
 226
 227   fHistSPDRAWMultvsZTriggEvts = new TH2F("fHistSPDRAWMultvsZTriggEvts", "",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 228   fHistSPDRAWMultvsZTriggEvts->Sumw2();
 229   fOutput->Add(fHistSPDRAWMultvsZTriggEvts);
 230
 231   fHistSPDRAWEtavsZ = new TH2F("fHistSPDRAWEtavsZ", "Tracklet pseudorapidity distribution", 120, -3.,3.,40,-20.,20.);
 232   fHistSPDRAWEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
 233   fHistSPDRAWEtavsZ->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 234   fHistSPDRAWEtavsZ->Sumw2();
 235   fOutput->Add(fHistSPDRAWEtavsZ);
 236
 237   fHistSPDmultEtacut = new TH1F("fHistSPDmultEtacut", "Tracklet multiplicity distribution",201,-0.5,200.5);
 238   fHistSPDmultEtacut->GetXaxis()->SetTitle("Reconstructed multiplicity (|#eta|<1.5)");
 239   fHistSPDmultEtacut->GetYaxis()->SetTitle("Entries");
 240   fOutput->Add(fHistSPDmultEtacut);
 241
 242   fHistSPDmult = new TH1F("fHistSPDmult", "Tracklet multiplicity distribution", 201,-0.5,200.5);
 243   fHistSPDmult->GetXaxis()->SetTitle("Reconstructed tracklet multiplicity");
 244   fHistSPDmult->GetYaxis()->SetTitle("Entries");
 245   fOutput->Add(fHistSPDmult);
 246
 247   fHistSPDeta = new TH1F("fHistSPDeta", "Tracklet pseudorapidity distribution", 120, -3.,3.);
 248   fHistSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 249   fHistSPDeta->GetYaxis()->SetTitle("Entries");
 250   fHistSPDeta->SetLineColor(kGreen);
 251   fHistSPDeta->SetLineWidth(3);
 252   fHistSPDeta->Sumw2();
 253   fOutput->Add(fHistSPDeta);
 254
 255   fHistSPDcl1multEtacutLay1 = new TH1F("fHistSPDcl1multEtacutLay1", "Cluster multiplicity (inner layer)",201,-0.5,200.5);
 256   fHistSPDcl1multEtacutLay1->GetXaxis()->SetTitle("Cluster multiplicity lay1 (|#eta|<2.)");
 257   fHistSPDcl1multEtacutLay1->GetYaxis()->SetTitle("Entries");
 258   fOutput->Add(fHistSPDcl1multEtacutLay1);
 259
 260   fHistSPDcl1mult = new TH1F("fHistSPDcl1mult", "Cluster multiplicity (inner layer)",201,-0.5,200.5);
 261   fHistSPDcl1mult->GetXaxis()->SetTitle("Cluster multiplicity lay1");
 262   fHistSPDcl1mult->GetYaxis()->SetTitle("Entries");
 263   fOutput->Add(fHistSPDcl1mult);
 264
 265   fHistSPDcl1eta  = new TH1F("fHistSPDcl1eta", "Cluster pseudorapidity (inner layer)", 120, -3.,3.);
 266   fHistSPDcl1eta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 267   fHistSPDcl1eta->GetYaxis()->SetTitle("Entries");
 268   fHistSPDcl1eta->Sumw2();
 269   fOutput->Add(fHistSPDcl1eta);
 270
 271   fHistSPDphi = new TH1F("fHistSPDphi", "Tracklet #phi  distribution", 360, 0.,2*TMath::Pi());
 272   fHistSPDphi->GetXaxis()->SetTitle("#varphi [rad]");
 273   fHistSPDphi->GetYaxis()->SetTitle("Entries");
 274   fOutput->Add(fHistSPDphi);
 275
 276   fHistSPDcl1phi= new TH1F("fHistSPDcl1phi", "Cluster #phi (inner layer) ", 360, 0.,2*TMath::Pi());
 277   fHistSPDcl1phi->GetXaxis()->SetTitle("#varphi [rad]");
 278   fHistSPDcl1phi->GetYaxis()->SetTitle("Entries");
 279   fOutput->Add(fHistSPDcl1phi);
 280
 281   fHistSPDtheta = new TH1F("fHistSPDtheta", "Tracklet #theta  distribution", 180, 0.,TMath::Pi());
 282   fHistSPDtheta->GetXaxis()->SetTitle("#theta [rad]");
 283   fHistSPDtheta->GetYaxis()->SetTitle("Entries");
 284   fOutput->Add(fHistSPDtheta);
 285
 286   fHistSPDcl1theta = new TH1F("fHistSPDcl1theta", "Cluster #theta (inner layer)", 180, 0.,TMath::Pi());
 287   fHistSPDcl1theta->GetXaxis()->SetTitle("#theta [rad]");
 288   fHistSPDcl1theta->GetYaxis()->SetTitle("Entries");
 289   fOutput->Add(fHistSPDcl1theta);
 290
 291   fHistSPDdePhi= new TH1F("fHistSPDdePhi", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
 292   fHistSPDdePhi->GetXaxis()->SetTitle("#Delta#varphi [rad]");
 293   fHistSPDdePhi->GetYaxis()->SetTitle("Entries");
 294   fOutput->Add(fHistSPDdePhi);
 295
 296   fHistSPDdePhiZ= new TH1F("fHistSPDdePhiZ", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
 297   fOutput->Add(fHistSPDdePhiZ);
 298
 299   fHistSPDdePhi3D= new TH1F("fHistSPDdePhi3D", "Tracklet #Delta#varphi  distribution",400,-0.1,.1);
 300   fOutput->Add(fHistSPDdePhi3D);
 301
 302   fHistSPDphivsSPDeta= new TH2F("fHistSPDphivsSPDeta", "Tracklets - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());
 303   fHistSPDphivsSPDeta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 304   fHistSPDphivsSPDeta->GetYaxis()->SetTitle("#varphi [rad]");
 305   fOutput->Add(fHistSPDphivsSPDeta);
 306
 307   fHistSPDcl1phivsSPDcl1eta= new TH2F("fHistSPDcl1phivsSPDcl1eta", "Clusters layer1 - #varphi vs #eta",120,-3.,3,360,0.,2*TMath::Pi());
 308   fHistSPDcl1phivsSPDcl1eta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 309   fHistSPDcl1phivsSPDcl1eta->GetYaxis()->SetTitle("#varphi [rad]");
 310   fOutput->Add(fHistSPDcl1phivsSPDcl1eta);
 311
 312   fHistSPDdeTheta= new TH1F("fHistSPDdeTheta", "Tracklet #Delta#theta distribution",100,-0.05,.05);
 313   fHistSPDdeTheta->GetXaxis()->SetTitle("#Delta#theta [rad]");
 314   fHistSPDdeTheta->GetYaxis()->SetTitle("Entries");
 315   fOutput->Add(fHistSPDdeTheta);
 316
 317
 318   fHistSPDvtxAnalysis = new TH1F("fHistSPDvtxAnalysis", "SPD vertex  distribution - all events",20,-20.,20.);
 319   fHistSPDvtxAnalysis->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 320   fHistSPDvtxAnalysis->GetYaxis()->SetTitle("Entries");
 321   fOutput->Add(fHistSPDvtxAnalysis);
 322
 323   fHistSPDvtx3D = new TH3F("fHistSPDvtx3D", "SPD vertex distribution",100,-1.,1.,100,-1.,1.,500,-50.,50.);
 324   fOutput->Add(fHistSPDvtx3D);
 325
 326   fHistSPDvtxZ = new TH1F("fHistSPDvtxZ", "SPD vertex distribution",500,-50.,50.);
 327   fOutput->Add(fHistSPDvtxZ);
 328
 329   fHistNcontribSPDvtxvsSPDvtx= new TH2F("fHistNcontribSPDvtxvsSPDvtx", " ",100,-50.,50.,10002,-2.,10000.);
 330   fHistNcontribSPDvtxvsSPDvtx->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 331   fHistNcontribSPDvtxvsSPDvtx->GetYaxis()->SetTitle("# contributors");
 332   fOutput->Add(fHistNcontribSPDvtxvsSPDvtx);
 333
 334   fHistNcontribSPDvtx3D= new TH1F("fHistNcontribSPDvtx3D", "SPD vtx 3D",10002,-2.,10000.);
 335   fHistNcontribSPDvtx3D->GetXaxis()->SetTitle("# contributors");
 336   fHistNcontribSPDvtx3D->GetYaxis()->SetTitle("Entries");
 337   fOutput->Add(fHistNcontribSPDvtx3D);
 338
 339   fHistNcontribSPDvtxZ= new TH1F("fHistNcontribSPDvtxZ", "SPD vtx Z",10002,-2.,10000.);
 340   fHistNcontribSPDvtxZ->GetXaxis()->SetTitle("# contributors");
 341   fHistNcontribSPDvtxZ->GetYaxis()->SetTitle("Entries");
 342   fOutput->Add(fHistNcontribSPDvtxZ);
 343
 344   fHistNcontribSPDvtxall= new TH1F("fHistNcontribSPDvtxall", "SPD vtx - all events",10002,-2.,10000.);
 345   fHistNcontribSPDvtxall->GetXaxis()->SetTitle("# contributors");
 346   fHistNcontribSPDvtxall->GetYaxis()->SetTitle("Entries");
 347   fOutput->Add(fHistNcontribSPDvtxall);
 348
 349   fHistSPDcl1multvsnFiredChipsLay1 = new TH2F("fHistSPDcl1multvsnFiredChipsLay1", "",401,0.,401.,201,-0.5,200.5);
 350   fHistSPDcl1multvsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
 351   fHistSPDcl1multvsnFiredChipsLay1->GetYaxis()->SetTitle("Cluster lay1 multiplicity");
 352   fOutput->Add(fHistSPDcl1multvsnFiredChipsLay1);
 353
 354   fHistSPDmultvsnFiredChipsLay1 = new TH2F("fHistSPDmultvsnFiredChipsLay1","",401,0.,401.,201,-0.5,200.5);
 355   fHistSPDmultvsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
 356   fHistSPDmultvsnFiredChipsLay1->GetYaxis()->SetTitle("Tracklet multiplicity");
 357   fOutput->Add(fHistSPDmultvsnFiredChipsLay1);
 358
 359   fHistSPDmultvsnFiredChipsLay2 = new TH2F("fHistSPDmultvsnFiredChipsLay2","",801,0.,801.,201,-0.5,200.5);
 360   fHistSPDmultvsnFiredChipsLay2->GetXaxis()->SetTitle("# fired chips lay2");
 361   fHistSPDmultvsnFiredChipsLay2->GetYaxis()->SetTitle("Tracklet multiplicity");
 362   fOutput->Add(fHistSPDmultvsnFiredChipsLay2);
 363
 364   fHistnFiredChipsLay2vsnFiredChipsLay1 = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1","",401,0.,401.,801,0.,801.);
 365   fHistnFiredChipsLay2vsnFiredChipsLay1->GetXaxis()->SetTitle("# fired chips lay1");
 366   fHistnFiredChipsLay2vsnFiredChipsLay1->GetYaxis()->SetTitle("# fired chip lay2");
 367   fOutput->Add(fHistnFiredChipsLay2vsnFiredChipsLay1);
 368
 369   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec = new TH2F("fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec","",401,0.,401.,801,0.,801.);
 370   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec->GetXaxis()->SetTitle("# fired chips lay1");
 371   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec->GetYaxis()->SetTitle("# fired chip lay2");
 372   fOutput->Add(fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec);
 373
 374   fHistSPDvtxRec = new TH1F("fHistSPDvtxRec", "SPD vertex  distribution",80,-20.,20.);
 375   fHistSPDvtxRec->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 376   fHistSPDvtxRec->GetYaxis()->SetTitle("Entries");
 377   fOutput->Add(fHistSPDvtxRec);
 378
 379
 380   if (fCorr) {
 381     fHistBkgCorrNum = new TH2F("fHistBkgCorrNum","",60,-3.00,3.00,20,-20.,20.);
 382     fOutput->Add(fHistBkgCorrNum);
 383
 384     fHistBkgCorrDen = new TH2F("fHistBkgCorrDen","",60,-3.00,3.00,20,-20.,20.);
 385     fOutput->Add(fHistBkgCorrDen);
 386
 387     fHistAlgEffNum = new TH2F("fHistAlgEffNum","",120,-3.00,3.00,40,-20.,20.);
 388     fOutput->Add(fHistAlgEffNum);
 389
 390     fHistNonDetectableCorrNum = new TH2F("fHistNonDetectableCorrNum","",60,-3.00,3.00,20,-20.,20.);
 391     fOutput->Add(fHistNonDetectableCorrNum);
 392
 393     fHistNonDetectableCorrDen = new TH2F("fHistNonDetectableCorrDen","",120,-3.00,3.00,40,-20.,20.);
 394     fOutput->Add(fHistNonDetectableCorrDen);
 395
 396     fHistTrackTrigVtxCorrNum = new TH2F("fHistTrackTrigVtxCorrNum","",60,-3.00,3.00,20,-20.,20.);
 397     fOutput->Add(fHistTrackTrigVtxCorrNum);
 398
 399     fHistTrackTrigCorrDen = new TH2F("fHistTrackTrigCorrDen","",60,-3.00,3.00,20,-20.,20.);
 400     fOutput->Add(fHistTrackTrigCorrDen);
 401
 402     fHistTrackTrigVtxCorrNumNSD = new TH2F("fHistTrackTrigVtxCorrNumNSD","",60,-3.00,3.00,20,-20.,20.);
 403     fOutput->Add(fHistTrackTrigVtxCorrNumNSD);
 404
 405     fHistTrackTrigNSD = new TH2F("fHistTrackTrigNSD","",60,-3.00,3.00,20,-20.,20.);
 406     fOutput->Add(fHistTrackTrigNSD);
 407
 408
 409     // Event level correction histograms
 410     fHistTrigVtxCorrNum = new TH2F("fHistTrigVtxCorrNum","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 411     fOutput->Add(fHistTrigVtxCorrNum);
 412
 413     fHistTrigVtxCorrDen = new TH2F("fHistTrigVtxCorrDen","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 414     fOutput->Add(fHistTrigVtxCorrDen);
 415
 416     fHistTrigCorrDen = new TH2F("fHistTrigCorrDen","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 417     fOutput->Add(fHistTrigCorrDen);
 418
 419     fHistTrigVtxCorrNumNSD = new TH2F("fHistTrigVtxCorrNumNSD","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 420     fOutput->Add(fHistTrigVtxCorrNumNSD);
 421
 422     fHistEvTrigNSD = new TH2F("fHistEvTrigNSD","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 423     fOutput->Add(fHistEvTrigNSD);
 424
 425     // MC distributions
 426     fHistMCEtavsZTriggMCvtxEvts = new TH2F("fHistMCEtavsZTriggMCvtxEvts","Generated pseudorapidity distribution",60,-3.00,3.00,20,-20.,20.);
 427     fHistMCEtavsZTriggMCvtxEvts->GetXaxis()->SetTitle("Pseudorapidity #eta");
 428     fHistMCEtavsZTriggMCvtxEvts->GetYaxis()->SetTitle("z_{MCvtx} [cm]");
 429     fHistMCEtavsZTriggMCvtxEvts->Sumw2();
 430     fOutput->Add(fHistMCEtavsZTriggMCvtxEvts);
 431
 432     fHistMCEtavsZTriggESDvtxEvts = new TH2F("fHistMCEtavsZTriggESDvtxEvts","Generated pseudorapidity distribution",60,-3.00,3.00,20,-20.,20.);
 433     fHistMCEtavsZTriggESDvtxEvts->GetXaxis()->SetTitle("Pseudorapidity #eta");
 434     fHistMCEtavsZTriggESDvtxEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 435     fHistMCEtavsZTriggESDvtxEvts->Sumw2();
 436     fOutput->Add(fHistMCEtavsZTriggESDvtxEvts);
 437
 438     fHistMCEtavsZ = new TH2F("fHistMCEtavsZ","Generated pseudorapidity distribution",60,-3.00,3.00,20,-20.,20.);
 439     fHistMCEtavsZ->GetXaxis()->SetTitle("Pseudorapidity #eta");
 440     fHistMCEtavsZ->GetYaxis()->SetTitle("z_{MCvtx} [cm]");
 441     fHistMCEtavsZ->Sumw2();
 442     fOutput->Add(fHistMCEtavsZ);
 443
 444     fHistMCEtaInel = new TH1F("fHistMCEtaInel","",7,-3.5,3.5);
 445     fHistMCEtaInel->GetXaxis()->SetTitle("Pseudorapidity #eta");
 446     fOutput->Add(fHistMCEtaInel);
 447
 448     fHistMCEtaNonDiffractive = new TH1F("fHistMCEtaNonDiffractive","",60,-3.,3.);
 449     fHistMCEtaNonDiffractive->GetXaxis()->SetTitle("Pseudorapidity #eta");
 450     fOutput->Add(fHistMCEtaNonDiffractive);
 451
 452     fHistMCEtaNonSingleDiffractive = new TH1F("fHistMCEtaNonSingleDiffractive","",60,-3.,3.);
 453     fHistMCEtaNonSingleDiffractive->GetXaxis()->SetTitle("Pseudorapidity #eta");
 454     fOutput->Add(fHistMCEtaNonSingleDiffractive);
 455
 456
 457
 458     fHistoProcessType = new TH1F("fHistoProcessType","",5,0.,5.);
 459     fOutput->Add(fHistoProcessType);
 460
 461     fHistoProcessTypeTriggered = new TH1F("fHistoProcessTypeTriggered","",5,0.,5);
 462     fOutput->Add(fHistoProcessTypeTriggered);
 463
 464     fHistContributorsvsMCVtx = new TH2F("fHistContributorsvsMCVtx","",200,-20.,20.,202,-2.,200.);
 465     fOutput->Add(fHistContributorsvsMCVtx);
 466
 467     fHistoDetectableNotr = new TH3F("fHistoDetectableNotr","",60,-3.00,3.00,20,-20.,20.,100,0.,10.);
 468     fOutput->Add(fHistoDetectableNotr);
 469
 470     fHistoDetectabletr = new TH2F("fHistoDetectabletr","",60,-3.00,3.00,20,-20.,20.);
 471     fOutput->Add(fHistoDetectabletr);
 472
 473     fHistoNonStoppingTracks = new TH2F("fHistoNonStoppingTracks","",60,-3.00,3.00,20,-20.,20.);
 474     fOutput->Add(fHistoNonStoppingTracks);
 475
 476     fHistoDetectedLay1 = new TH2F("fHistoDetectedLay1","",60,-3.00,3.00,20,-20.,20.);
 477     fOutput->Add(fHistoDetectedLay1);
 478
 479     fHistoDetectedLay2 = new TH2F("fHistoDetectedLay2","",60,-3.00,3.00,20,-20.,20.);
 480     fOutput->Add(fHistoDetectedLay2);
 481
 482     fHistoPt = new TH1F("fHistoPt","",100,.0,10.);
 483     fOutput->Add(fHistoPt);
 484
 485     fHistoRTRm1 = new TH1F("fHistoRTRm1","",10000,0.,5000);
 486     fOutput->Add(fHistoRTRm1);
 487
 488     fHistMCvtx = new TH3F("fHistMCvtx", "MC vertex distribution",100,-.5,.5,100,-.5,.5,500,-50.,50.);
 489     fOutput->Add(fHistMCvtx);
 490
 491     fHistMultAllNonDiff  = new TH1F("fHistMultAllNonDiff","",20,-0.5,19.5);
 492     fOutput->Add(fHistMultAllNonDiff);
 493     fHistMultAllSingleDiff = new TH1F("fHistMultAllSingleDiff","",20,-0.5,19.5);
 494     fHistMultAllSingleDiff->Sumw2();
 495     fOutput->Add(fHistMultAllSingleDiff);
 496     fHistMultAllDoubleDiff = new TH1F("fHistMultAllDoubleDiff","",20,-0.5,19.5);
 497     fHistMultAllDoubleDiff->Sumw2();
 498     fOutput->Add(fHistMultAllDoubleDiff);
 499     fHistMultTrVtxNonDiff = new TH1F("fHistMultTrVtxNonDiff","",20,-0.5,19.5);
 500     fHistMultTrVtxNonDiff->Sumw2();
 501     fOutput->Add(fHistMultTrVtxNonDiff);
 502     fHistMultTrVtxSingleDiff = new TH1F("fHistMultTrVtxSingleDiff","",20,-0.5,19.5);
 503     fHistMultTrVtxSingleDiff->Sumw2();
 504     fOutput->Add(fHistMultTrVtxSingleDiff);
 505     fHistMultTrVtxDoubleDiff = new TH1F("fHistMultTrVtxDoubleDiff","",20,-0.5,19.5);
 506     fHistMultTrVtxDoubleDiff->Sumw2();
 507     fOutput->Add(fHistMultTrVtxDoubleDiff);
 508
 509     fHistMCEtaNonSingleDiffractiveLargeBin = new TH1F("fHistMCEtaNonSingleDiffractiveLargeBin","",7,-3.5,3.5);
 510     fHistMCEtaNonSingleDiffractiveLargeBin->GetXaxis()->SetTitle("Pseudorapidity #eta");
 511     fOutput->Add(fHistMCEtaNonSingleDiffractiveLargeBin);
 512
 513   }
 514 }
 515
 516 //________________________________________________________________________
 517 void AliAnalysisTaskSPDdNdEta::Exec(Option_t *)
 518 {
 519   // Main loop
 520   // Called for each event
 521
 522
 523   if (!fESD) {
 524     Printf("ERROR: fESD not available");
 525     return;
 526   }
 527
 528   // ESD vertex
 529   Bool_t eventWithVertex = kFALSE;
 530   const AliESDVertex* vtxESD = fESD->GetVertex();
 531   Double_t esdvtx[3];
 532   vtxESD->GetXYZ(esdvtx);
 533   Int_t nContrib = vtxESD->GetNContributors();
 534   //... check resolution
 535   Double_t zRes = vtxESD->GetZRes();
 536   const AliMultiplicity* multESD = fESD->GetMultiplicity();
 537
 538   // Loading tracklets...
 539   Int_t multSPD = 0;
 540   multSPD = multESD->GetNumberOfTracklets();
 541
 542   AliMultiplicity * mult1 = (AliMultiplicity*)multESD;
 543   Short_t nFiredChipsLay1 = mult1->GetNumberOfFiredChips(0);
 544   Short_t nFiredChipsLay2 = mult1->GetNumberOfFiredChips(1);
 545   Int_t multSPDEtacut = 0;
 546   Int_t multSPDcl1 = 0;
 547   Int_t nSingleCl1 = 0;
 548   Int_t multSPDcl1EtacutLay1 = 0;
 549   nSingleCl1 = multESD->GetNumberOfSingleClusters();
 550   multSPDcl1 = nSingleCl1 + multSPD;
 551   Float_t* recEtaSPDcl1 = new Float_t[multSPD+nSingleCl1];
 552   if (esdvtx[2]!=0.) fHistSPDvtxRec->Fill(esdvtx[2]);
 553
 554   if (esdvtx[2]!=0.&&zRes<0.1) eventWithVertex = kTRUE;
 555   else {
 556     esdvtx[2] = 0.;
 557     multSPD = 0;
 558   }
 559
 560 //  Printf("There are %d tracklets in this event", multSPD);
 561
 562   // Trigger
 563   Bool_t eventTriggered = kFALSE;
 564   ULong64_t triggerMask;
 565   ULong64_t spdFO   = (1 << 14);
 566   ULong64_t v0left  = (1 << 11);
 567   ULong64_t v0right = (1 << 12);
 568
 569   triggerMask=fESD->GetTriggerMask();
 570
 571   // No trigger
 572   if (fTrigger==0) eventTriggered = kTRUE;
 573   // MB1: GFO || V0OR
 574   if (fTrigger==1) eventTriggered = triggerMask&spdFO || ((triggerMask&v0left) || (triggerMask&v0right));
 575   // MB2: GFO && V0OR
 576   if (fTrigger==2) eventTriggered = triggerMask&spdFO && ((triggerMask&v0left) || (triggerMask&v0right));
 577
 578   PostData(0, fOutput);
 579
 580
 581   // Selected events: triggered with vertex
 582   if (eventTriggered&&eventWithVertex) {
 583
 584     if (multSPD!=0) fHistSPDvtxAnalysis->Fill(esdvtx[2]);
 585     if (strcmp(vtxESD->GetTitle(),"vertexer: 3D") == 0) fHistSPDvtx3D->Fill(esdvtx[0],esdvtx[1],esdvtx[2]);
 586     else fHistSPDvtxZ->Fill(esdvtx[2]);
 587
 588     for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 589       Float_t thetaTr= multESD->GetTheta(itracklet);
 590       Float_t phiTr= multESD->GetPhi(itracklet);
 591       Float_t dePhiTr= multESD->GetDeltaPhi(itracklet);
 592 //      Double_t deThetaTr= multESD->GetDeltaTheta(itracklet);
 593       Float_t recEtaSPD =multESD->GetEta(itracklet);
 594       recEtaSPDcl1[itracklet] = recEtaSPD;
 595
 596       fHistSPDeta->Fill(recEtaSPD);
 597       fHistSPDRAWEtavsZ->Fill(recEtaSPD,esdvtx[2]);
 598       fHistSPDcl1eta->Fill(recEtaSPD);
 599       fHistSPDphi->Fill(phiTr);
 600       fHistSPDcl1phi->Fill(phiTr);
 601       fHistSPDtheta->Fill(thetaTr);
 602       fHistSPDcl1theta->Fill(thetaTr);
 603       fHistSPDdePhi->Fill(dePhiTr);
 604 //      fHistSPDdeTheta->Fill(deThetaTr);
 605
 606       if (strcmp(vtxESD->GetTitle(),"vertexer: Z") == 0) fHistSPDdePhiZ->Fill(dePhiTr);
 607       if (strcmp(vtxESD->GetTitle(),"vertexer: 3D") == 0) fHistSPDdePhi3D->Fill(dePhiTr);
 608       fHistSPDphivsSPDeta->Fill(recEtaSPD,phiTr);
 609       fHistSPDcl1phivsSPDcl1eta->Fill(recEtaSPD,phiTr);
 610
 611       // Calculate multiplicity in etacut
 612       if (TMath::Abs(recEtaSPD)<1.5) multSPDEtacut++;
 613       if (TMath::Abs(recEtaSPD)<2.)  multSPDcl1EtacutLay1++;
 614     }
 615
 616     for (Int_t iCl1=0; iCl1<nSingleCl1; ++iCl1) {
 617       Float_t thetaSingleCl1 = multESD->GetThetaSingle(iCl1);
 618       // Calculate eta
 619       Float_t etaSingleCl1 = -TMath::Log(TMath::Tan(thetaSingleCl1/2.));
 620       Float_t phiSingleCl1 = multESD->GetPhiSingle(iCl1);
 621       recEtaSPDcl1[iCl1+multSPD] = etaSingleCl1;
 622
 623       fHistSPDcl1eta->Fill(etaSingleCl1);
 624       fHistSPDcl1phi->Fill(phiSingleCl1);
 625       fHistSPDcl1theta->Fill(thetaSingleCl1);
 626       fHistSPDcl1phivsSPDcl1eta->Fill(etaSingleCl1,phiSingleCl1);
 627       if (TMath::Abs(etaSingleCl1)<2.) multSPDcl1EtacutLay1++;
 628     }
 629
 630     fHistSPDmultEtacut->Fill(multSPDEtacut);
 631     fHistSPDmult->Fill(multSPD);
 632     fHistSPDcl1multEtacutLay1->Fill(multSPDcl1EtacutLay1);
 633     if (strcmp(vtxESD->GetTitle(),"vertexer: 3D") == 0) {
 634       fHistNcontribSPDvtx3D->Fill(nContrib);
 635     }
 636     if (strcmp(vtxESD->GetTitle(),"vertexer: Z") == 0)  {
 637      fHistNcontribSPDvtxZ->Fill(nContrib);
 638     }
 639     fHistNcontribSPDvtxvsSPDvtx->Fill(esdvtx[2],nContrib);
 640     fHistSPDRAWMultvsZ->Fill(multSPD,esdvtx[2]);
 641     fHistSPDmultvsnFiredChipsLay1->Fill(nFiredChipsLay1,multSPD);
 642     fHistSPDmultvsnFiredChipsLay2->Fill(nFiredChipsLay2,multSPD);
 643
 644   } // End selected events
 645
 646   if (eventTriggered) {
 647     fHistSPDRAWMultvsZTriggEvts->Fill(multSPD,esdvtx[2]);
 648   }
 649
 650   fHistSPDcl1mult->Fill(multSPDcl1);
 651   fHistNcontribSPDvtxall->Fill(nContrib);
 652
 653   fHistSPDcl1multvsnFiredChipsLay1->Fill(nFiredChipsLay1,multSPDcl1);
 654
 655   fHistnFiredChipsLay2vsnFiredChipsLay1->Fill(nFiredChipsLay1,nFiredChipsLay2);
 656
 657   if (esdvtx[2]==0.) {
 658     fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec->Fill(nFiredChipsLay1,nFiredChipsLay2);
 659   }
 660
 661   delete[] recEtaSPDcl1;
 662
 663   if (fCorr) {
 664     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 665     if (!eventHandler) {
 666       Printf("ERROR: Could not retrieve MC event handler");
 667       return;
 668     }
 669
 670     AliMCEvent* mcEvent = eventHandler->MCEvent();
 671     if (!mcEvent) {
 672        Printf("ERROR: Could not retrieve MC event");
 673        return;
 674     }
 675
 676     AliStack* stack = mcEvent->Stack();
 677     if (!stack) {
 678       AliDebug(AliLog::kError, "Stack not available");
 679       return;
 680     }
 681
 682     AliHeader* header = mcEvent->Header();
 683     if (!header) {
 684       AliDebug(AliLog::kError, "Header not available");
 685       return;
 686     }
 687     AliGenEventHeader* genHeader = header->GenEventHeader();
 688     // MC vertex
 689     TArrayF vtxMC(3);
 690     genHeader->PrimaryVertex(vtxMC);
 691
 692     fHistMCvtx->Fill(vtxMC[0],vtxMC[1],vtxMC[2]);
 693
 694     //Adding process type selection
 695     Int_t processType;
 696     Bool_t nsdEv = kFALSE;
 697     Bool_t ndEv = kFALSE;
 698     Bool_t sdEv = kFALSE;
 699     Bool_t ddEv = kFALSE;
 700     Bool_t inelEv = kFALSE;
 701
 702     AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 703
 704     AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 705
 706     if (fpythia&&pythiaGenHeader) {
 707       processType = pythiaGenHeader->ProcessType();
 708       if (processType!=91&&processType!=92&&processType!=93&&processType!=94)  ndEv = kTRUE; // non difffractive
 709       if (processType!=91&&processType!=92&&processType!=93)  nsdEv = kTRUE;                 // non single diffractive
 710       if (processType==92||processType==93)  sdEv = kTRUE;                                   //single diffractive
 711       if (processType==94)  ddEv = kTRUE;                                                    //double diffractive
 712       if (processType!=91)  inelEv = kTRUE;                                                  //inelastic
 713
 714     } else if (!fpythia&&dpmHeader) {
 715       processType = dpmHeader->ProcessType();
 716       if (processType==1)  ndEv = kTRUE;                                         // non diffractive
 717       if (processType!=2&&processType!=5&&processType!=6)  nsdEv = kTRUE;        // non single diffractive
 718       if (processType==5||processType==6)  sdEv = kTRUE;                         // single diffractive
 719       if (processType==4||processType==7)  ddEv = kTRUE;                         // double diffractive
 720       if (processType!=2)  inelEv = kTRUE;                                       // inelastic
 721     } else if (!pythiaGenHeader&&!dpmHeader) {
 722       printf("Unknown header type: neither DPMjet nor Pythia. \n");
 723       return ;
 724     }
 725
 726     if (ndEv) {
 727       fHistoProcessType->Fill(0);
 728       if (eventTriggered) fHistoProcessTypeTriggered->Fill(0);
 729     }
 730     if (nsdEv) {
 731       fHistoProcessType->Fill(1);
 732       if (eventTriggered) fHistoProcessTypeTriggered->Fill(1);
 733     } else if (sdEv) {
 734       fHistoProcessType->Fill(3);
 735       if (eventTriggered) fHistoProcessTypeTriggered->Fill(3);
 736     }
 737     if (ddEv) {
 738       fHistoProcessType->Fill(4);
 739       if (eventTriggered) fHistoProcessTypeTriggered->Fill(4);
 740     }
 741     if (inelEv) {
 742       fHistoProcessType->Fill(2);                          // inel is always true
 743       if (eventTriggered) fHistoProcessTypeTriggered->Fill(2);
 744     }
 745
 746     // Tracks from MC
 747     Int_t  multMCCharged = 0;
 748     Int_t  multMCChargedEtacut = 0;
 749     Int_t  nMCPart = stack->GetNprimary();
 750     Float_t* etagen = new Float_t[nMCPart];
 751     Float_t* ptgen = new Float_t[nMCPart];
 752     Int_t* stackIndexOfPrimaryParts = new Int_t[nMCPart];
 753     Int_t* reconstructedPrimaryPart = new Int_t[nMCPart];
 754     Int_t* detectedPrimaryPartLay1 = new Int_t[nMCPart];
 755     Int_t* detectedPrimaryPartLay2 = new Int_t[nMCPart];
 756     Int_t* detectablePrimaryPart = new Int_t[nMCPart];
 757
 758     // Loading track references...
 759     Float_t rminL1 = 3.4;
 760     Float_t rmaxL1 = 4.4;
 761     Float_t rminL2 = 6.9;
 762     Float_t rmaxL2 = 7.9;
 763
 764     TTree* tRefTree = eventHandler->TreeTR();
 765
 766     AliTrackReference *tref=0x0;
 767     mcEvent->ConnectTreeTR(tRefTree);
 768
 769     // Loop over MC particles
 770     for (Int_t imc=0; imc<nMCPart; imc++) {
 771       TParticle* part = stack->Particle(imc);
 772       Bool_t isPrimary = stack->IsPhysicalPrimary(imc);
 773       if (!isPrimary)                        continue;
 774       TParticlePDG* pdgPart = part->GetPDG();
 775       if (TMath::Abs(pdgPart->Charge())!=3)  continue;
 776       Float_t theta = part->Theta();
 777       if (theta==0 || theta==TMath::Pi())    continue;
 778       Float_t eta = part->Eta();
 779       Float_t pt = part->Pt();
 780       etagen[multMCCharged] = eta;
 781       ptgen[multMCCharged] = pt;
 782       stackIndexOfPrimaryParts[multMCCharged] = imc;
 783
 784       reconstructedPrimaryPart[multMCCharged]=kFALSE;
 785       detectedPrimaryPartLay1[multMCCharged]=kFALSE;
 786       detectedPrimaryPartLay2[multMCCharged]=kFALSE;
 787       detectablePrimaryPart[multMCCharged]=kFALSE;
 788
 789       fHistoPt->Fill(ptgen[multMCCharged]);
 790       if (ndEv)                                                    // non difffractive
 791         fHistMCEtaNonDiffractive->Fill(etagen[multMCCharged]);
 792       if (nsdEv) {                                                 // non single diffractive
 793         fHistMCEtaNonSingleDiffractive->Fill(etagen[multMCCharged]);
 794         fHistMCEtaNonSingleDiffractiveLargeBin->Fill(etagen[multMCCharged]);
 795       }
 796       // inel is always true
 797       fHistMCEtaInel->Fill(etagen[multMCCharged]);
 798
 799
 800       AliMCParticle* mcpart = (AliMCParticle*) mcEvent->GetTrack(imc);
 801       Int_t nref = mcpart->GetNumberOfTrackReferences();
 802
 803       // Detectable primaries
 804       if (nref==0) {
 805         detectablePrimaryPart[multMCCharged]=kTRUE;
 806         fHistoDetectableNotr->Fill(etagen[multMCCharged],vtxMC[2],ptgen[multMCCharged]);
 807       } else if (nref>0) {
 808         tref = mcpart->GetTrackReference(nref-1);
 809         if (tref->DetectorId()!=-1) {
 810           detectablePrimaryPart[multMCCharged]=kTRUE;
 811           fHistoNonStoppingTracks->Fill(etagen[multMCCharged],vtxMC[2]);
 812           for (Int_t iref=0;iref<nref;iref++) { //since it is detectable, check if it is also detected
 813             tref = mcpart->GetTrackReference(iref);
 814             if (tref) {
 815               if (tref->R()>rminL2&&tref->R()<rmaxL2) {
 816                 if (tref->DetectorId()==0) {
 817                   detectedPrimaryPartLay2[multMCCharged]=kTRUE;
 818                   fHistoDetectedLay2->Fill(etagen[multMCCharged],vtxMC[2]);
 819                   break;
 820                 }
 821               }
 822             }
 823           }
 824         } else { //last is -1 -> particle disappeared. Where?
 825           tref = mcpart->GetTrackReference(nref-1);
 826           fHistoRTRm1->Fill(tref->R());
 827           if (tref->R()>rmaxL2) {
 828             detectablePrimaryPart[multMCCharged]=kTRUE;
 829             fHistoDetectabletr->Fill(etagen[multMCCharged],vtxMC[2]);
 830             for (Int_t iref=0;iref<nref;iref++) { //since it is detectable, check if it is also detected
 831               tref = mcpart->GetTrackReference(iref);
 832               if (tref) {
 833                 if (tref->R()>rminL2&&tref->R()<rmaxL2) {
 834                   if (tref->DetectorId()==0) {
 835                     detectedPrimaryPartLay2[multMCCharged]=kTRUE;
 836                     fHistoDetectedLay2->Fill(etagen[multMCCharged],vtxMC[2]);
 837                     break;
 838                   }
 839                 }
 840               }
 841             }
 842           } else if (tref->R()>=rminL2&&tref->R()<=rmaxL2) {
 843             for (Int_t iref=0;iref<nref;iref++) {
 844               tref = mcpart->GetTrackReference(iref);
 845               if (tref) {
 846                 if (tref->R()>rminL2&&tref->R()<rmaxL2) {
 847                   if (tref->DetectorId()==0) {
 848                     detectablePrimaryPart[multMCCharged]=kTRUE;
 849                     detectedPrimaryPartLay2[multMCCharged]=kTRUE;
 850                     fHistoDetectedLay2->Fill(etagen[multMCCharged],vtxMC[2]);
 851                     fHistoDetectabletr->Fill(etagen[multMCCharged],vtxMC[2]);
 852                     break;
 853                   }
 854                 }
 855               }
 856             }
 857           }
 858         }
 859         // Find out detected prims on each layer
 860         for (Int_t iref=0; iref<nref; iref++) {
 861           tref = mcpart->GetTrackReference(iref);
 862           if (tref->R()>rminL1&&tref->R()<rmaxL1&&tref->DetectorId()==0) {
 863             detectedPrimaryPartLay1[multMCCharged] = kTRUE;
 864             fHistoDetectedLay1->Fill(etagen[multMCCharged],vtxMC[2]);
 865           }
 866         }
 867       }
 868
 869       multMCCharged++;
 870       if (TMath::Abs(eta)<1.4) multMCChargedEtacut++;
 871     } // End of MC particle loop
 872
 873
 874     if (ndEv) { // non diffractive
 875       fHistMultAllNonDiff->Fill(multMCChargedEtacut);
 876       if (eventTriggered&&esdvtx[2]!=0.) fHistMultTrVtxNonDiff->Fill(multMCChargedEtacut);
 877
 878     }
 879     if (sdEv) { // single diffractive
 880       fHistMultAllSingleDiff->Fill(multMCChargedEtacut);
 881       if (eventTriggered&&esdvtx[2]!=0.) fHistMultTrVtxSingleDiff->Fill(multMCChargedEtacut);
 882     }
 883     if (ddEv) { // double diffractive
 884       fHistMultAllDoubleDiff->Fill(multMCChargedEtacut);
 885       if (eventTriggered&&esdvtx[2]!=0.) fHistMultTrVtxDoubleDiff->Fill(multMCChargedEtacut);
 886     }
 887
 888     if (esdvtx[2]==0.) {
 889       fHistContributorsvsMCVtx->Fill(vtxMC[2],nContrib);
 890     }
 891     // Event selection
 892     if (eventTriggered&&eventWithVertex) {
 893
 894       for (Int_t itracklet=0; itracklet<multSPD; ++itracklet) {
 895
 896         Int_t labL1 = multESD->GetLabel(itracklet,0);
 897         Int_t labL2 = multESD->GetLabel(itracklet,1);
 898
 899         fHistBkgCorrDen->Fill(multESD->GetEta(itracklet),esdvtx[2]);
 900
 901         if (labL1==labL2) {
 902           for (Int_t imc=0; imc<multMCCharged; imc++) {
 903             if (labL1==stackIndexOfPrimaryParts[imc]) {
 904               if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) {
 905                 reconstructedPrimaryPart[imc]=kTRUE;
 906                 break;
 907               }
 908             }
 909           }
 910         }
 911       }
 912
 913       for (Int_t imc=0; imc<multMCCharged; imc++) {
 914         if (reconstructedPrimaryPart[imc]) {
 915           fHistBkgCorrNum->Fill(etagen[imc],vtxMC[2]);
 916         }
 917         if (detectedPrimaryPartLay1[imc]&&detectedPrimaryPartLay2[imc]) fHistAlgEffNum->Fill(etagen[imc],vtxMC[2]);
 918
 919         fHistNonDetectableCorrNum->Fill(etagen[imc],vtxMC[2]);
 920
 921         if (detectablePrimaryPart[imc]) fHistNonDetectableCorrDen->Fill(etagen[imc],vtxMC[2]);
 922
 923         fHistMCEtavsZTriggESDvtxEvts->Fill(etagen[imc],esdvtx[2]);
 924         fHistMCEtavsZTriggMCvtxEvts->Fill(etagen[imc],vtxMC[2]);
 925       }
 926
 927       fHistTrigVtxCorrDen->Fill(multSPD,vtxMC[2]);
 928     } // End of selected events
 929
 930     if (eventTriggered) {
 931         fHistTrigCorrDen->Fill(multSPD,vtxMC[2]);
 932         fHistTrigVtxCorrNum->Fill(multSPD,vtxMC[2]);
 933       if (nsdEv) {
 934         fHistEvTrigNSD->Fill(multSPD,vtxMC[2]); //to compute errors
 935         fHistTrigVtxCorrNumNSD->Fill(multSPD,vtxMC[2]);
 936       }
 937       for (Int_t imc=0; imc<multMCCharged; imc++) {
 938         fHistTrackTrigCorrDen->Fill(etagen[imc],vtxMC[2]);
 939         if (nsdEv) fHistTrackTrigNSD->Fill(etagen[imc],vtxMC[2]); //to compute errors
 940       }
 941     } else {
 942        fHistTrigVtxCorrNum->Fill(0.,vtxMC[2]);
 943        if (nsdEv) fHistTrigVtxCorrNumNSD->Fill(0.,vtxMC[2]);
 944
 945     }
 946
 947     for (Int_t imc=0; imc<multMCCharged; imc++) {
 948       fHistTrackTrigVtxCorrNum->Fill(etagen[imc],vtxMC[2]);
 949       fHistMCEtavsZ->Fill(etagen[imc],vtxMC[2]);
 950       if (nsdEv)
 951         fHistTrackTrigVtxCorrNumNSD->Fill(etagen[imc],vtxMC[2]);
 952     }
 953
 954     delete[] etagen;
 955     delete[] ptgen;
 956     delete[] stackIndexOfPrimaryParts;
 957     delete[] reconstructedPrimaryPart;
 958     delete[] detectedPrimaryPartLay1;
 959     delete[] detectedPrimaryPartLay2;
 960     delete[] detectablePrimaryPart;
 961   }
 962 }
 963
 964 //________________________________________________________________________
 965 void AliAnalysisTaskSPDdNdEta::Terminate(Option_t *)
 966 {
 967   // Called once at the end of the query
 968   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 969   if (!fOutput) {     Printf("ERROR: fOutput not available");
 970     return;
 971   }
 972
 973   fHistSPDRAWMultvsZ= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDRAWMultvsZ"));
 974   fHistSPDRAWMultvsZTriggEvts = dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDRAWMultvsZTriggEvts"));
 975   fHistSPDRAWEtavsZ = dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDRAWEtavsZ"));
 976
 977   fHistSPDmultEtacut = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDmultEtacut"));
 978   fHistSPDmult = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDmult"));
 979   fHistSPDeta = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDeta"));
 980   fHistSPDcl1multEtacutLay1 = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDcl1multEtacutLay1"));
 981   fHistSPDcl1mult = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDcl1mult"));
 982   fHistSPDcl1eta = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDcl1eta"));
 983   fHistSPDphi = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDphi"));
 984   fHistSPDcl1phi = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDcl1phi"));
 985   fHistSPDtheta = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDtheta"));
 986   fHistSPDcl1theta = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDcl1theta"));
 987   fHistSPDdePhi = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDdePhi"));
 988   fHistSPDdePhiZ = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDdePhiZ"));
 989   fHistSPDdePhi3D = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDdePhi3D"));
 990   fHistSPDphivsSPDeta= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDphivsSPDeta"));
 991   fHistSPDcl1phivsSPDcl1eta= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDcl1phivsSPDcl1eta"));
 992   fHistSPDdeTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDdeTheta"));
 993
 994   fHistSPDvtxAnalysis = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDvtxAnalysis"));
 995   fHistSPDvtx3D = dynamic_cast<TH3F*> (fOutput->FindObject("fHistSPDvtx3D"));
 996   fHistSPDvtxZ = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDvtxZ"));
 997   fHistNcontribSPDvtxvsSPDvtx = dynamic_cast<TH2F*> (fOutput->FindObject("fHistNcontribSPDvtxvsSPDvtx"));
 998   fHistNcontribSPDvtx3D = dynamic_cast<TH1F*> (fOutput->FindObject("fHistNcontribSPDvtx3D"));
 999   fHistNcontribSPDvtxZ = dynamic_cast<TH1F*> (fOutput->FindObject("fHistNcontribSPDvtxZ"));
1000   fHistNcontribSPDvtxall = dynamic_cast<TH1F*> (fOutput->FindObject("fHistNcontribSPDvtxall"));
1001
1002   fHistSPDcl1multvsnFiredChipsLay1= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDcl1multvsnFiredChipsLay1"));
1003   fHistSPDmultvsnFiredChipsLay1= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDmultvsnFiredChipsLay1"));
1004   fHistSPDmultvsnFiredChipsLay2= dynamic_cast<TH2F*> (fOutput->FindObject("fHistSPDmultvsnFiredChipsLay2"));
1005   fHistnFiredChipsLay2vsnFiredChipsLay1= dynamic_cast<TH2F*> (fOutput->FindObject("fHistnFiredChipsLay2vsnFiredChipsLay1"));
1006   fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec= dynamic_cast<TH2F*> (fOutput->FindObject("fHistnFiredChipsLay2vsnFiredChipsLay1novtxrec"));
1007
1008   fHistSPDvtxRec = dynamic_cast<TH1F*> (fOutput->FindObject("fHistSPDvtxRec"));
1009
1010   if (fCorr) {
1011
1012     fHistBkgCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistBkgCorrNum"));
1013     fHistBkgCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistBkgCorrDen"));
1014
1015     fHistAlgEffNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistAlgEffNum"));
1016
1017     fHistNonDetectableCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistNonDetectableCorrNum"));
1018     fHistNonDetectableCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistNonDetectableCorrDen"));
1019
1020     fHistTrackTrigVtxCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrackTrigVtxCorrNum"));
1021
1022     fHistTrackTrigCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrackTrigCorrDen"));
1023     fHistTrackTrigVtxCorrNumNSD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrackTrigVtxCorrNumNSD"));
1024     fHistTrackTrigNSD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrackTrigNSD"));
1025
1026
1027     fHistTrigVtxCorrNum = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigVtxCorrNum"));
1028     fHistTrigVtxCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigVtxCorrDen"));
1029
1030     fHistTrigCorrDen = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigCorrDen"));
1031
1032     fHistTrigVtxCorrNumNSD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistTrigVtxCorrNumNSD"));
1033     fHistEvTrigNSD = dynamic_cast<TH2F*> (fOutput->FindObject("fHistEvTrigNSD"));
1034
1035     fHistMCEtavsZTriggMCvtxEvts = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZTriggMCvtxEvts"));
1036     fHistMCEtavsZTriggESDvtxEvts = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZTriggESDvtxEvts"));
1037     fHistMCEtavsZ = dynamic_cast<TH2F*> (fOutput->FindObject("fHistMCEtavsZ"));
1038
1039     fHistMCEtaInel = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMCEtaInel"));
1040     fHistMCEtaNonDiffractive = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMCEtaNonDiffractive"));
1041     fHistMCEtaNonSingleDiffractive = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMCEtaNonSingleDiffractive"));
1042
1043     fHistoProcessType = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoProcessType"));
1044     fHistoProcessTypeTriggered = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoProcessTypeTriggered"));
1045
1046     fHistContributorsvsMCVtx = dynamic_cast<TH2F*> (fOutput->FindObject("fHistContributorsvsMCVtx"));
1047
1048     fHistoDetectableNotr = dynamic_cast<TH3F*> (fOutput->FindObject("fHistoDetectableNotr"));
1049     fHistoDetectabletr = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectabletr"));
1050
1051     fHistoNonStoppingTracks = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoNonStoppingTracks"));
1052
1053     fHistoDetectedLay1 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectedLay1"));
1054     fHistoDetectedLay2 = dynamic_cast<TH2F*> (fOutput->FindObject("fHistoDetectedLay2"));
1055
1056     fHistoPt = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoPt"));
1057
1058     fHistoRTRm1 = dynamic_cast<TH1F*> (fOutput->FindObject("fHistoRTRm1"));
1059     fHistMCvtx = dynamic_cast<TH3F*> (fOutput->FindObject("fHistMCvtx"));
1060
1061     fHistMultAllNonDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultAllNonDiff"));
1062     fHistMultAllSingleDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultAllSingleDiff"));
1063     fHistMultAllDoubleDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultAllDoubleDiff"));
1064     fHistMultTrVtxNonDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultTrVtxNonDiff"));
1065     fHistMultTrVtxSingleDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultTrVtxSingleDiff"));
1066     fHistMultTrVtxDoubleDiff = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMultTrVtxDoubleDiff"));
1067
1068     fHistMCEtaNonSingleDiffractiveLargeBin = dynamic_cast<TH1F*> (fOutput->FindObject("fHistMCEtaNonSingleDiffractiveLargeBin"));
1069   }
1070 }