PWG0/dNdEta/AlidNdEtaTask.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaTask.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TVector3.h>
   9 #include <TChain.h>
  10 #include <TFile.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TH3F.h>
  14 #include <TParticle.h>
  15 #include <TRandom.h>
  16 #include <TNtuple.h>
  17 #include <TObjString.h>
  18 #include <TF1.h>
  19 #include <TGraph.h>
  20
  21 #include <AliLog.h>
  22 #include <AliESDVertex.h>
  23 #include <AliESDEvent.h>
  24 #include <AliStack.h>
  25 #include <AliHeader.h>
  26 #include <AliGenEventHeader.h>
  27 #include <AliMultiplicity.h>
  28 #include <AliAnalysisManager.h>
  29 #include <AliMCEventHandler.h>
  30 #include <AliMCEvent.h>
  31 #include <AliESDInputHandler.h>
  32 #include <AliESDHeader.h>
  33
  34 #include "AliESDtrackCuts.h"
  35 #include "AliPWG0Helper.h"
  36 #include "AliCorrection.h"
  37 #include "AliCorrectionMatrix3D.h"
  38 #include "dNdEta/dNdEtaAnalysis.h"
  39 #include "AliTriggerAnalysis.h"
  40
  41 ClassImp(AlidNdEtaTask)
  42
  43 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
  44   AliAnalysisTask("AlidNdEtaTask", ""),
  45   fESD(0),
  46   fOutput(0),
  47   fOption(opt),
  48   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
  49   fTrigger(AliTriggerAnalysis::kMB1),
  50   fFillPhi(kFALSE),
  51   fDeltaPhiCut(-1),
  52   fReadMC(kFALSE),
  53   fUseMCVertex(kFALSE),
  54   fOnlyPrimaries(kFALSE),
  55   fUseMCKine(kFALSE),
  56   fEsdTrackCuts(0),
  57   fdNdEtaAnalysisESD(0),
  58   fMult(0),
  59   fMultVtx(0),
  60   fEvents(0),
  61   fVertexResolution(0),
  62   fdNdEtaAnalysis(0),
  63   fdNdEtaAnalysisND(0),
  64   fdNdEtaAnalysisNSD(0),
  65   fdNdEtaAnalysisTr(0),
  66   fdNdEtaAnalysisTrVtx(0),
  67   fdNdEtaAnalysisTracks(0),
  68   fPartPt(0),
  69   fVertex(0),
  70   fPhi(0),
  71   fRawPt(0),
  72   fEtaPhi(0),
  73   fDeltaPhi(0),
  74   fDeltaTheta(0),
  75   fFiredChips(0),
  76   fTriggerVsTime(0),
  77   fStats(0)
  78 {
  79   //
  80   // Constructor. Initialization of pointers
  81   //
  82
  83   // Define input and output slots here
  84   DefineInput(0, TChain::Class());
  85   DefineOutput(0, TList::Class());
  86
  87   fZPhi[0] = 0;
  88   fZPhi[1] = 0;
  89
  90   AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
  91 }
  92
  93 AlidNdEtaTask::~AlidNdEtaTask()
  94 {
  95   //
  96   // Destructor
  97   //
  98
  99   // histograms are in the output list and deleted when the output
 100   // list is deleted by the TSelector dtor
 101
 102   if (fOutput) {
 103     delete fOutput;
 104     fOutput = 0;
 105   }
 106 }
 107
 108 Bool_t AlidNdEtaTask::Notify()
 109 {
 110   static Int_t count = 0;
 111   count++;
 112   Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
 113   return kTRUE;
 114 }
 115
 116 //________________________________________________________________________
 117 void AlidNdEtaTask::ConnectInputData(Option_t *)
 118 {
 119   // Connect ESD
 120   // Called once
 121
 122   Printf("AlidNdEtaTask::ConnectInputData called");
 123
 124   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 125
 126   if (!esdH) {
 127     Printf("ERROR: Could not get ESDInputHandler");
 128   } else {
 129     fESD = esdH->GetEvent();
 130
 131     TString branches("AliESDHeader Vertex ");
 132
 133     if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger & AliTriggerAnalysis::kOfflineFlag)
 134       branches += "AliMultiplicity ";
 135
 136     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 137       branches += "Tracks ";
 138
 139     // Enable only the needed branches
 140     esdH->SetActiveBranches(branches);
 141   }
 142
 143   // disable info messages of AliMCEvent (per event)
 144   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 145 }
 146
 147 void AlidNdEtaTask::CreateOutputObjects()
 148 {
 149   // create result objects and add to output list
 150
 151   Printf("AlidNdEtaTask::CreateOutputObjects");
 152
 153   if (fOnlyPrimaries)
 154     Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
 155
 156   if (fUseMCKine)
 157     Printf("WARNING: Using MC kine information. This is for systematical checks only.");
 158
 159   if (fUseMCVertex)
 160     Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 161
 162   fOutput = new TList;
 163   fOutput->SetOwner();
 164
 165   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
 166   fOutput->Add(fdNdEtaAnalysisESD);
 167
 168   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
 169   fOutput->Add(fMult);
 170
 171   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
 172   fOutput->Add(fMultVtx);
 173
 174   for (Int_t i=0; i<3; ++i)
 175   {
 176     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
 177     fPartEta[i]->Sumw2();
 178     fOutput->Add(fPartEta[i]);
 179   }
 180
 181   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
 182   fOutput->Add(fEvents);
 183
 184   Float_t resMax = (fAnalysisMode & AliPWG0Helper::kSPD) ? 0.2 : 2;
 185   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, resMax);
 186   fOutput->Add(fVertexResolution);
 187
 188   fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
 189   fOutput->Add(fPhi);
 190
 191   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*5, 0, 2 * TMath::Pi());
 192   fOutput->Add(fEtaPhi);
 193
 194   fTriggerVsTime = new TGraph; //TH1F("fTriggerVsTime", "fTriggerVsTime;trigger time;count", 100, 0, 100);
 195   fTriggerVsTime->SetName("fTriggerVsTime");
 196   fTriggerVsTime->GetXaxis()->SetTitle("trigger time");
 197   fTriggerVsTime->GetYaxis()->SetTitle("count");
 198   fOutput->Add(fTriggerVsTime);
 199
 200   fStats = new TH1F("fStats", "fStats", 3, 0.5, 3.5);
 201   fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
 202   fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
 203   fStats->GetXaxis()->SetBinLabel(3, "trigger");
 204   fOutput->Add(fStats);
 205
 206   if (fAnalysisMode & AliPWG0Helper::kSPD)
 207   {
 208     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.2, 0.2);
 209     fOutput->Add(fDeltaPhi);
 210     fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.2, 0.2);
 211     fOutput->Add(fDeltaTheta);
 212     fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
 213     fOutput->Add(fFiredChips);
 214     for (Int_t i=0; i<2; i++)
 215     {
 216       fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -20, 20, 180, 0, TMath::Pi() * 2);
 217       fOutput->Add(fZPhi[i]);
 218     }
 219   }
 220
 221   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 222   {
 223     fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
 224     fOutput->Add(fRawPt);
 225   }
 226
 227   fVertex = new TH3F("vertex_check", "vertex_check", 100, -1, 1, 100, -1, 1, 100, -30, 30);
 228   fOutput->Add(fVertex);
 229
 230   if (fReadMC)
 231   {
 232     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
 233     fOutput->Add(fdNdEtaAnalysis);
 234
 235     fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
 236     fOutput->Add(fdNdEtaAnalysisND);
 237
 238     fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
 239     fOutput->Add(fdNdEtaAnalysisNSD);
 240
 241     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
 242     fOutput->Add(fdNdEtaAnalysisTr);
 243
 244     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
 245     fOutput->Add(fdNdEtaAnalysisTrVtx);
 246
 247     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
 248     fOutput->Add(fdNdEtaAnalysisTracks);
 249
 250     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
 251     fPartPt->Sumw2();
 252     fOutput->Add(fPartPt);
 253   }
 254
 255   if (fEsdTrackCuts)
 256   {
 257     fEsdTrackCuts->SetName("fEsdTrackCuts");
 258     fOutput->Add(fEsdTrackCuts);
 259   }
 260 }
 261
 262 void AlidNdEtaTask::Exec(Option_t*)
 263 {
 264   // process the event
 265
 266   // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
 267   Bool_t eventTriggered = kFALSE;
 268   const AliESDVertex* vtxESD = 0;
 269
 270   // post the data already here
 271   PostData(0, fOutput);
 272
 273   // ESD analysis
 274   if (fESD)
 275   {
 276     // check event type (should be PHYSICS = 7)
 277     AliESDHeader* esdHeader = fESD->GetHeader();
 278     if (!esdHeader)
 279     {
 280       Printf("ERROR: esdHeader could not be retrieved");
 281       return;
 282     }
 283
 284     Printf("Trigger classes: %s:", fESD->GetFiredTriggerClasses().Data());
 285
 286 //    UInt_t eventType = esdHeader->GetEventType();
 287 //     if (eventType != 7)
 288 //     {
 289 //       Printf("Skipping event because it is of type %d", eventType);
 290 //       return;
 291 //     }
 292
 293     // trigger definition
 294     static AliTriggerAnalysis* triggerAnalysis = new AliTriggerAnalysis;
 295     eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
 296     if (eventTriggered)
 297       fStats->Fill(3);
 298
 299     // get the ESD vertex
 300     vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 301
 302     //Printf("vertex is: %s", fESD->GetPrimaryVertex()->GetTitle());
 303
 304     Double_t vtx[3];
 305
 306     // fill z vertex resolution
 307     if (vtxESD)
 308     {
 309       fVertexResolution->Fill(vtxESD->GetZRes());
 310       fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
 311
 312       if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 313       {
 314           vtxESD->GetXYZ(vtx);
 315
 316           // vertex stats
 317           if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
 318           {
 319             fStats->Fill(1);
 320           }
 321           else if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
 322             fStats->Fill(2);
 323       }
 324       else
 325         vtxESD = 0;
 326     }
 327
 328     // needed for syst. studies
 329     AliStack* stack = 0;
 330     TArrayF vtxMC(3);
 331
 332     if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
 333       if (!fReadMC) {
 334         Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
 335         return;
 336       }
 337
 338       AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 339       if (!eventHandler) {
 340         Printf("ERROR: Could not retrieve MC event handler");
 341         return;
 342       }
 343
 344       AliMCEvent* mcEvent = eventHandler->MCEvent();
 345       if (!mcEvent) {
 346         Printf("ERROR: Could not retrieve MC event");
 347         return;
 348       }
 349
 350       AliHeader* header = mcEvent->Header();
 351       if (!header)
 352       {
 353         AliDebug(AliLog::kError, "Header not available");
 354         return;
 355       }
 356
 357       // get the MC vertex
 358       AliGenEventHeader* genHeader = header->GenEventHeader();
 359       if (!genHeader)
 360       {
 361         AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 362         return;
 363       }
 364       genHeader->PrimaryVertex(vtxMC);
 365
 366       if (fUseMCVertex)
 367         vtx[2] = vtxMC[2];
 368
 369       stack = mcEvent->Stack();
 370       if (!stack)
 371       {
 372         AliDebug(AliLog::kError, "Stack not available");
 373         return;
 374       }
 375     }
 376
 377     // create list of (label, eta, pt) tuples
 378     Int_t inputCount = 0;
 379     Int_t* labelArr = 0;
 380     Float_t* etaArr = 0;
 381     Float_t* thirdDimArr = 0;
 382     if (fAnalysisMode & AliPWG0Helper::kSPD)
 383     {
 384       // get tracklets
 385       const AliMultiplicity* mult = fESD->GetMultiplicity();
 386       if (!mult)
 387       {
 388         AliDebug(AliLog::kError, "AliMultiplicity not available");
 389         return;
 390       }
 391
 392       labelArr = new Int_t[mult->GetNumberOfTracklets()];
 393       etaArr = new Float_t[mult->GetNumberOfTracklets()];
 394       thirdDimArr = new Float_t[mult->GetNumberOfTracklets()];
 395
 396       // get multiplicity from SPD tracklets
 397       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 398       {
 399         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 400
 401         if (fOnlyPrimaries)
 402           if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
 403             continue;
 404
 405         Float_t deltaPhi = mult->GetDeltaPhi(i);
 406         // prevent values to be shifted by 2 Pi()
 407         if (deltaPhi < -TMath::Pi())
 408           deltaPhi += TMath::Pi() * 2;
 409         if (deltaPhi > TMath::Pi())
 410           deltaPhi -= TMath::Pi() * 2;
 411
 412         if (TMath::Abs(deltaPhi) > 1)
 413           printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
 414
 415         Int_t label = mult->GetLabel(i, 0);
 416         Float_t eta = mult->GetEta(i);
 417
 418         // control histograms
 419         Float_t phi = mult->GetPhi(i);
 420         if (phi < 0)
 421           phi += TMath::Pi() * 2;
 422         fPhi->Fill(phi);
 423         fEtaPhi->Fill(eta, phi);
 424
 425         if (deltaPhi < 0.01)
 426         {
 427           // layer 0
 428           Float_t z = vtx[2] + 3.9 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
 429           fZPhi[0]->Fill(z, phi);
 430           // layer 1
 431           z = vtx[2] + 7.6 / TMath::Tan(2 * TMath::ATan(TMath::Exp(- eta)));
 432           fZPhi[1]->Fill(z, phi);
 433         }
 434
 435         fDeltaPhi->Fill(deltaPhi);
 436         fDeltaTheta->Fill(mult->GetDeltaTheta(i));
 437
 438         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
 439           continue;
 440
 441         if (fUseMCKine)
 442         {
 443           if (label > 0)
 444           {
 445             TParticle* particle = stack->Particle(label);
 446             eta = particle->Eta();
 447             phi = particle->Phi();
 448           }
 449           else
 450             Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
 451         }
 452
 453         etaArr[inputCount] = eta;
 454         labelArr[inputCount] = label;
 455         thirdDimArr[inputCount] = phi;
 456         ++inputCount;
 457       }
 458
 459       if (!fFillPhi)
 460       {
 461         // fill multiplicity in third axis
 462         for (Int_t i=0; i<inputCount; ++i)
 463           thirdDimArr[i] = inputCount;
 464       }
 465
 466       Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
 467       fFiredChips->Fill(firedChips, inputCount);
 468       Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
 469     }
 470     else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
 471     {
 472       if (!fEsdTrackCuts)
 473       {
 474         AliDebug(AliLog::kError, "fESDTrackCuts not available");
 475         return;
 476       }
 477
 478       if (vtxESD)
 479       {
 480         // get multiplicity from ESD tracks
 481         TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
 482         Int_t nGoodTracks = list->GetEntries();
 483         Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
 484
 485         labelArr = new Int_t[nGoodTracks];
 486         etaArr = new Float_t[nGoodTracks];
 487         thirdDimArr = new Float_t[nGoodTracks];
 488
 489         // loop over esd tracks
 490         for (Int_t i=0; i<nGoodTracks; i++)
 491         {
 492           AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 493           if (!esdTrack)
 494           {
 495             AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 496             continue;
 497           }
 498
 499           Float_t phi = esdTrack->Phi();
 500           if (phi < 0)
 501             phi += TMath::Pi() * 2;
 502
 503           Float_t eta = esdTrack->Eta();
 504           Int_t label = TMath::Abs(esdTrack->GetLabel());
 505           Float_t pT  = esdTrack->Pt();
 506
 507           // force pT to fixed value without B field
 508           if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
 509             pT = 1;
 510
 511           fPhi->Fill(phi);
 512           fEtaPhi->Fill(eta, phi);
 513           if (eventTriggered && vtxESD)
 514             fRawPt->Fill(pT);
 515
 516           if (fOnlyPrimaries)
 517           {
 518             if (label == 0)
 519               continue;
 520
 521             if (stack->IsPhysicalPrimary(label) == kFALSE)
 522               continue;
 523           }
 524
 525           if (fUseMCKine)
 526           {
 527             if (label > 0)
 528             {
 529               TParticle* particle = stack->Particle(label);
 530               eta = particle->Eta();
 531               pT = particle->Pt();
 532             }
 533             else
 534               Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
 535           }
 536
 537           etaArr[inputCount] = eta;
 538           labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 539           thirdDimArr[inputCount] = pT;
 540           ++inputCount;
 541         }
 542
 543         // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
 544
 545         delete list;
 546       }
 547     }
 548     else
 549       return;
 550
 551     // Processing of ESD information (always)
 552     if (eventTriggered)
 553     {
 554       // control hist
 555       fMult->Fill(inputCount);
 556       fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
 557
 558       fTriggerVsTime->SetPoint(fTriggerVsTime->GetN(), fESD->GetTimeStamp(), 1);
 559
 560       if (vtxESD)
 561       {
 562         // control hist
 563         fMultVtx->Fill(inputCount);
 564
 565         for (Int_t i=0; i<inputCount; ++i)
 566         {
 567           Float_t eta = etaArr[i];
 568           Float_t thirdDim = thirdDimArr[i];
 569
 570           fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
 571
 572           if (TMath::Abs(vtx[2]) < 10)
 573           {
 574             fPartEta[0]->Fill(eta);
 575
 576             if (vtx[2] < 0)
 577               fPartEta[1]->Fill(eta);
 578             else
 579               fPartEta[2]->Fill(eta);
 580           }
 581         }
 582
 583         // for event count per vertex
 584         fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
 585
 586         // control hist
 587         fEvents->Fill(vtx[2]);
 588
 589         if (fReadMC)
 590         {
 591           // from tracks is only done for triggered and vertex reconstructed events
 592           for (Int_t i=0; i<inputCount; ++i)
 593           {
 594             Int_t label = labelArr[i];
 595
 596             if (label < 0)
 597               continue;
 598
 599             //Printf("Getting particle of track %d", label);
 600             TParticle* particle = stack->Particle(label);
 601
 602             if (!particle)
 603             {
 604               AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
 605               continue;
 606             }
 607
 608             Float_t thirdDim = -1;
 609             if (fAnalysisMode & AliPWG0Helper::kSPD)
 610             {
 611               if (fFillPhi)
 612               {
 613                 thirdDim = particle->Phi();
 614               }
 615               else
 616                 thirdDim = inputCount;
 617             }
 618             else
 619               thirdDim = particle->Pt();
 620
 621             fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), thirdDim);
 622           } // end of track loop
 623
 624           // for event count per vertex
 625           fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
 626         }
 627       }
 628     }
 629
 630     if (etaArr)
 631       delete[] etaArr;
 632     if (labelArr)
 633       delete[] labelArr;
 634     if (thirdDimArr)
 635       delete[] thirdDimArr;
 636   }
 637
 638   if (fReadMC)   // Processing of MC information (optional)
 639   {
 640     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 641     if (!eventHandler) {
 642       Printf("ERROR: Could not retrieve MC event handler");
 643       return;
 644     }
 645
 646     AliMCEvent* mcEvent = eventHandler->MCEvent();
 647     if (!mcEvent) {
 648       Printf("ERROR: Could not retrieve MC event");
 649       return;
 650     }
 651
 652     AliStack* stack = mcEvent->Stack();
 653     if (!stack)
 654     {
 655       AliDebug(AliLog::kError, "Stack not available");
 656       return;
 657     }
 658
 659     AliHeader* header = mcEvent->Header();
 660     if (!header)
 661     {
 662       AliDebug(AliLog::kError, "Header not available");
 663       return;
 664     }
 665
 666     // get the MC vertex
 667     AliGenEventHeader* genHeader = header->GenEventHeader();
 668     if (!genHeader)
 669     {
 670       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 671       return;
 672     }
 673
 674     TArrayF vtxMC(3);
 675     genHeader->PrimaryVertex(vtxMC);
 676
 677     // get process type
 678     Int_t processType = AliPWG0Helper::GetEventProcessType(header);
 679     AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
 680
 681     if (processType==AliPWG0Helper::kInvalidProcess)
 682       AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
 683
 684     // loop over mc particles
 685     Int_t nPrim  = stack->GetNprimary();
 686
 687     Int_t nAcceptedParticles = 0;
 688
 689     // count particles first, then fill
 690     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 691     {
 692       if (!stack->IsPhysicalPrimary(iMc))
 693         continue;
 694
 695       //Printf("Getting particle %d", iMc);
 696       TParticle* particle = stack->Particle(iMc);
 697
 698       if (!particle)
 699         continue;
 700
 701       if (particle->GetPDG()->Charge() == 0)
 702         continue;
 703
 704       //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 705       Float_t eta = particle->Eta();
 706
 707        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
 708       if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
 709         nAcceptedParticles++;
 710     }
 711
 712     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 713     {
 714       if (!stack->IsPhysicalPrimary(iMc))
 715         continue;
 716
 717       //Printf("Getting particle %d", iMc);
 718       TParticle* particle = stack->Particle(iMc);
 719
 720       if (!particle)
 721         continue;
 722
 723       if (particle->GetPDG()->Charge() == 0)
 724         continue;
 725
 726       Float_t eta = particle->Eta();
 727       Float_t thirdDim = -1;
 728
 729       if (fAnalysisMode & AliPWG0Helper::kSPD)
 730       {
 731         if (fFillPhi)
 732         {
 733           thirdDim = particle->Phi();
 734         }
 735         else
 736           thirdDim = nAcceptedParticles;
 737       }
 738       else
 739         thirdDim = particle->Pt();
 740
 741       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
 742
 743       if (processType != AliPWG0Helper::kSD)
 744         fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
 745
 746       if (processType == AliPWG0Helper::kND)
 747         fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
 748
 749       if (eventTriggered)
 750       {
 751         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
 752         if (vtxESD)
 753           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
 754       }
 755
 756       if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
 757       {
 758         Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
 759         fPartPt->Fill(particle->Pt(), value);
 760       }
 761     }
 762
 763     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
 764     if (processType != AliPWG0Helper::kSD)
 765       fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
 766     if (processType == AliPWG0Helper::kND)
 767       fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
 768
 769     if (eventTriggered)
 770     {
 771       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
 772       if (vtxESD)
 773         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
 774     }
 775   }
 776 }
 777
 778 void AlidNdEtaTask::Terminate(Option_t *)
 779 {
 780   // The Terminate() function is the last function to be called during
 781   // a query. It always runs on the client, it can be used to present
 782   // the results graphically or save the results to file.
 783
 784   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 785   if (!fOutput)
 786     Printf("ERROR: fOutput not available");
 787
 788   if (fOutput)
 789   {
 790     fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
 791     fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
 792     fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
 793     for (Int_t i=0; i<3; ++i)
 794       fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 795     fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
 796     fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
 797
 798     fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
 799     fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
 800     fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
 801     for (Int_t i=0; i<2; ++i)
 802       fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
 803     fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
 804     fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
 805     fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
 806     fTriggerVsTime = dynamic_cast<TGraph*> (fOutput->FindObject("fTriggerVsTime"));
 807     fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
 808
 809     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
 810   }
 811
 812   if (!fdNdEtaAnalysisESD)
 813   {
 814     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
 815   }
 816   else
 817   {
 818     if (fMult && fMultVtx)
 819     {
 820       new TCanvas;
 821       fMult->Draw();
 822       fMultVtx->SetLineColor(2);
 823       fMultVtx->Draw("SAME");
 824     }
 825
 826     if (fFiredChips)
 827     {
 828       new TCanvas;
 829       fFiredChips->Draw("COLZ");
 830     }
 831
 832     if (fPartEta[0])
 833     {
 834       Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
 835       Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
 836
 837       Printf("%d events with vertex used", events1 + events2);
 838
 839       if (events1 > 0 && events2 > 0)
 840       {
 841         fPartEta[0]->Scale(1.0 / (events1 + events2));
 842         fPartEta[1]->Scale(1.0 / events1);
 843         fPartEta[2]->Scale(1.0 / events2);
 844
 845         for (Int_t i=0; i<3; ++i)
 846           fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 847
 848         new TCanvas("control", "control", 500, 500);
 849         for (Int_t i=0; i<3; ++i)
 850         {
 851           fPartEta[i]->SetLineColor(i+1);
 852           fPartEta[i]->Draw((i==0) ? "" : "SAME");
 853         }
 854       }
 855     }
 856
 857     if (fEvents)
 858     {
 859         new TCanvas("control3", "control3", 500, 500);
 860         fEvents->Draw();
 861     }
 862
 863     TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
 864
 865     if (fdNdEtaAnalysisESD)
 866       fdNdEtaAnalysisESD->SaveHistograms();
 867
 868     if (fEsdTrackCuts)
 869       fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 870
 871     if (fMult)
 872       fMult->Write();
 873
 874     if (fMultVtx)
 875       fMultVtx->Write();
 876
 877     for (Int_t i=0; i<3; ++i)
 878       if (fPartEta[i])
 879         fPartEta[i]->Write();
 880
 881     if (fEvents)
 882       fEvents->Write();
 883
 884     if (fVertexResolution)
 885       fVertexResolution->Write();
 886
 887     if (fDeltaPhi)
 888       fDeltaPhi->Write();
 889
 890     if (fDeltaTheta)
 891       fDeltaTheta->Write();
 892
 893     if (fPhi)
 894       fPhi->Write();
 895
 896     if (fRawPt)
 897       fRawPt->Write();
 898
 899     if (fEtaPhi)
 900       fEtaPhi->Write();
 901
 902     for (Int_t i=0; i<2; ++i)
 903       if (fZPhi[i])
 904         fZPhi[i]->Write();
 905
 906     if (fFiredChips)
 907       fFiredChips->Write();
 908
 909     if (fTriggerVsTime)
 910       fTriggerVsTime->Write();
 911
 912     if (fStats)
 913       fStats->Write();
 914
 915     fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
 916     if (fVertex)
 917       fVertex->Write();
 918
 919     fout->Write();
 920     fout->Close();
 921
 922     Printf("Writting result to analysis_esd_raw.root");
 923   }
 924
 925   if (fReadMC)
 926   {
 927     if (fOutput)
 928     {
 929       fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 930       fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
 931       fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
 932       fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
 933       fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
 934       fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
 935       fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
 936     }
 937
 938     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
 939     {
 940       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
 941       return;
 942     }
 943
 944     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
 945     fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
 946     fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
 947     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
 948     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
 949     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
 950
 951     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
 952     fPartPt->Scale(1.0/events);
 953     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 954
 955     TFile* fout = new TFile("analysis_mc.root","RECREATE");
 956
 957     fdNdEtaAnalysis->SaveHistograms();
 958     fdNdEtaAnalysisND->SaveHistograms();
 959     fdNdEtaAnalysisNSD->SaveHistograms();
 960     fdNdEtaAnalysisTr->SaveHistograms();
 961     fdNdEtaAnalysisTrVtx->SaveHistograms();
 962     fdNdEtaAnalysisTracks->SaveHistograms();
 963
 964     if (fPartPt)
 965       fPartPt->Write();
 966
 967     fout->Write();
 968     fout->Close();
 969
 970     Printf("Writting result to analysis_mc.root");
 971
 972     if (fPartPt)
 973     {
 974       new TCanvas("control2", "control2", 500, 500);
 975       fPartPt->Draw();
 976     }
 977   }
 978 }