Fix for CMake
[u/mrichter/AliRoot.git] / PWG0 / multiplicity / AliMultiplicityTask.cxx
CommitLineData
a9017e49 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>
69b09e3b 11#include <TH1D.h>
a9017e49 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
745d6088 32#include "AliESDtrackCuts.h"
a9017e49 33#include "AliPWG0Helper.h"
0a0f4adb 34#include "multiplicity/AliMultiplicityCorrection.h"
a9017e49 35#include "AliCorrection.h"
36#include "AliCorrectionMatrix3D.h"
37
a9017e49 38ClassImp(AliMultiplicityTask)
39
40AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
41 AliAnalysisTask("AliMultiplicityTask", ""),
42 fESD(0),
43 fOption(opt),
745d6088 44 fAnalysisMode(AliPWG0Helper::kSPD),
a2cdb42c 45 fTrigger(AliPWG0Helper::kMB1),
69b09e3b 46 fDeltaPhiCut(-1),
2fa65f52 47 fReadMC(kFALSE),
f3eb27f6 48 fUseMCVertex(kFALSE),
a9017e49 49 fMultiplicity(0),
50 fEsdTrackCuts(0),
51 fSystSkipParticles(kFALSE),
52 fSelectProcessType(0),
53 fParticleSpecies(0),
69b09e3b 54 fdNdpT(0),
a9017e49 55 fPtSpectrum(0),
56 fOutput(0)
57{
58 //
59 // Constructor. Initialization of pointers
60 //
61
69b09e3b 62 for (Int_t i = 0; i<8; i++)
a9017e49 63 fParticleCorrection[i] = 0;
64
65 // Define input and output slots here
66 DefineInput(0, TChain::Class());
67 DefineOutput(0, TList::Class());
68}
69
70AliMultiplicityTask::~AliMultiplicityTask()
71{
72 //
73 // Destructor
74 //
75
76 // histograms are in the output list and deleted when the output
77 // list is deleted by the TSelector dtor
78
79 if (fOutput) {
80 delete fOutput;
81 fOutput = 0;
82 }
83}
84
85//________________________________________________________________________
86void AliMultiplicityTask::ConnectInputData(Option_t *)
87{
88 // Connect ESD
89 // Called once
90
2fa65f52 91 Printf("AliMultiplicityTask::ConnectInputData called");
92
a9017e49 93 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
94 if (!tree) {
95 Printf("ERROR: Could not read tree from input slot 0");
96 } else {
97 // Disable all branches and enable only the needed ones
69b09e3b 98 /*
f3eb27f6 99 tree->SetBranchStatus("*", 0);
2fa65f52 100
f3eb27f6 101 tree->SetBranchStatus("AliESDHeader*", 1);
102 tree->SetBranchStatus("*Vertex*", 1);
a9017e49 103
f3eb27f6 104 if (fAnalysisMode == AliPWG0Helper::kSPD) {
105 tree->SetBranchStatus("AliMultiplicity*", 1);
106 }
a9017e49 107
745d6088 108 if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
f3eb27f6 109 //AliESDtrackCuts::EnableNeededBranches(tree);
110 tree->SetBranchStatus("Tracks*", 1);
2fa65f52 111 }
69b09e3b 112 */
a9017e49 113
a9017e49 114 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
115
116 if (!esdH) {
117 Printf("ERROR: Could not get ESDInputHandler");
118 } else
119 fESD = esdH->GetEvent();
120 }
f3eb27f6 121
122 // disable info messages of AliMCEvent (per event)
123 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
a9017e49 124}
125
126void AliMultiplicityTask::CreateOutputObjects()
127{
128 // create result objects and add to output list
129
130 fOutput = new TList;
131 fOutput->SetOwner();
132
133 fMultiplicity = new AliMultiplicityCorrection("Multiplicity", "Multiplicity");
134 fOutput->Add(fMultiplicity);
69b09e3b 135
136 fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
137 fdNdpT->Sumw2();
138 fOutput->Add(fdNdpT);
a9017e49 139
140 if (fOption.Contains("skip-particles"))
141 {
142 fSystSkipParticles = kTRUE;
143 AliInfo("WARNING: Systematic study enabled. Particles will be skipped.");
144 }
145
146 if (fOption.Contains("particle-efficiency"))
69b09e3b 147 for (Int_t i = 0; i<8; i++)
a9017e49 148 {
149 fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
150 fOutput->Add(fParticleCorrection[i]);
151 }
152
153 if (fOption.Contains("only-process-type-nd"))
154 fSelectProcessType = 1;
155
156 if (fOption.Contains("only-process-type-sd"))
157 fSelectProcessType = 2;
158
159 if (fOption.Contains("only-process-type-dd"))
160 fSelectProcessType = 3;
161
162 if (fSelectProcessType != 0)
163 AliInfo(Form("WARNING: Systematic study enabled. Only considering process type %d", fSelectProcessType));
164
165 if (fOption.Contains("pt-spectrum-hist"))
166 {
167 TFile* file = TFile::Open("ptspectrum_fit.root");
168 if (file)
169 {
170 TString subStr(fOption(fOption.Index("pt-spectrum")+17, 3));
171 TString histName(Form("ptspectrum_%s", subStr.Data()));
172 AliInfo(Form("Pt-Spectrum modification. Using %s.", histName.Data()));
69b09e3b 173 fPtSpectrum = (TH1D*) file->Get(histName);
174 if (!fPtSpectrum)
175 AliError("Histogram not found");
a9017e49 176 }
177 else
178 AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
a9017e49 179 }
180
181 if (fOption.Contains("pt-spectrum-func"))
182 {
183 if (fPtSpectrum)
184 {
69b09e3b 185 Printf("Using function for systematic p_t study");
a9017e49 186 }
187 else
188 {
69b09e3b 189 Printf("ERROR: Could not find function for systematic p_t study");
190 fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
a9017e49 191 fPtSpectrum->SetBinContent(1, 1);
192 }
a9017e49 193 }
194
69b09e3b 195 if (fPtSpectrum)
196 Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
197
a9017e49 198 if (fOption.Contains("particle-species")) {
199 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");
200 fOutput->Add(fParticleSpecies);
201 }
202
203 // TODO set seed for random generator
204}
205
206void AliMultiplicityTask::Exec(Option_t*)
207{
208 // process the event
209
210 // Check prerequisites
211 if (!fESD)
212 {
745d6088 213 AliDebug(AliLog::kError, "ESD not available");
a9017e49 214 return;
215 }
216
69b09e3b 217 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD, fTrigger);
7dcf959e 218 //Printf("%lld", fESD->GetTriggerMask());
a9017e49 219
f3eb27f6 220 const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
221 Bool_t eventVertex = (vtxESD != 0);
a9017e49 222
a9017e49 223 Double_t vtx[3];
f3eb27f6 224 if (vtxESD)
225 vtxESD->GetXYZ(vtx);
a9017e49 226
2fa65f52 227 // post the data already here
228 PostData(0, fOutput);
69b09e3b 229
a9017e49 230 //const Float_t kPtCut = 0.3;
231
2fa65f52 232 // create list of (label, eta) tuples
233 Int_t inputCount = 0;
234 Int_t* labelArr = 0;
235 Float_t* etaArr = 0;
745d6088 236 if (fAnalysisMode == AliPWG0Helper::kSPD)
a9017e49 237 {
2fa65f52 238 // get tracklets
239 const AliMultiplicity* mult = fESD->GetMultiplicity();
240 if (!mult)
a9017e49 241 {
2fa65f52 242 AliDebug(AliLog::kError, "AliMultiplicity not available");
243 return;
a9017e49 244 }
245
2fa65f52 246 labelArr = new Int_t[mult->GetNumberOfTracklets()];
247 etaArr = new Float_t[mult->GetNumberOfTracklets()];
248
249 // get multiplicity from ITS tracklets
250 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
a9017e49 251 {
2fa65f52 252 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
a9017e49 253
69b09e3b 254 Float_t deltaPhi = mult->GetDeltaPhi(i);
255
256 if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
2fa65f52 257 continue;
a9017e49 258
2fa65f52 259 etaArr[inputCount] = mult->GetEta(i);
69b09e3b 260 if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
261 {
262 labelArr[inputCount] = mult->GetLabel(i, 0);
263 }
264 else
265 labelArr[inputCount] = -1;
266
2fa65f52 267 ++inputCount;
268 }
269 }
745d6088 270 else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
2fa65f52 271 {
272 if (!fEsdTrackCuts)
273 {
274 AliDebug(AliLog::kError, "fESDTrackCuts not available");
275 return;
276 }
a9017e49 277
69b09e3b 278 if (vtxESD)
a9017e49 279 {
69b09e3b 280 // get multiplicity from ESD tracks
281 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, (fAnalysisMode == AliPWG0Helper::kTPC));
282 Int_t nGoodTracks = list->GetEntries();
283
284 labelArr = new Int_t[nGoodTracks];
285 etaArr = new Float_t[nGoodTracks];
286
287 // loop over esd tracks
288 for (Int_t i=0; i<nGoodTracks; i++)
a9017e49 289 {
69b09e3b 290 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
291 if (!esdTrack)
292 {
293 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
294 continue;
295 }
296
297 etaArr[inputCount] = esdTrack->Eta();
298 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
299 ++inputCount;
a9017e49 300 }
69b09e3b 301
302 delete list;
a9017e49 303 }
2fa65f52 304 }
a9017e49 305
2fa65f52 306 // eta range for nMCTracksSpecies, nESDTracksSpecies
307 Float_t etaRange = -1;
308 switch (fAnalysisMode) {
773f84dc 309 case AliPWG0Helper::kInvalid: break;
745d6088 310 case AliPWG0Helper::kTPC:
311 case AliPWG0Helper::kTPCITS:
69b09e3b 312 etaRange = 1.0; break;
313 case AliPWG0Helper::kSPD: etaRange = 1.0; break;
eedcd663 314 default: break;
2fa65f52 315 }
a9017e49 316
2fa65f52 317 if (!fReadMC) // Processing of ESD information
318 {
319 if (eventTriggered && eventVertex)
320 {
321 Int_t nESDTracks05 = 0;
69b09e3b 322 Int_t nESDTracks10 = 0;
b3b260d1 323 Int_t nESDTracks14 = 0;
a9017e49 324
2fa65f52 325 for (Int_t i=0; i<inputCount; ++i)
326 {
327 Float_t eta = etaArr[i];
a9017e49 328
2fa65f52 329 if (TMath::Abs(eta) < 0.5)
330 nESDTracks05++;
a9017e49 331
69b09e3b 332 if (TMath::Abs(eta) < 1.0)
333 nESDTracks10++;
a9017e49 334
b3b260d1 335 if (TMath::Abs(eta) < 1.4)
336 nESDTracks14++;
a9017e49 337 }
2fa65f52 338
b3b260d1 339 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 340 }
341 }
342 else if (fReadMC) // Processing of MC information
343 {
344 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
345 if (!eventHandler) {
346 Printf("ERROR: Could not retrieve MC event handler");
347 return;
a9017e49 348 }
a9017e49 349
2fa65f52 350 AliMCEvent* mcEvent = eventHandler->MCEvent();
351 if (!mcEvent) {
352 Printf("ERROR: Could not retrieve MC event");
353 return;
354 }
a9017e49 355
2fa65f52 356 AliStack* stack = mcEvent->Stack();
357 if (!stack)
358 {
359 AliDebug(AliLog::kError, "Stack not available");
360 return;
361 }
69b09e3b 362
2fa65f52 363 AliHeader* header = mcEvent->Header();
364 if (!header)
365 {
366 AliDebug(AliLog::kError, "Header not available");
367 return;
368 }
a9017e49 369
f3eb27f6 370 if (fUseMCVertex)
371 {
372 Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
373 // get the MC vertex
374 AliGenEventHeader* genHeader = header->GenEventHeader();
375 TArrayF vtxMC(3);
376 genHeader->PrimaryVertex(vtxMC);
377
378 vtx[2] = vtxMC[2];
379 }
69b09e3b 380
381 // get process information
382 AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(header);
f3eb27f6 383
2fa65f52 384 Bool_t processEvent = kTRUE;
385 if (fSelectProcessType > 0)
386 {
2fa65f52 387 processEvent = kFALSE;
a9017e49 388
2fa65f52 389 // non diffractive
69b09e3b 390 if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
2fa65f52 391 processEvent = kTRUE;
a9017e49 392
2fa65f52 393 // single diffractive
69b09e3b 394 if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
2fa65f52 395 processEvent = kTRUE;
a9017e49 396
2fa65f52 397 // double diffractive
69b09e3b 398 if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
2fa65f52 399 processEvent = kTRUE;
a9017e49 400
2fa65f52 401 if (!processEvent)
69b09e3b 402 Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
2fa65f52 403 }
a9017e49 404
2fa65f52 405 if (processEvent)
406 {
407 // get the MC vertex
408 AliGenEventHeader* genHeader = header->GenEventHeader();
409 TArrayF vtxMC(3);
410 genHeader->PrimaryVertex(vtxMC);
411
412 // get number of tracks from MC
413 // loop over mc particles
414 Int_t nPrim = stack->GetNprimary();
415 Int_t nMCPart = stack->GetNtrack();
416
417 // tracks in different eta ranges
418 Int_t nMCTracks05 = 0;
69b09e3b 419 Int_t nMCTracks10 = 0;
b3b260d1 420 Int_t nMCTracks14 = 0;
2fa65f52 421 Int_t nMCTracksAll = 0;
422
423 // tracks per particle species, in |eta| < 2 (systematic study)
424 Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
425 for (Int_t i = 0; i<4; ++i)
426 nMCTracksSpecies[i] = 0;
427
428 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
429 {
430 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
a9017e49 431
2fa65f52 432 TParticle* particle = stack->Particle(iMc);
a9017e49 433
2fa65f52 434 if (!particle)
435 {
436 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
437 continue;
438 }
a9017e49 439
2fa65f52 440 Bool_t debug = kFALSE;
441 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
442 {
443 //printf("%d) DROPPED ", iMC);
444 //particle->Print();
445 continue;
446 }
a9017e49 447
2fa65f52 448 //printf("%d) OK ", iMC);
449 //particle->Print();
a9017e49 450
2fa65f52 451 //if (particle->Pt() < kPtCut)
452 // continue;
a9017e49 453
2fa65f52 454 Int_t particleWeight = 1;
a9017e49 455
2fa65f52 456 Float_t pt = particle->Pt();
a9017e49 457
2fa65f52 458 // in case of systematic study, weight according to the change of the pt spectrum
459 // it cannot be just multiplied because we cannot count "half of a particle"
460 // instead a random generator decides if the particle is counted twice (if value > 1)
461 // or not (if value < 0)
462 if (fPtSpectrum)
a9017e49 463 {
2fa65f52 464 Int_t bin = fPtSpectrum->FindBin(pt);
465 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
a9017e49 466 {
2fa65f52 467 Float_t factor = fPtSpectrum->GetBinContent(bin);
468 if (factor > 0)
469 {
470 Float_t random = gRandom->Uniform();
471 if (factor > 1 && random < factor - 1)
472 {
473 particleWeight = 2;
474 }
475 else if (factor < 1 && random < 1 - factor)
476 particleWeight = 0;
477 }
a9017e49 478 }
a9017e49 479 }
a9017e49 480
2fa65f52 481 //Printf("MC weight is: %d", particleWeight);
a9017e49 482
2fa65f52 483 if (TMath::Abs(particle->Eta()) < 0.5)
484 nMCTracks05 += particleWeight;
a9017e49 485
69b09e3b 486 if (TMath::Abs(particle->Eta()) < 1.0)
487 nMCTracks10 += particleWeight;
a9017e49 488
b3b260d1 489 if (TMath::Abs(particle->Eta()) < 1.4)
490 nMCTracks14 += particleWeight;
a9017e49 491
2fa65f52 492 nMCTracksAll += particleWeight;
69b09e3b 493
494 if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
495 fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
a9017e49 496
2fa65f52 497 if (fParticleCorrection[0] || fParticleSpecies)
498 {
499 Int_t id = -1;
500 switch (TMath::Abs(particle->GetPdgCode()))
501 {
502 case 211: id = 0; break;
503 case 321: id = 1; break;
504 case 2212: id = 2; break;
505 default: id = 3; break;
506 }
a9017e49 507
2fa65f52 508 if (TMath::Abs(particle->Eta()) < etaRange)
509 nMCTracksSpecies[id]++;
510
511 if (fParticleCorrection[id])
512 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
513 }
514 } // end of mc particle
a9017e49 515
b3b260d1 516 fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
a9017e49 517
2fa65f52 518 if (eventTriggered && eventVertex)
a9017e49 519 {
2fa65f52 520 Int_t nESDTracks05 = 0;
69b09e3b 521 Int_t nESDTracks10 = 0;
b3b260d1 522 Int_t nESDTracks14 = 0;
a9017e49 523
2fa65f52 524 // tracks per particle species, in |eta| < 2 (systematic study)
525 Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
526 for (Int_t i = 0; i<7; ++i)
527 nESDTracksSpecies[i] = 0;
a9017e49 528
2fa65f52 529 Bool_t* foundPrimaries = new Bool_t[nPrim]; // to prevent double counting
530 for (Int_t i=0; i<nPrim; i++)
531 foundPrimaries[i] = kFALSE;
a9017e49 532
69b09e3b 533 Bool_t* foundPrimaries2 = new Bool_t[nPrim]; // to prevent double counting
534 for (Int_t i=0; i<nPrim; i++)
535 foundPrimaries2[i] = kFALSE;
536
2fa65f52 537 Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
538 for (Int_t i=0; i<nMCPart; i++)
539 foundTracks[i] = kFALSE;
a9017e49 540
2fa65f52 541 for (Int_t i=0; i<inputCount; ++i)
542 {
543 Float_t eta = etaArr[i];
544 Int_t label = labelArr[i];
a9017e49 545
2fa65f52 546 Int_t particleWeight = 1;
547
69b09e3b 548 // systematic study: 5% lower efficiency
549 if (fSystSkipParticles && (gRandom->Uniform() < 0.05))
550 continue;
551
2fa65f52 552 // in case of systematic study, weight according to the change of the pt spectrum
553 if (fPtSpectrum)
554 {
555 TParticle* mother = 0;
556
557 // preserve label for later
558 Int_t labelCopy = label;
559 if (labelCopy >= 0)
560 labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
561 if (labelCopy >= 0)
562 mother = stack->Particle(labelCopy);
563
564 // in case of pt study we do not count particles w/o label, because they cannot be scaled
565 if (!mother)
566 continue;
567
568 // it cannot be just multiplied because we cannot count "half of a particle"
569 // instead a random generator decides if the particle is counted twice (if value > 1)
570 // or not (if value < 0)
571 Int_t bin = fPtSpectrum->FindBin(mother->Pt());
572 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
573 {
574 Float_t factor = fPtSpectrum->GetBinContent(bin);
575 if (factor > 0)
576 {
577 Float_t random = gRandom->Uniform();
578 if (factor > 1 && random < factor - 1)
579 {
580 particleWeight = 2;
581 }
582 else if (factor < 1 && random < 1 - factor)
583 particleWeight = 0;
584 }
585 }
586 }
587
588 //Printf("ESD weight is: %d", particleWeight);
589
590 if (TMath::Abs(eta) < 0.5)
591 nESDTracks05 += particleWeight;
a9017e49 592
69b09e3b 593 if (TMath::Abs(eta) < 1.0)
594 nESDTracks10 += particleWeight;
a9017e49 595
b3b260d1 596 if (TMath::Abs(eta) < 1.4)
597 nESDTracks14 += particleWeight;
a9017e49 598
69b09e3b 599 if (fParticleSpecies)
2fa65f52 600 {
601 Int_t motherLabel = -1;
602 TParticle* mother = 0;
603
604 // find mother
605 if (label >= 0)
606 motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
607 if (motherLabel >= 0)
608 mother = stack->Particle(motherLabel);
609
610 if (!mother)
611 {
612 // count tracks that did not have a label
613 if (TMath::Abs(eta) < etaRange)
614 nESDTracksSpecies[4]++;
2fa65f52 615 }
69b09e3b 616 else
2fa65f52 617 {
69b09e3b 618 // get particle type (pion, proton, kaon, other) of mother
619 Int_t idMother = -1;
620 switch (TMath::Abs(mother->GetPdgCode()))
621 {
622 case 211: idMother = 0; break;
623 case 321: idMother = 1; break;
624 case 2212: idMother = 2; break;
625 default: idMother = 3; break;
626 }
627
628 // double counting is ok for particle ratio study
629 if (TMath::Abs(eta) < etaRange)
630 nESDTracksSpecies[idMother]++;
631
632 // double counting is not ok for efficiency study
633
634 // 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)
635 if (foundTracks[label])
636 {
637 if (TMath::Abs(eta) < etaRange)
638 nESDTracksSpecies[6]++;
639 }
640 else
641 {
642 foundTracks[label] = kTRUE;
643
644 // particle (primary) already counted?
645 if (foundPrimaries[motherLabel])
646 {
647 if (TMath::Abs(eta) < etaRange)
648 nESDTracksSpecies[5]++;
649 }
650 else
651 foundPrimaries[motherLabel] = kTRUE;
652 }
2fa65f52 653 }
69b09e3b 654 }
655
656 if (fParticleCorrection[0])
657 {
658 if (label >= 0 && stack->IsPhysicalPrimary(label))
659 {
660 TParticle* particle = stack->Particle(label);
2fa65f52 661
69b09e3b 662 // get particle type (pion, proton, kaon, other)
663 Int_t id = -1;
664 switch (TMath::Abs(particle->GetPdgCode()))
665 {
666 case 211: id = 0; break;
667 case 321: id = 1; break;
668 case 2212: id = 2; break;
669 default: id = 3; break;
670 }
2fa65f52 671
69b09e3b 672 // todo check if values are not completely off??
2fa65f52 673
69b09e3b 674 // particle (primary) already counted?
675 if (!foundPrimaries2[label])
676 {
677 foundPrimaries2[label] = kTRUE;
678 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
679 }
680 }
681 }
682 }
683
684 if (fParticleCorrection[0])
685 {
686 // if the particle decays/stops before this radius we do not see it
687 // 8cm larger than SPD layer 2
688 // 123cm TPC radius where a track has about 50 clusters (cut limit)
689 const Float_t endRadius = (fAnalysisMode == AliPWG0Helper::kSPD) ? 8. : 123;
690
691 // loop over all primaries that have not been found
692 for (Int_t i=0; i<nPrim; i++)
693 {
694 // already found
695 if (foundPrimaries2[i])
696 continue;
697
698 TParticle* particle = 0;
699 TClonesArray* trackrefs = 0;
700 mcEvent->GetParticleAndTR(i, particle, trackrefs);
701
702 // true primary and charged
703 if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
704 continue;
705
706 //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
707 if (TMath::Abs(particle->Eta()) > 3)
708 continue;
709
710 // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
711
712 // get particle type (pion, proton, kaon, other)
713 Int_t id = -1;
714 switch (TMath::Abs(particle->GetPdgCode()))
2fa65f52 715 {
69b09e3b 716 case 211: id = 4; break;
717 case 321: id = 5; break;
718 case 2212: id = 6; break;
719 default: id = 7; break;
720 }
721
722 if (!fParticleCorrection[id])
723 continue;
724
725 // get last track reference
726 AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
727
728 if (!trackref)
729 {
730 Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
731 particle->Print();
2fa65f52 732 continue;
733 }
69b09e3b 734
735 // particle in tracking volume long enough...
736 if (trackref->R() > endRadius)
737 continue;
738
739 if (particle->GetLastDaughter() >= 0)
2fa65f52 740 {
69b09e3b 741 Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
742 //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
743 if (uID == kPDecay)
744 {
745 // decayed
746
747 Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
748 particle->Print();
749 Printf("Daughers:");
750 for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
751 stack->Particle(d)->Print();
752 Printf("");
753
754 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
755 continue;
756 }
757 }
758
759 if (trackref->DetectorId() == -1)
760 {
761 // stopped
762 Printf("Particle %d stopped at %f:", i, trackref->R());
763 particle->Print();
764 Printf("");
765
766 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
2fa65f52 767 continue;
768 }
69b09e3b 769
770 Printf("Particle %d simply not tracked", i);
771 particle->Print();
772 Printf("");
2fa65f52 773 }
774 }
69b09e3b 775
2fa65f52 776 delete[] foundTracks;
777 delete[] foundPrimaries;
69b09e3b 778 delete[] foundPrimaries2;
a9017e49 779
b3b260d1 780 if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
f3eb27f6 781 {
782 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
b3b260d1 783 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);
f3eb27f6 784 }
a9017e49 785
2fa65f52 786 // fill response matrix using vtxMC (best guess)
b3b260d1 787 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks14, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks14);
a9017e49 788
b3b260d1 789 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 790
791 if (fParticleSpecies)
792 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]);
793 }
794 }
795 }
a9017e49 796
69b09e3b 797 if (etaArr)
798 delete[] etaArr;
799 if (labelArr)
800 delete[] labelArr;
a9017e49 801}
802
803void AliMultiplicityTask::Terminate(Option_t *)
804{
805 // The Terminate() function is the last function to be called during
806 // a query. It always runs on the client, it can be used to present
807 // the results graphically or save the results to file.
808
809 fOutput = dynamic_cast<TList*> (GetOutputData(0));
810 if (!fOutput) {
811 Printf("ERROR: fOutput not available");
812 return;
813 }
814
815 fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
69b09e3b 816 for (Int_t i = 0; i < 8; ++i)
a9017e49 817 fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
818 fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
69b09e3b 819
820 fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
a9017e49 821
822 if (!fMultiplicity)
823 {
824 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
825 return;
826 }
827
dca331bb 828 TFile* file = TFile::Open("multiplicity.root", "RECREATE");
a9017e49 829
830 fMultiplicity->SaveHistograms();
69b09e3b 831 for (Int_t i = 0; i < 8; ++i)
a9017e49 832 if (fParticleCorrection[i])
833 fParticleCorrection[i]->SaveHistograms();
834 if (fParticleSpecies)
835 fParticleSpecies->Write();
69b09e3b 836 if (fdNdpT)
837 fdNdpT->Write();
a9017e49 838
839 TObjString option(fOption);
840 option.Write();
841
842 file->Close();
843
7dcf959e 844 Printf("Written result to multiplicity.root");
a9017e49 845}