]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG0/multiplicity/AliMultiplicityTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Implemented new multiplicity estimator (tracks+tracklets-made-of-clusters-not-used...
[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"
eb9356d5 37#include "AliPhysicsSelection.h"
70fdd197 38#include "AliTriggerAnalysis.h"
2701e5ae 39#include "AliVEvent.h"
a9017e49 40
a9017e49 41ClassImp(AliMultiplicityTask)
42
43AliMultiplicityTask::AliMultiplicityTask(const char* opt) :
44 AliAnalysisTask("AliMultiplicityTask", ""),
45 fESD(0),
46 fOption(opt),
70fdd197 47 fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kSPD | AliPWG0Helper::kFieldOn)),
48 fTrigger(AliTriggerAnalysis::kMB1),
69b09e3b 49 fDeltaPhiCut(-1),
eb9356d5 50 fDiffTreatment(AliPWG0Helper::kMCFlags),
2fa65f52 51 fReadMC(kFALSE),
f3eb27f6 52 fUseMCVertex(kFALSE),
a9017e49 53 fMultiplicity(0),
54 fEsdTrackCuts(0),
55 fSystSkipParticles(kFALSE),
56 fSelectProcessType(0),
57 fParticleSpecies(0),
69b09e3b 58 fdNdpT(0),
a9017e49 59 fPtSpectrum(0),
eb9356d5 60 fTemp1(0),
61 fTemp2(0),
f55b30d2 62 fVertex(0),
63 fEtaPhi(0),
a9017e49 64 fOutput(0)
65{
66 //
67 // Constructor. Initialization of pointers
68 //
69
69b09e3b 70 for (Int_t i = 0; i<8; i++)
a9017e49 71 fParticleCorrection[i] = 0;
eb9356d5 72
73 for (Int_t i=0; i<3; i++)
74 fEta[i] = 0;
a9017e49 75
76 // Define input and output slots here
77 DefineInput(0, TChain::Class());
78 DefineOutput(0, TList::Class());
eb9356d5 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));
a9017e49 91}
92
93AliMultiplicityTask::~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
2701e5ae 102 if (fOutput&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
a9017e49 103 delete fOutput;
104 fOutput = 0;
105 }
106}
107
108//________________________________________________________________________
109void AliMultiplicityTask::ConnectInputData(Option_t *)
110{
111 // Connect ESD
112 // Called once
113
2fa65f52 114 Printf("AliMultiplicityTask::ConnectInputData called");
115
a9017e49 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
69b09e3b 121 /*
f3eb27f6 122 tree->SetBranchStatus("*", 0);
2fa65f52 123
f3eb27f6 124 tree->SetBranchStatus("AliESDHeader*", 1);
125 tree->SetBranchStatus("*Vertex*", 1);
a9017e49 126
70fdd197 127 if (fAnalysisMode & AliPWG0Helper::kSPD) {
f3eb27f6 128 tree->SetBranchStatus("AliMultiplicity*", 1);
129 }
a9017e49 130
70fdd197 131 if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS) {
f3eb27f6 132 //AliESDtrackCuts::EnableNeededBranches(tree);
133 tree->SetBranchStatus("Tracks*", 1);
2fa65f52 134 }
69b09e3b 135 */
a9017e49 136
a9017e49 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 }
f3eb27f6 144
145 // disable info messages of AliMCEvent (per event)
146 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
a9017e49 147}
148
149void 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);
69b09e3b 158
159 fdNdpT = new TH1F("fdNdpT", "fdNdpT", 1000, 0, 10);
160 fdNdpT->Sumw2();
161 fOutput->Add(fdNdpT);
a9017e49 162
a9017e49 163 if (fOption.Contains("particle-efficiency"))
69b09e3b 164 for (Int_t i = 0; i<8; i++)
a9017e49 165 {
166 fParticleCorrection[i] = new AliCorrection(Form("correction_%d", i), Form("correction_%d", i));
167 fOutput->Add(fParticleCorrection[i]);
168 }
169
a9017e49 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()));
69b09e3b 178 fPtSpectrum = (TH1D*) file->Get(histName);
179 if (!fPtSpectrum)
180 AliError("Histogram not found");
a9017e49 181 }
182 else
183 AliError("Could not open ptspectrum_fit.root. Pt Spectrum study could not be enabled.");
a9017e49 184 }
185
186 if (fOption.Contains("pt-spectrum-func"))
187 {
188 if (fPtSpectrum)
189 {
69b09e3b 190 Printf("Using function for systematic p_t study");
a9017e49 191 }
192 else
193 {
69b09e3b 194 Printf("ERROR: Could not find function for systematic p_t study");
195 fPtSpectrum = new TH1D("fPtSpectrum", "fPtSpectrum", 1, 0, 100);
a9017e49 196 fPtSpectrum->SetBinContent(1, 1);
197 }
a9017e49 198 }
199
69b09e3b 200 if (fPtSpectrum)
201 Printf("WARNING: Systematic study enabled. Pt spectrum will be modified");
202
a9017e49 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
eb9356d5 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
f55b30d2 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
a9017e49 223 // TODO set seed for random generator
224}
225
226void AliMultiplicityTask::Exec(Option_t*)
227{
228 // process the event
229
230 // Check prerequisites
231 if (!fESD)
232 {
745d6088 233 AliDebug(AliLog::kError, "ESD not available");
a9017e49 234 return;
235 }
236
eb9356d5 237 AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
238 if (!inputHandler)
239 {
240 Printf("ERROR: Could not receive input handler");
241 return;
242 }
243
2701e5ae 244 Bool_t eventTriggered = inputHandler->IsEventSelected() & AliVEvent::kMB;
a9017e49 245
eb9356d5 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
f3eb27f6 255 const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
eb9356d5 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
f3eb27f6 263 Bool_t eventVertex = (vtxESD != 0);
a9017e49 264
a9017e49 265 Double_t vtx[3];
f3eb27f6 266 if (vtxESD)
f55b30d2 267 {
f3eb27f6 268 vtxESD->GetXYZ(vtx);
f55b30d2 269 fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
270 }
271
2fa65f52 272 // post the data already here
273 PostData(0, fOutput);
69b09e3b 274
a9017e49 275 //const Float_t kPtCut = 0.3;
276
2fa65f52 277 // create list of (label, eta) tuples
278 Int_t inputCount = 0;
279 Int_t* labelArr = 0;
280 Float_t* etaArr = 0;
2701e5ae 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)
a9017e49 372 {
2701e5ae 373
eb9356d5 374 if (vtxESD)
a9017e49 375 {
2701e5ae 376 AliDebug(AliLog::kInfo, "Processing tracklets");
eb9356d5 377 // get tracklets
378 const AliMultiplicity* mult = fESD->GetMultiplicity();
379 if (!mult)
69b09e3b 380 {
eb9356d5 381 AliDebug(AliLog::kError, "AliMultiplicity not available");
382 return;
69b09e3b 383 }
eb9356d5 384
eb9356d5 385 Bool_t foundInEta10 = kFALSE;
2701e5ae 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.
eb9356d5 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))
2701e5ae 405 {
eb9356d5 406 Printf("Skipped tracklet!");
407 continue;
2701e5ae 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
eb9356d5 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
d29b31c5 433 for (Int_t j=0; j<3; j++)
eb9356d5 434 {
d29b31c5 435 if (vtx[2] > fMultiplicity->GetVertexBegin(j) && vtx[2] < fMultiplicity->GetVertexEnd(j))
436 fEta[j]->Fill(etaArr[inputCount]);
eb9356d5 437 }
69b09e3b 438
f55b30d2 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
eb9356d5 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;
2fa65f52 451 }
452 }
a9017e49 453
2fa65f52 454 // eta range for nMCTracksSpecies, nESDTracksSpecies
70fdd197 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// }
a9017e49 464
2fa65f52 465 if (!fReadMC) // Processing of ESD information
466 {
eb9356d5 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)
2fa65f52 472 {
eb9356d5 473 Float_t eta = etaArr[i];
a9017e49 474
eb9356d5 475 if (TMath::Abs(eta) < 0.5)
476 nESDTracks05++;
a9017e49 477
eb9356d5 478 if (TMath::Abs(eta) < 1.0)
479 nESDTracks10++;
a9017e49 480
eb9356d5 481 if (TMath::Abs(eta) < 1.3)
482 nESDTracks14++;
483 }
484
eb9356d5 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)
b3b260d1 492 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 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;
a9017e49 500 }
a9017e49 501
2fa65f52 502 AliMCEvent* mcEvent = eventHandler->MCEvent();
503 if (!mcEvent) {
504 Printf("ERROR: Could not retrieve MC event");
505 return;
506 }
a9017e49 507
2fa65f52 508 AliStack* stack = mcEvent->Stack();
509 if (!stack)
510 {
511 AliDebug(AliLog::kError, "Stack not available");
512 return;
513 }
69b09e3b 514
2fa65f52 515 AliHeader* header = mcEvent->Header();
516 if (!header)
517 {
518 AliDebug(AliLog::kError, "Header not available");
519 return;
520 }
a9017e49 521
f3eb27f6 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 }
69b09e3b 532
533 // get process information
eb9356d5 534 AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
f3eb27f6 535
2fa65f52 536 Bool_t processEvent = kTRUE;
537 if (fSelectProcessType > 0)
538 {
2fa65f52 539 processEvent = kFALSE;
a9017e49 540
2fa65f52 541 // non diffractive
69b09e3b 542 if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
2fa65f52 543 processEvent = kTRUE;
a9017e49 544
2fa65f52 545 // single diffractive
69b09e3b 546 if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
2fa65f52 547 processEvent = kTRUE;
a9017e49 548
2fa65f52 549 // double diffractive
69b09e3b 550 if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
2fa65f52 551 processEvent = kTRUE;
a9017e49 552
2fa65f52 553 if (!processEvent)
69b09e3b 554 Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
2fa65f52 555 }
a9017e49 556
2fa65f52 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;
69b09e3b 571 Int_t nMCTracks10 = 0;
b3b260d1 572 Int_t nMCTracks14 = 0;
2fa65f52 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));
a9017e49 583
2fa65f52 584 TParticle* particle = stack->Particle(iMc);
a9017e49 585
2fa65f52 586 if (!particle)
587 {
588 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
589 continue;
590 }
a9017e49 591
2fa65f52 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 }
a9017e49 599
2fa65f52 600 //printf("%d) OK ", iMC);
601 //particle->Print();
a9017e49 602
2fa65f52 603 //if (particle->Pt() < kPtCut)
604 // continue;
a9017e49 605
2fa65f52 606 Int_t particleWeight = 1;
a9017e49 607
2fa65f52 608 Float_t pt = particle->Pt();
a9017e49 609
2fa65f52 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)
a9017e49 615 {
2fa65f52 616 Int_t bin = fPtSpectrum->FindBin(pt);
617 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
a9017e49 618 {
2fa65f52 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 }
a9017e49 630 }
a9017e49 631 }
a9017e49 632
2fa65f52 633 //Printf("MC weight is: %d", particleWeight);
a9017e49 634
2fa65f52 635 if (TMath::Abs(particle->Eta()) < 0.5)
636 nMCTracks05 += particleWeight;
a9017e49 637
69b09e3b 638 if (TMath::Abs(particle->Eta()) < 1.0)
639 nMCTracks10 += particleWeight;
a9017e49 640
eb9356d5 641 if (TMath::Abs(particle->Eta()) < 1.3)
b3b260d1 642 nMCTracks14 += particleWeight;
a9017e49 643
2fa65f52 644 nMCTracksAll += particleWeight;
69b09e3b 645
646 if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
647 fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
a9017e49 648
2fa65f52 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 }
a9017e49 659
2fa65f52 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
a9017e49 667
b3b260d1 668 fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
a9017e49 669
eb9356d5 670 // ESD processing
671 Int_t nESDTracks05 = 0;
672 Int_t nESDTracks10 = 0;
673 Int_t nESDTracks14 = 0;
a9017e49 674
eb9356d5 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;
a9017e49 679
eb9356d5 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;
a9017e49 683
eb9356d5 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;
69b09e3b 687
eb9356d5 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;
a9017e49 691
eb9356d5 692 for (Int_t i=0; i<inputCount; ++i)
693 {
694 Float_t eta = etaArr[i];
695 Int_t label = labelArr[i];
a9017e49 696
eb9356d5 697 Int_t particleWeight = 1;
2fa65f52 698
eb9356d5 699 // in case of systematic study, weight according to the change of the pt spectrum
700 if (fPtSpectrum)
701 {
702 TParticle* mother = 0;
2fa65f52 703
eb9356d5 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);
2fa65f52 710
eb9356d5 711 // in case of pt study we do not count particles w/o label, because they cannot be scaled
712 if (!mother)
713 continue;
2fa65f52 714
eb9356d5 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)
2fa65f52 723 {
eb9356d5 724 Float_t random = gRandom->Uniform();
725 if (factor > 1 && random < factor - 1)
2fa65f52 726 {
eb9356d5 727 particleWeight = 2;
2fa65f52 728 }
eb9356d5 729 else if (factor < 1 && random < 1 - factor)
730 particleWeight = 0;
2fa65f52 731 }
732 }
eb9356d5 733 }
2fa65f52 734
eb9356d5 735 //Printf("ESD weight is: %d", particleWeight);
2fa65f52 736
eb9356d5 737 if (TMath::Abs(eta) < 0.5)
738 nESDTracks05 += particleWeight;
a9017e49 739
eb9356d5 740 if (TMath::Abs(eta) < 1.0)
741 nESDTracks10 += particleWeight;
a9017e49 742
eb9356d5 743 if (TMath::Abs(eta) < 1.3)
744 nESDTracks14 += particleWeight;
a9017e49 745
eb9356d5 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
2fa65f52 764 {
eb9356d5 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]++;
2fa65f52 778
eb9356d5 779 // double counting is not ok for efficiency study
2fa65f52 780
eb9356d5 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])
2fa65f52 783 {
2fa65f52 784 if (TMath::Abs(eta) < etaRange)
eb9356d5 785 nESDTracksSpecies[6]++;
2fa65f52 786 }
69b09e3b 787 else
2fa65f52 788 {
eb9356d5 789 foundTracks[label] = kTRUE;
69b09e3b 790
eb9356d5 791 // particle (primary) already counted?
792 if (foundPrimaries[motherLabel])
69b09e3b 793 {
794 if (TMath::Abs(eta) < etaRange)
eb9356d5 795 nESDTracksSpecies[5]++;
69b09e3b 796 }
797 else
eb9356d5 798 foundPrimaries[motherLabel] = kTRUE;
2fa65f52 799 }
69b09e3b 800 }
eb9356d5 801 }
69b09e3b 802
eb9356d5 803 if (fParticleCorrection[0])
804 {
805 if (label >= 0 && stack->IsPhysicalPrimary(label))
69b09e3b 806 {
eb9356d5 807 TParticle* particle = stack->Particle(label);
2fa65f52 808
eb9356d5 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 }
2fa65f52 818
eb9356d5 819 // todo check if values are not completely off??
2fa65f52 820
eb9356d5 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());
69b09e3b 826 }
827 }
828 }
eb9356d5 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++)
69b09e3b 840 {
eb9356d5 841 // already found
842 if (foundPrimaries2[i])
843 continue;
69b09e3b 844
eb9356d5 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;
69b09e3b 871
eb9356d5 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 }
69b09e3b 881
eb9356d5 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)
69b09e3b 891 {
eb9356d5 892 // decayed
69b09e3b 893
eb9356d5 894 Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
69b09e3b 895 particle->Print();
2a311c4d 896 Printf("Daugthers:");
eb9356d5 897 for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
898 stack->Particle(d)->Print();
2a311c4d 899 Printf(" ");
69b09e3b 900
eb9356d5 901 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
2fa65f52 902 continue;
903 }
eb9356d5 904 }
905
906 if (trackref->DetectorId() == -1)
907 {
908 // stopped
909 Printf("Particle %d stopped at %f:", i, trackref->R());
69b09e3b 910 particle->Print();
2a311c4d 911 Printf(" ");
eb9356d5 912
913 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
914 continue;
2fa65f52 915 }
eb9356d5 916
917 Printf("Particle %d simply not tracked", i);
918 particle->Print();
2a311c4d 919 Printf(" ");
2fa65f52 920 }
eb9356d5 921 }
922
923 delete[] foundTracks;
924 delete[] foundPrimaries;
925 delete[] foundPrimaries2;
a9017e49 926
eb9356d5 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// }
a9017e49 932
eb9356d5 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 {
2fa65f52 946 // fill response matrix using vtxMC (best guess)
eb9356d5 947 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks14, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks14);
a9017e49 948
b3b260d1 949 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 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 }
a9017e49 956
69b09e3b 957 if (etaArr)
958 delete[] etaArr;
959 if (labelArr)
960 delete[] labelArr;
2701e5ae 961
a9017e49 962}
963
964void 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"));
69b09e3b 977 for (Int_t i = 0; i < 8; ++i)
a9017e49 978 fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
979 fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
69b09e3b 980
981 fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
a9017e49 982
983 if (!fMultiplicity)
984 {
985 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
986 return;
987 }
988
eb9356d5 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");
a9017e49 998
999 fMultiplicity->SaveHistograms();
69b09e3b 1000 for (Int_t i = 0; i < 8; ++i)
a9017e49 1001 if (fParticleCorrection[i])
1002 fParticleCorrection[i]->SaveHistograms();
1003 if (fParticleSpecies)
1004 fParticleSpecies->Write();
69b09e3b 1005 if (fdNdpT)
1006 fdNdpT->Write();
a9017e49 1007
eb9356d5 1008 fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1009 if (fTemp1)
1010 fTemp1->Write();
1011
f55b30d2 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
eb9356d5 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
a9017e49 1034 TObjString option(fOption);
1035 option.Write();
1036
1037 file->Close();
1038
eb9356d5 1039 Printf("Written result to %s", fileName.Data());
a9017e49 1040}