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