PWG0/multiplicity/AliMultiplicityTask.cxx

   1 /* $Id$ */
   2
   3 #include "AliMultiplicityTask.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 <TH1D.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
  20 #include <AliLog.h>
  21 #include <AliESDVertex.h>
  22 #include <AliESDEvent.h>
  23 #include <AliStack.h>
  24 #include <AliHeader.h>
  25 #include <AliGenEventHeader.h>
  26 #include <AliMultiplicity.h>
  27 #include <AliAnalysisManager.h>
  28 #include <AliMCEventHandler.h>
  29 #include <AliMCEvent.h>
  30 #include <AliESDInputHandler.h>
  31
  32 #include "AliESDtrackCuts.h"
  33 #include "AliPWG0Helper.h"
  34 #include "multiplicity/AliMultiplicityCorrection.h"
  35 #include "AliCorrection.h"
  36 #include "AliCorrectionMatrix3D.h"
  37 #include "AliPhysicsSelection.h"
  38 #include "AliTriggerAnalysis.h"
  39 #include "AliVEvent.h"
  40
  41 ClassImp(AliMultiplicityTask)
  42
  43 AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
  44   AliAnalysisTask("AliMultiplicityTask", ""),
  45   fESD(0),
  46   fOption(opt),
  47   fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
  48   fTrigger(AliTriggerAnalysis::kMB1),
  49   fDeltaPhiCut(-1),
  50   fDiffTreatment(AliPWG0Helper::kMCFlags),
  51   fReadMC(kFALSE),
  52   fUseMCVertex(kFALSE),
  53   fMultiplicity(0),
  54   fEsdTrackCuts(0),
  55   fSystSkipParticles(kFALSE),
  56   fSelectProcessType(0),
  57   fParticleSpecies(0),
  58   fdNdpT(0),
  59   fPtSpectrum(0),
  60   fTemp1(0),
  61   fTemp2(0),
  62   fVertex(0),
  63   fEtaPhi(0),
  64   fOutput(0)
  65 {
  66   //
  67   // Constructor. Initialization of pointers
  68   //
  69
  70   for (Int_t i = 0; i<8; i++)
  71     fParticleCorrection[i] = 0;
  72
  73   for (Int_t i=0; i<3; i++)
  74     fEta[i] = 0;
  75
  76   // Define input and output slots here
  77   DefineInput(0, TChain::Class());
  78   DefineOutput(0, TList::Class());
  79
  80   if (fOption.Contains("only-process-type-nd"))
  81     fSelectProcessType = 1;
  82
  83   if (fOption.Contains("only-process-type-sd"))
  84     fSelectProcessType = 2;
  85
  86   if (fOption.Contains("only-process-type-dd"))
  87     fSelectProcessType = 3;
  88
  89   if (fSelectProcessType != 0)
  90     AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
  91 }
  92
  93 AliMultiplicityTask::~AliMultiplicityTask()
  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&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 103     delete fOutput;
 104     fOutput = 0;
 105   }
 106 }
 107
 108 //________________________________________________________________________
 109 void AliMultiplicityTask::ConnectInputData(Option_t *)
 110 {
 111   // Connect ESD
 112   // Called once
 113
 114   Printf("AliMultiplicityTask::ConnectInputData called");
 115
 116   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 117   if (!tree) {
 118     Printf("ERROR: Could not read tree from input slot 0");
 119   } else {
 120     // Disable all branches and enable only the needed ones
 121     /*
 122     tree->SetBranchStatus("*", 0);
 123
 124     tree->SetBranchStatus("AliESDHeader*", 1);
 125     tree->SetBranchStatus("*Vertex*", 1);
 126
 127     if (fAnalysisMode & AliPWG0Helper::kSPD) {
 128       tree->SetBranchStatus("AliMultiplicity*", 1);
 129     }
 130
 131     if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
 132       //AliESDtrackCuts::EnableNeededBranches(tree);
 133       tree->SetBranchStatus("Tracks*", 1);
 134     }
 135     */
 136
 137     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 138
 139     if (!esdH) {
 140       Printf("ERROR: Could not get ESDInputHandler");
 141     } else
 142       fESD = esdH->GetEvent();
 143   }
 144
 145   // disable info messages of AliMCEvent (per event)
 146   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 147 }
 148
 149 void AliMultiplicityTask::CreateOutputObjects()
 150 {
 151   // create result objects and add to output list
 152
 153   fOutput = new TList;
 154   fOutput->SetOwner();
 155
 156   fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
 157   fOutput->Add(fMultiplicity);
 158
 159   fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
 160   fdNdpT->Sumw2();
 161   fOutput->Add(fdNdpT);
 162
 163   if (fOption.Contains("particle-efficiency"))
 164     for (Int_t i = 0; i<8; i++)
 165     {
 166       fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
 167       fOutput->Add(fParticleCorrection[i]);
 168     }
 169
 170   if (fOption.Contains("pt-spectrum-hist"))
 171   {
 172     TFile* file = TFile::Open("ptspectrum_fit.root");
 173     if (file)
 174     {
 175       TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
 176       TString histName(Form("ptspectrum_%s", subStr.Data()));
 177       AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
 178       fPtSpectrum = (TH1D*) file->Get(histName);
 179       if (!fPtSpectrum)
 180         AliError("Histogram not found");
 181     }
 182     else
 183       AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
 184   }
 185
 186   if (fOption.Contains("pt-spectrum-func"))
 187   {
 188     if (fPtSpectrum)
 189     {
 190       Printf("Using function for systematic p_t study");
 191     }
 192     else
 193     {
 194       Printf("ERROR: Could not find function for systematic p_t study");
 195       fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
 196       fPtSpectrum->SetBinContent(1, 1);
 197     }
 198   }
 199
 200   if (fPtSpectrum)
 201     Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
 202
 203   if (fOption.Contains("particle-species")) {
 204     fParticleSpecies = new TNtuple("fParticleSpecies", "fParticleSpecies", "vtx:Pi_True:K_True:p_True:o_True:Pi_Rec:K_Rec:p_Rec:o_Rec:nolabel:doublePrim:doubleCount");
 205     fOutput->Add(fParticleSpecies);
 206   }
 207
 208   fTemp1 = new TH2F("fTemp1", "fTemp1", 100, -0.5, 99.5, 100, -0.5, 99.5);
 209   fOutput->Add(fTemp1);
 210
 211   for (Int_t i=0; i<3; i++)
 212   {
 213     fEta[i] = new TH1F(Form("fEta_%d", i), ";#eta", 100, -2, 2);
 214     fOutput->Add(fEta[i]);
 215   }
 216
 217   fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
 218   fOutput->Add(fVertex);
 219
 220   fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
 221   fOutput->Add(fEtaPhi);
 222
 223   // TODO set seed for random generator
 224 }
 225
 226 void AliMultiplicityTask::Exec(Option_t*)
 227 {
 228   // process the event
 229
 230   // Check prerequisites
 231   if (!fESD)
 232   {
 233     AliDebug(AliLog::kError, "ESD not available");
 234     return;
 235   }
 236
 237   AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 238   if (!inputHandler)
 239   {
 240     Printf("ERROR: Could not receive input handler");
 241     return;
 242   }
 243
 244   Bool_t eventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
 245
 246   static AliTriggerAnalysis* triggerAnalysis = 0;
 247   if (!triggerAnalysis)
 248   {
 249     AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
 250     triggerAnalysis = physicsSelection->GetTriggerAnalysis();
 251   }
 252   if (eventTriggered)
 253     eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
 254
 255   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
 256   if (vtxESD && !AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
 257     vtxESD = 0;
 258
 259   // remove vertices outside +- 15 cm
 260   if (vtxESD && TMath::Abs(vtxESD->GetZv()) > 15)
 261     vtxESD = 0;
 262
 263   Bool_t eventVertex = (vtxESD != 0);
 264
 265   Double_t vtx[3];
 266   if (vtxESD)
 267   {
 268     vtxESD->GetXYZ(vtx);
 269     fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
 270   }
 271
 272   // post the data already here
 273   PostData(0, fOutput);
 274
 275   //const Float_t kPtCut = 0.3;
 276
 277   // create list of (label, eta) tuples
 278   Int_t inputCount = 0;
 279   Int_t* labelArr = 0;
 280   Float_t* etaArr = 0;
 281   Int_t nGoodTracks = 0;  // number of total good tracks is needed both in SPD and TPC loops if the mode is kTPCSPD
 282   Int_t nESDTracks  = 0; // Total number of ESD tracks (including those not passing the cuts);
 283   const int kRejBit = BIT(15); // set this bit in ESD tracks if it is rejected by a cut
 284
 285   if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)
 286   {
 287     if (!fEsdTrackCuts)
 288     {
 289       AliDebug(AliLog::kError, "fESDTrackCuts not available");
 290       return;
 291     }
 292
 293     if (vtxESD)
 294     {
 295       // get multiplicity from ESD tracks
 296       TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode & AliPWG0Helper::kTPC));
 297       nGoodTracks = list->GetEntries();
 298
 299       Int_t arraySize = nGoodTracks; // if we're also using tracklets, we need to increase this
 300       if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 301         // if I have to also count clusters later on, I have to make sure the arrays are big enough
 302         // Moreover, we need to keep track of rejected tracks to avoid double-counting
 303         Printf( "Processing tracks");
 304         const AliMultiplicity* mult = fESD->GetMultiplicity();
 305         if (!mult)
 306           {
 307             AliDebug(AliLog::kError, "AliMultiplicity not available");
 308             return;
 309           }
 310         arraySize += mult->GetNumberOfTracklets();
 311         // allocate and fill array for rejected tracks
 312         nESDTracks = fESD->GetNumberOfTracks();
 313         for(Int_t itracks = 0; itracks < nESDTracks; itracks++){
 314           AliESDtrack* track = fESD->GetTrack(itracks);
 315           if( itracks!=track->GetID() ){
 316             AliFatal("Track ID not matching track index!");
 317           }
 318           if (!fEsdTrackCuts->AcceptTrack(track)) {
 319             track->SetBit(kRejBit);
 320           }
 321         }
 322       }
 323
 324       labelArr = new Int_t[arraySize];
 325       etaArr = new Float_t[arraySize];
 326
 327       // loop over esd tracks
 328       for (Int_t i=0; i<nGoodTracks; i++)
 329       {
 330         AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
 331         if (!esdTrack)
 332         {
 333           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
 334           continue;
 335         }
 336
 337         if (esdTrack->Pt() < 0.15){
 338           esdTrack->SetBit(kRejBit);
 339           continue;
 340         }
 341
 342         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 343           // Cuts by ruben to flag secondaries
 344           if (esdTrack->IsOn(AliESDtrack::kMultSec) ) continue;
 345           if (esdTrack->IsOn(AliESDtrack::kMultInV0)) continue;
 346         }
 347
 348         Float_t d0z0[2],covd0z0[3];
 349         esdTrack->GetImpactParameters(d0z0,covd0z0);
 350         Float_t sigma= 0.0050+0.0060/TMath::Power(esdTrack->Pt(),0.9);
 351         Float_t d0max = 7.*sigma;
 352         if (TMath::Abs(d0z0[0]) > d0max) {
 353           //      rejLabelArr[irejCount++] = esdTrack->GetID(); // We
 354           //      don't count the tracklet if the track is a
 355           //      secondary, so this must be commented out
 356           continue;
 357         }
 358
 359         if (vtxESD && TMath::Abs(vtx[2]) < 10)
 360           fEtaPhi->Fill(esdTrack->Eta(), esdTrack->Phi());
 361
 362         etaArr[inputCount] = esdTrack->Eta();
 363         labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
 364         ++inputCount;
 365       }
 366
 367       delete list;
 368     }
 369   }
 370   // SPD tracklets are used either in SPD standalone mode or in TPC+SPD mode (tracks + tracklets not using any cluster already used in tracks)
 371   if (fAnalysisMode & AliPWG0Helper::kSPD || fAnalysisMode & AliPWG0Helper::kTPCSPD)
 372   {
 373
 374     if (vtxESD)
 375     {
 376       AliDebug(AliLog::kInfo, "Processing tracklets");
 377       // get tracklets
 378       const AliMultiplicity* mult = fESD->GetMultiplicity();
 379       if (!mult)
 380       {
 381         AliDebug(AliLog::kError, "AliMultiplicity not available");
 382         return;
 383       }
 384
 385       Bool_t foundInEta10 = kFALSE;
 386       if (!(fAnalysisMode & AliPWG0Helper::kTPCSPD)) {
 387         // if we are counting both tracks and tracklets, these arrays were already initialized above
 388         AliDebug(AliLog::kInfo, "Booking arrays");
 389         labelArr = new Int_t[mult->GetNumberOfTracklets()];
 390         etaArr = new Float_t[mult->GetNumberOfTracklets()];
 391       }
 392       if (inputCount) foundInEta10 = kTRUE; // by definition, if we found a track.
 393
 394       // get multiplicity from ITS tracklets
 395       for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
 396       {
 397         //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
 398
 399         Float_t deltaPhi = mult->GetDeltaPhi(i);
 400
 401         if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
 402           continue;
 403
 404         if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
 405           {
 406           Printf("Skipped tracklet!");
 407           continue;
 408           }
 409
 410         // if counting tracks+tracklets, check if clusters where already used in tracks
 411         if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
 412             Int_t id1,id2;
 413             mult->GetTrackletTrackIDs(i,0,id1,id2);
 414             if ( (id1>=0&& !fESD->GetTrack(id1)->TestBit(kRejBit) ) ||
 415                  (id2>=0&& !fESD->GetTrack(id2)->TestBit(kRejBit) )
 416                  ) {
 417               printf ("Skipping tracklet: already used in track");
 418               continue;
 419             }
 420         // Ruben also had this, but we're not using pure ITSSA tracks here:
 421         // if (fSPDMult->FreeClustersTracklet(itr,1)) trITSSApure++; // not used in ITS_SA_Pure track
 422         }
 423
 424
 425         etaArr[inputCount] = mult->GetEta(i);
 426         if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
 427         {
 428           labelArr[inputCount] = mult->GetLabel(i, 0);
 429         }
 430         else
 431           labelArr[inputCount] = -1;
 432
 433         for (Int_t j=0; j<3; j++)
 434         {
 435           if (vtx[2] > fMultiplicity->GetVertexBegin(j) && vtx[2] < fMultiplicity->GetVertexEnd(j))
 436             fEta[j]->Fill(etaArr[inputCount]);
 437         }
 438
 439         if (vtxESD && TMath::Abs(vtx[2]) < 10)
 440           fEtaPhi->Fill(etaArr[inputCount], mult->GetPhi(i));
 441
 442           // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
 443         if (TMath::Abs(etaArr[inputCount]) < 1)
 444           foundInEta10 = kTRUE;
 445
 446         ++inputCount;
 447       }
 448
 449       if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
 450         eventTriggered = kFALSE;
 451     }
 452   }
 453
 454   // eta range for nMCTracksSpecies, nESDTracksSpecies
 455   Float_t etaRange = 1.0;
 456 //   switch (fAnalysisMode) {
 457 //     case AliPWG0Helper::kInvalid: break;
 458 //     case AliPWG0Helper::kTPC:
 459 //     case AliPWG0Helper::kTPCITS:
 460 //      etaRange = 1.0; break;
 461 //     case AliPWG0Helper::kSPD: etaRange = 1.0; break;
 462 //     default: break;
 463 //   }
 464
 465   if (!fReadMC) // Processing of ESD information
 466   {
 467     Int_t nESDTracks05 = 0;
 468     Int_t nESDTracks10 = 0;
 469     Int_t nESDTracks14 = 0;
 470
 471     for (Int_t i=0; i<inputCount; ++i)
 472     {
 473       Float_t eta = etaArr[i];
 474
 475       if (TMath::Abs(eta) < 0.5)
 476         nESDTracks05++;
 477
 478       if (TMath::Abs(eta) < 1.0)
 479         nESDTracks10++;
 480
 481       if (TMath::Abs(eta) < 1.3)
 482         nESDTracks14++;
 483     }
 484
 485     //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
 486     //  Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d; Tracks: %d %d %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber(), nESDTracks05, nESDTracks10, nESDTracks14);
 487
 488     if (eventTriggered)
 489       fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
 490
 491     if (eventTriggered && eventVertex)
 492       fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 493   }
 494   else if (fReadMC)   // Processing of MC information
 495   {
 496     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 497     if (!eventHandler) {
 498       Printf("ERROR: Could not retrieve MC event handler");
 499       return;
 500     }
 501
 502     AliMCEvent* mcEvent = eventHandler->MCEvent();
 503     if (!mcEvent) {
 504       Printf("ERROR: Could not retrieve MC event");
 505       return;
 506     }
 507
 508     AliStack* stack = mcEvent->Stack();
 509     if (!stack)
 510     {
 511       AliDebug(AliLog::kError, "Stack not available");
 512       return;
 513     }
 514
 515     AliHeader* header = mcEvent->Header();
 516     if (!header)
 517     {
 518       AliDebug(AliLog::kError, "Header not available");
 519       return;
 520     }
 521
 522     if (fUseMCVertex)
 523     {
 524       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 525       // get the MC vertex
 526       AliGenEventHeader* genHeader = header->GenEventHeader();
 527       TArrayF vtxMC(3);
 528       genHeader->PrimaryVertex(vtxMC);
 529
 530       vtx[2] = vtxMC[2];
 531     }
 532
 533     // get process information
 534     AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
 535
 536     Bool_t processEvent = kTRUE;
 537     if (fSelectProcessType > 0)
 538     {
 539       processEvent = kFALSE;
 540
 541       // non diffractive
 542       if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
 543         processEvent = kTRUE;
 544
 545       // single diffractive
 546       if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
 547         processEvent = kTRUE;
 548
 549       // double diffractive
 550       if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
 551         processEvent = kTRUE;
 552
 553       if (!processEvent)
 554         Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
 555     }
 556
 557     if (processEvent)
 558     {
 559       // get the MC vertex
 560       AliGenEventHeader* genHeader = header->GenEventHeader();
 561       TArrayF vtxMC(3);
 562       genHeader->PrimaryVertex(vtxMC);
 563
 564       // get number of tracks from MC
 565       // loop over mc particles
 566       Int_t nPrim  = stack->GetNprimary();
 567       Int_t nMCPart = stack->GetNtrack();
 568
 569       // tracks in different eta ranges
 570       Int_t nMCTracks05 = 0;
 571       Int_t nMCTracks10 = 0;
 572       Int_t nMCTracks14 = 0;
 573       Int_t nMCTracksAll = 0;
 574
 575       // tracks per particle species, in |eta| < 2 (systematic study)
 576       Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
 577       for (Int_t i = 0; i<4; ++i)
 578         nMCTracksSpecies[i] = 0;
 579
 580       for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 581       {
 582         AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 583
 584         TParticle* particle = stack->Particle(iMc);
 585
 586         if (!particle)
 587         {
 588           AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 589           continue;
 590         }
 591
 592         Bool_t debug = kFALSE;
 593         if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
 594         {
 595           //printf("%d) DROPPED ", iMC);
 596           //particle->Print();
 597           continue;
 598         }
 599
 600         //printf("%d) OK      ", iMC);
 601         //particle->Print();
 602
 603         //if (particle->Pt() < kPtCut)
 604         //  continue;
 605
 606         Int_t particleWeight = 1;
 607
 608         Float_t pt = particle->Pt();
 609
 610         // in case of systematic study, weight according to the change of the pt spectrum
 611         // it cannot be just multiplied because we cannot count "half of a particle"
 612         // instead a random generator decides if the particle is counted twice (if value > 1)
 613         // or not (if value < 0)
 614         if (fPtSpectrum)
 615         {
 616           Int_t bin = fPtSpectrum->FindBin(pt);
 617           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 618           {
 619             Float_t factor = fPtSpectrum->GetBinContent(bin);
 620             if (factor > 0)
 621             {
 622               Float_t random = gRandom->Uniform();
 623               if (factor > 1 && random < factor - 1)
 624               {
 625                 particleWeight = 2;
 626               }
 627               else if (factor < 1 && random < 1 - factor)
 628                 particleWeight = 0;
 629             }
 630           }
 631         }
 632
 633         //Printf("MC weight is: %d", particleWeight);
 634
 635         if (TMath::Abs(particle->Eta()) < 0.5)
 636           nMCTracks05 += particleWeight;
 637
 638         if (TMath::Abs(particle->Eta()) < 1.0)
 639           nMCTracks10 += particleWeight;
 640
 641         if (TMath::Abs(particle->Eta()) < 1.3)
 642           nMCTracks14 += particleWeight;
 643
 644         nMCTracksAll += particleWeight;
 645
 646         if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
 647           fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
 648
 649         if (fParticleCorrection[0] || fParticleSpecies)
 650         {
 651           Int_t id = -1;
 652           switch (TMath::Abs(particle->GetPdgCode()))
 653           {
 654             case 211: id = 0; break;
 655             case 321: id = 1; break;
 656             case 2212: id = 2; break;
 657             default: id = 3; break;
 658           }
 659
 660           if (TMath::Abs(particle->Eta()) < etaRange)
 661             nMCTracksSpecies[id]++;
 662
 663           if (fParticleCorrection[id])
 664             fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 665         }
 666       } // end of mc particle
 667
 668       fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
 669
 670       // ESD processing
 671       Int_t nESDTracks05 = 0;
 672       Int_t nESDTracks10 = 0;
 673       Int_t nESDTracks14 = 0;
 674
 675       // tracks per particle species, in |eta| < 2 (systematic study)
 676       Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
 677       for (Int_t i = 0; i<7; ++i)
 678         nESDTracksSpecies[i] = 0;
 679
 680       Bool_t* foundPrimaries = new Bool_t[nPrim];   // to prevent double counting
 681       for (Int_t i=0; i<nPrim; i++)
 682         foundPrimaries[i] = kFALSE;
 683
 684       Bool_t* foundPrimaries2 = new Bool_t[nPrim];   // to prevent double counting
 685       for (Int_t i=0; i<nPrim; i++)
 686         foundPrimaries2[i] = kFALSE;
 687
 688       Bool_t* foundTracks = new Bool_t[nMCPart];    // to prevent double counting
 689       for (Int_t i=0; i<nMCPart; i++)
 690         foundTracks[i] = kFALSE;
 691
 692       for (Int_t i=0; i<inputCount; ++i)
 693       {
 694         Float_t eta = etaArr[i];
 695         Int_t label = labelArr[i];
 696
 697         Int_t particleWeight = 1;
 698
 699         // in case of systematic study, weight according to the change of the pt spectrum
 700         if (fPtSpectrum)
 701         {
 702           TParticle* mother = 0;
 703
 704           // preserve label for later
 705           Int_t labelCopy = label;
 706           if (labelCopy >= 0)
 707             labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
 708           if (labelCopy >= 0)
 709             mother = stack->Particle(labelCopy);
 710
 711           // in case of pt study we do not count particles w/o label, because they cannot be scaled
 712           if (!mother)
 713             continue;
 714
 715           // it cannot be just multiplied because we cannot count "half of a particle"
 716           // instead a random generator decides if the particle is counted twice (if value > 1)
 717           // or not (if value < 0)
 718           Int_t bin = fPtSpectrum->FindBin(mother->Pt());
 719           if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
 720           {
 721             Float_t factor = fPtSpectrum->GetBinContent(bin);
 722             if (factor > 0)
 723             {
 724               Float_t random = gRandom->Uniform();
 725               if (factor > 1 && random < factor - 1)
 726               {
 727                 particleWeight = 2;
 728               }
 729               else if (factor < 1 && random < 1 - factor)
 730                 particleWeight = 0;
 731             }
 732           }
 733         }
 734
 735         //Printf("ESD weight is: %d", particleWeight);
 736
 737         if (TMath::Abs(eta) < 0.5)
 738           nESDTracks05 += particleWeight;
 739
 740         if (TMath::Abs(eta) < 1.0)
 741           nESDTracks10 += particleWeight;
 742
 743         if (TMath::Abs(eta) < 1.3)
 744           nESDTracks14 += particleWeight;
 745
 746         if (fParticleSpecies)
 747         {
 748           Int_t motherLabel = -1;
 749           TParticle* mother = 0;
 750
 751           // find mother
 752           if (label >= 0)
 753             motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
 754           if (motherLabel >= 0)
 755             mother = stack->Particle(motherLabel);
 756
 757           if (!mother)
 758           {
 759             // count tracks that did not have a label
 760             if (TMath::Abs(eta) < etaRange)
 761               nESDTracksSpecies[4]++;
 762           }
 763           else
 764           {
 765             // get particle type (pion, proton, kaon, other) of mother
 766             Int_t idMother = -1;
 767             switch (TMath::Abs(mother->GetPdgCode()))
 768             {
 769               case 211: idMother = 0; break;
 770               case 321: idMother = 1; break;
 771               case 2212: idMother = 2; break;
 772               default: idMother = 3; break;
 773             }
 774
 775             // double counting is ok for particle ratio study
 776             if (TMath::Abs(eta) < etaRange)
 777               nESDTracksSpecies[idMother]++;
 778
 779             // double counting is not ok for efficiency study
 780
 781             // check if we already counted this particle, this way distinguishes double counted particles (bug/artefact in tracking) or double counted primaries due to secondaries (physics)
 782             if (foundTracks[label])
 783             {
 784               if (TMath::Abs(eta) < etaRange)
 785                 nESDTracksSpecies[6]++;
 786             }
 787             else
 788             {
 789               foundTracks[label] = kTRUE;
 790
 791               // particle (primary) already counted?
 792               if (foundPrimaries[motherLabel])
 793               {
 794                 if (TMath::Abs(eta) < etaRange)
 795                   nESDTracksSpecies[5]++;
 796               }
 797               else
 798                 foundPrimaries[motherLabel] = kTRUE;
 799             }
 800           }
 801         }
 802
 803         if (fParticleCorrection[0])
 804         {
 805           if (label >= 0 && stack->IsPhysicalPrimary(label))
 806           {
 807             TParticle* particle = stack->Particle(label);
 808
 809             // get particle type (pion, proton, kaon, other)
 810             Int_t id = -1;
 811             switch (TMath::Abs(particle->GetPdgCode()))
 812             {
 813               case 211: id = 0; break;
 814               case 321: id = 1; break;
 815               case 2212: id = 2; break;
 816               default: id = 3; break;
 817             }
 818
 819             // todo check if values are not completely off??
 820
 821             // particle (primary) already counted?
 822             if (!foundPrimaries2[label])
 823             {
 824               foundPrimaries2[label] = kTRUE;
 825               fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 826             }
 827           }
 828         }
 829       }
 830
 831       if (fParticleCorrection[0])
 832       {
 833         // if the particle decays/stops before this radius we do not see it
 834         // 8cm larger than SPD layer 2
 835         // 123cm TPC radius where a track has about 50 clusters (cut limit)
 836         const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
 837
 838         // loop over all primaries that have not been found
 839         for (Int_t i=0; i<nPrim; i++)
 840         {
 841           // already found
 842           if (foundPrimaries2[i])
 843             continue;
 844
 845           TParticle* particle = 0;
 846           TClonesArray* trackrefs = 0;
 847           mcEvent->GetParticleAndTR(i, particle, trackrefs);
 848
 849           // true primary and charged
 850           if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
 851             continue;
 852
 853           //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
 854           if (TMath::Abs(particle->Eta()) > 3)
 855             continue;
 856
 857           // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
 858
 859           // get particle type (pion, proton, kaon, other)
 860           Int_t id = -1;
 861           switch (TMath::Abs(particle->GetPdgCode()))
 862           {
 863             case 211: id = 4; break;
 864             case 321: id = 5; break;
 865             case 2212: id = 6; break;
 866             default: id = 7; break;
 867           }
 868
 869           if (!fParticleCorrection[id])
 870             continue;
 871
 872           // get last track reference
 873           AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
 874
 875           if (!trackref)
 876           {
 877             Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
 878             particle->Print();
 879             continue;
 880           }
 881
 882           // particle in tracking volume long enough...
 883           if (trackref->R() > endRadius)
 884             continue;
 885
 886           if (particle->GetLastDaughter() >= 0)
 887           {
 888             Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
 889             //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
 890             if (uID == kPDecay)
 891             {
 892               // decayed
 893
 894               Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
 895               particle->Print();
 896               Printf("Daugthers:");
 897               for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
 898                 stack->Particle(d)->Print();
 899               Printf(" ");
 900
 901               fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
 902               continue;
 903             }
 904           }
 905
 906           if (trackref->DetectorId() == -1)
 907           {
 908             // stopped
 909             Printf("Particle %d stopped at %f:", i, trackref->R());
 910             particle->Print();
 911             Printf(" ");
 912
 913             fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
 914             continue;
 915           }
 916
 917           Printf("Particle %d simply not tracked", i);
 918           particle->Print();
 919           Printf(" ");
 920         }
 921       }
 922
 923       delete[] foundTracks;
 924       delete[] foundPrimaries;
 925       delete[] foundPrimaries2;
 926
 927 //         if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
 928 //         {
 929 //             TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 930 //             printf("WARNING: Event %lld %s (vtx-z = %f, recon: %f, contrib: %d, res: %f) has %d generated and %d reconstructed...\n", tree->GetReadEntry(), tree->GetCurrentFile()->GetName(), vtxMC[2], vtx[2], vtxESD->GetNContributors(), vtxESD->GetZRes(), nMCTracks14, nESDTracks14);
 931 //         }
 932
 933       if (eventTriggered)
 934       {
 935         fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
 936         fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
 937 //         if (!eventVertex)
 938 //         {
 939 //           if (nESDTracks05 == 0)
 940 //             fTemp1->Fill(nMCTracks05, nESDTracks05);
 941 //         }
 942       }
 943
 944       if (eventTriggered && eventVertex)
 945       {
 946         // fill response matrix using vtxMC (best guess)
 947         fMultiplicity->FillCorrection(vtxMC[2],  nMCTracks05,  nMCTracks10,  nMCTracks14,  nMCTracksAll,  nESDTracks05, nESDTracks10, nESDTracks14);
 948
 949         fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
 950
 951         if (fParticleSpecies)
 952           fParticleSpecies->Fill(vtxMC[2], nMCTracksSpecies[0], nMCTracksSpecies[1], nMCTracksSpecies[2], nMCTracksSpecies[3], nESDTracksSpecies[0], nESDTracksSpecies[1], nESDTracksSpecies[2], nESDTracksSpecies[3], nESDTracksSpecies[4], nESDTracksSpecies[5], nESDTracksSpecies[6]);
 953       }
 954     }
 955   }
 956
 957   if (etaArr)
 958     delete[] etaArr;
 959   if (labelArr)
 960     delete[] labelArr;
 961
 962 }
 963
 964 void AliMultiplicityTask::Terminate(Option_t *)
 965 {
 966   // The Terminate() function is the last function to be called during
 967   // a query. It always runs on the client, it can be used to present
 968   // the results graphically or save the results to file.
 969
 970   fOutput = dynamic_cast<TList*> (GetOutputData(0));
 971   if (!fOutput) {
 972     Printf("ERROR: fOutput not available");
 973     return;
 974   }
 975
 976   fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
 977   for (Int_t i = 0; i < 8; ++i)
 978     fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
 979   fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
 980
 981   fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
 982
 983   if (!fMultiplicity)
 984   {
 985     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
 986     return;
 987   }
 988
 989   TString fileName("multiplicity");
 990   if (fSelectProcessType == 1)
 991     fileName += "ND";
 992   if (fSelectProcessType == 2)
 993     fileName += "SD";
 994   if (fSelectProcessType == 3)
 995     fileName += "DD";
 996   fileName += ".root";
 997   TFile* file = TFile::Open(fileName, "RECREATE");
 998
 999   fMultiplicity->SaveHistograms();
1000   for (Int_t i = 0; i < 8; ++i)
1001     if (fParticleCorrection[i])
1002       fParticleCorrection[i]->SaveHistograms();
1003   if (fParticleSpecies)
1004     fParticleSpecies->Write();
1005   if (fdNdpT)
1006     fdNdpT->Write();
1007
1008   fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1009   if (fTemp1)
1010     fTemp1->Write();
1011
1012   fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1013   if (fEtaPhi)
1014     fEtaPhi->Write();
1015
1016   fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1017   if (fVertex)
1018     fVertex->Write();
1019
1020   for (Int_t i=0; i<3; i++)
1021   {
1022     fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
1023     if (fEta[i])
1024     {
1025       fEta[i]->Sumw2();
1026       Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
1027       if (events > 0)
1028         fEta[i]->Scale(1.0 / events);
1029       fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
1030       fEta[i]->Write();
1031     }
1032   }
1033
1034   TObjString option(fOption);
1035   option.Write();
1036
1037   file->Close();
1038
1039   Printf("Written result to %s", fileName.Data());
1040 }