Transition PWG0 -> PWGUD
[u/mrichter/AliRoot.git] / PWGUD / selectors / 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"
87368ddd 34#include "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 {
335a5e1b 249 AliPhysicsSelection* physicsSelection = static_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
eb9356d5 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");
fd396b1a 382 if (labelArr) delete[] labelArr;
383 if (etaArr) delete[] etaArr;
eb9356d5 384 return;
69b09e3b 385 }
eb9356d5 386
eb9356d5 387 Bool_t foundInEta10 = kFALSE;
ee6ee43d 388 if ((fAnalysisMode & AliPWG0Helper::kSPD) && !(fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS || fAnalysisMode & AliPWG0Helper::kTPCSPD)) {
2701e5ae 389 // if we are counting both tracks and tracklets, these arrays were already initialized above
390 AliDebug(AliLog::kInfo, "Booking arrays");
fd396b1a 391 if (!labelArr) labelArr = new Int_t[mult->GetNumberOfTracklets()];
392 if (!etaArr) etaArr = new Float_t[mult->GetNumberOfTracklets()];
2701e5ae 393 }
394 if (inputCount) foundInEta10 = kTRUE; // by definition, if we found a track.
eb9356d5 395
396 // get multiplicity from ITS tracklets
397 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
398 {
399 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
400
401 Float_t deltaPhi = mult->GetDeltaPhi(i);
402
403 if (fDeltaPhiCut > 0 && TMath::Abs(deltaPhi) > fDeltaPhiCut)
404 continue;
405
406 if (fSystSkipParticles && gRandom->Uniform() < (0.0153))
2701e5ae 407 {
eb9356d5 408 Printf("Skipped tracklet!");
409 continue;
2701e5ae 410 }
411
412 // if counting tracks+tracklets, check if clusters where already used in tracks
413 if (fAnalysisMode & AliPWG0Helper::kTPCSPD) {
414 Int_t id1,id2;
415 mult->GetTrackletTrackIDs(i,0,id1,id2);
416 if ( (id1>=0&& !fESD->GetTrack(id1)->TestBit(kRejBit) ) ||
417 (id2>=0&& !fESD->GetTrack(id2)->TestBit(kRejBit) )
418 ) {
419 printf ("Skipping tracklet: already used in track");
420 continue;
421 }
422 // Ruben also had this, but we're not using pure ITSSA tracks here:
423 // if (fSPDMult->FreeClustersTracklet(itr,1)) trITSSApure++; // not used in ITS_SA_Pure track
424 }
425
426
eb9356d5 427 etaArr[inputCount] = mult->GetEta(i);
428 if (mult->GetLabel(i, 0) == mult->GetLabel(i, 1))
429 {
430 labelArr[inputCount] = mult->GetLabel(i, 0);
431 }
432 else
433 labelArr[inputCount] = -1;
434
d29b31c5 435 for (Int_t j=0; j<3; j++)
eb9356d5 436 {
d29b31c5 437 if (vtx[2] > fMultiplicity->GetVertexBegin(j) && vtx[2] < fMultiplicity->GetVertexEnd(j))
438 fEta[j]->Fill(etaArr[inputCount]);
eb9356d5 439 }
69b09e3b 440
f55b30d2 441 if (vtxESD && TMath::Abs(vtx[2]) < 10)
442 fEtaPhi->Fill(etaArr[inputCount], mult->GetPhi(i));
443
444 // we have to repeat the trigger here, because the tracklet might have been kicked out fSystSkipParticles
eb9356d5 445 if (TMath::Abs(etaArr[inputCount]) < 1)
446 foundInEta10 = kTRUE;
447
448 ++inputCount;
449 }
450
451 if (fSystSkipParticles && (fTrigger & AliTriggerAnalysis::kOneParticle) && !foundInEta10)
452 eventTriggered = kFALSE;
2fa65f52 453 }
454 }
a9017e49 455
2fa65f52 456 // eta range for nMCTracksSpecies, nESDTracksSpecies
70fdd197 457 Float_t etaRange = 1.0;
458// switch (fAnalysisMode) {
459// case AliPWG0Helper::kInvalid: break;
460// case AliPWG0Helper::kTPC:
461// case AliPWG0Helper::kTPCITS:
462// etaRange = 1.0; break;
463// case AliPWG0Helper::kSPD: etaRange = 1.0; break;
464// default: break;
465// }
a9017e49 466
2fa65f52 467 if (!fReadMC) // Processing of ESD information
468 {
eb9356d5 469 Int_t nESDTracks05 = 0;
470 Int_t nESDTracks10 = 0;
471 Int_t nESDTracks14 = 0;
472
473 for (Int_t i=0; i<inputCount; ++i)
2fa65f52 474 {
eb9356d5 475 Float_t eta = etaArr[i];
a9017e49 476
eb9356d5 477 if (TMath::Abs(eta) < 0.5)
478 nESDTracks05++;
a9017e49 479
eb9356d5 480 if (TMath::Abs(eta) < 1.0)
481 nESDTracks10++;
a9017e49 482
eb9356d5 483 if (TMath::Abs(eta) < 1.3)
484 nESDTracks14++;
485 }
486
eb9356d5 487 //if (nESDTracks05 >= 20 || nESDTracks10 >= 30 || nESDTracks14 >= 32)
488 // 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);
489
490 if (eventTriggered)
491 fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
492
493 if (eventTriggered && eventVertex)
b3b260d1 494 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 495 }
496 else if (fReadMC) // Processing of MC information
497 {
498 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
499 if (!eventHandler) {
500 Printf("ERROR: Could not retrieve MC event handler");
fd396b1a 501 if (labelArr) delete[] labelArr;
502 if (etaArr) delete[] etaArr;
2fa65f52 503 return;
a9017e49 504 }
a9017e49 505
2fa65f52 506 AliMCEvent* mcEvent = eventHandler->MCEvent();
507 if (!mcEvent) {
508 Printf("ERROR: Could not retrieve MC event");
fd396b1a 509 if (labelArr) delete[] labelArr;
510 if (etaArr) delete[] etaArr;
2fa65f52 511 return;
512 }
a9017e49 513
2fa65f52 514 AliStack* stack = mcEvent->Stack();
515 if (!stack)
516 {
517 AliDebug(AliLog::kError, "Stack not available");
fd396b1a 518 if (labelArr) delete[] labelArr;
519 if (etaArr) delete[] etaArr;
2fa65f52 520 return;
521 }
69b09e3b 522
2fa65f52 523 AliHeader* header = mcEvent->Header();
524 if (!header)
525 {
526 AliDebug(AliLog::kError, "Header not available");
fd396b1a 527 if (labelArr) delete[] labelArr;
528 if (etaArr) delete[] etaArr;
2fa65f52 529 return;
530 }
a9017e49 531
f3eb27f6 532 if (fUseMCVertex)
533 {
534 Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
535 // get the MC vertex
536 AliGenEventHeader* genHeader = header->GenEventHeader();
537 TArrayF vtxMC(3);
538 genHeader->PrimaryVertex(vtxMC);
539
540 vtx[2] = vtxMC[2];
541 }
69b09e3b 542
543 // get process information
eb9356d5 544 AliPWG0Helper::MCProcessType processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
f3eb27f6 545
2fa65f52 546 Bool_t processEvent = kTRUE;
547 if (fSelectProcessType > 0)
548 {
2fa65f52 549 processEvent = kFALSE;
a9017e49 550
2fa65f52 551 // non diffractive
69b09e3b 552 if (fSelectProcessType == 1 && processType == AliPWG0Helper::kND)
2fa65f52 553 processEvent = kTRUE;
a9017e49 554
2fa65f52 555 // single diffractive
69b09e3b 556 if (fSelectProcessType == 2 && processType == AliPWG0Helper::kSD)
2fa65f52 557 processEvent = kTRUE;
a9017e49 558
2fa65f52 559 // double diffractive
69b09e3b 560 if (fSelectProcessType == 3 && processType == AliPWG0Helper::kDD)
2fa65f52 561 processEvent = kTRUE;
a9017e49 562
2fa65f52 563 if (!processEvent)
69b09e3b 564 Printf("Skipping this event, because it is not of the requested process type (%d)", (Int_t) processType);
2fa65f52 565 }
a9017e49 566
2fa65f52 567 if (processEvent)
568 {
569 // get the MC vertex
570 AliGenEventHeader* genHeader = header->GenEventHeader();
571 TArrayF vtxMC(3);
572 genHeader->PrimaryVertex(vtxMC);
573
574 // get number of tracks from MC
575 // loop over mc particles
576 Int_t nPrim = stack->GetNprimary();
577 Int_t nMCPart = stack->GetNtrack();
578
579 // tracks in different eta ranges
580 Int_t nMCTracks05 = 0;
69b09e3b 581 Int_t nMCTracks10 = 0;
b3b260d1 582 Int_t nMCTracks14 = 0;
2fa65f52 583 Int_t nMCTracksAll = 0;
584
585 // tracks per particle species, in |eta| < 2 (systematic study)
586 Int_t nMCTracksSpecies[4]; // (pi, K, p, other)
587 for (Int_t i = 0; i<4; ++i)
588 nMCTracksSpecies[i] = 0;
589
590 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
591 {
592 AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
a9017e49 593
2fa65f52 594 TParticle* particle = stack->Particle(iMc);
a9017e49 595
2fa65f52 596 if (!particle)
597 {
598 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
599 continue;
600 }
a9017e49 601
2fa65f52 602 Bool_t debug = kFALSE;
603 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim, debug) == kFALSE)
604 {
605 //printf("%d) DROPPED ", iMC);
606 //particle->Print();
607 continue;
608 }
a9017e49 609
2fa65f52 610 //printf("%d) OK ", iMC);
611 //particle->Print();
a9017e49 612
2fa65f52 613 //if (particle->Pt() < kPtCut)
614 // continue;
a9017e49 615
2fa65f52 616 Int_t particleWeight = 1;
a9017e49 617
2fa65f52 618 Float_t pt = particle->Pt();
a9017e49 619
2fa65f52 620 // in case of systematic study, weight according to the change of the pt spectrum
621 // it cannot be just multiplied because we cannot count "half of a particle"
622 // instead a random generator decides if the particle is counted twice (if value > 1)
623 // or not (if value < 0)
624 if (fPtSpectrum)
a9017e49 625 {
2fa65f52 626 Int_t bin = fPtSpectrum->FindBin(pt);
627 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
a9017e49 628 {
2fa65f52 629 Float_t factor = fPtSpectrum->GetBinContent(bin);
630 if (factor > 0)
631 {
632 Float_t random = gRandom->Uniform();
633 if (factor > 1 && random < factor - 1)
634 {
635 particleWeight = 2;
636 }
637 else if (factor < 1 && random < 1 - factor)
638 particleWeight = 0;
639 }
a9017e49 640 }
a9017e49 641 }
a9017e49 642
2fa65f52 643 //Printf("MC weight is: %d", particleWeight);
a9017e49 644
2fa65f52 645 if (TMath::Abs(particle->Eta()) < 0.5)
646 nMCTracks05 += particleWeight;
a9017e49 647
69b09e3b 648 if (TMath::Abs(particle->Eta()) < 1.0)
649 nMCTracks10 += particleWeight;
a9017e49 650
eb9356d5 651 if (TMath::Abs(particle->Eta()) < 1.3)
b3b260d1 652 nMCTracks14 += particleWeight;
a9017e49 653
2fa65f52 654 nMCTracksAll += particleWeight;
69b09e3b 655
656 if (particle->Pt() > 0 && TMath::Abs(particle->Eta()) < 1.0)
657 fdNdpT->Fill(particle->Pt(), 1.0 / particle->Pt());
a9017e49 658
2fa65f52 659 if (fParticleCorrection[0] || fParticleSpecies)
660 {
661 Int_t id = -1;
662 switch (TMath::Abs(particle->GetPdgCode()))
663 {
664 case 211: id = 0; break;
665 case 321: id = 1; break;
666 case 2212: id = 2; break;
667 default: id = 3; break;
668 }
a9017e49 669
2fa65f52 670 if (TMath::Abs(particle->Eta()) < etaRange)
671 nMCTracksSpecies[id]++;
672
673 if (fParticleCorrection[id])
674 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
675 }
676 } // end of mc particle
a9017e49 677
b3b260d1 678 fMultiplicity->FillGenerated(vtxMC[2], eventTriggered, eventVertex, processType, (Int_t) nMCTracks05, (Int_t) nMCTracks10, (Int_t) nMCTracks14, (Int_t) nMCTracksAll);
a9017e49 679
eb9356d5 680 // ESD processing
681 Int_t nESDTracks05 = 0;
682 Int_t nESDTracks10 = 0;
683 Int_t nESDTracks14 = 0;
a9017e49 684
eb9356d5 685 // tracks per particle species, in |eta| < 2 (systematic study)
686 Int_t nESDTracksSpecies[7]; // (pi, K, p, other, nolabel, doublecount_prim, doublecount_all)
687 for (Int_t i = 0; i<7; ++i)
688 nESDTracksSpecies[i] = 0;
a9017e49 689
7f86bd47 690 Bool_t* foundPrimaries = new Bool_t[nMCPart]; // to prevent double counting
eb9356d5 691 for (Int_t i=0; i<nPrim; i++)
692 foundPrimaries[i] = kFALSE;
a9017e49 693
7f86bd47 694 Bool_t* foundPrimaries2 = new Bool_t[nMCPart]; // to prevent double counting
eb9356d5 695 for (Int_t i=0; i<nPrim; i++)
696 foundPrimaries2[i] = kFALSE;
69b09e3b 697
eb9356d5 698 Bool_t* foundTracks = new Bool_t[nMCPart]; // to prevent double counting
699 for (Int_t i=0; i<nMCPart; i++)
700 foundTracks[i] = kFALSE;
a9017e49 701
eb9356d5 702 for (Int_t i=0; i<inputCount; ++i)
703 {
704 Float_t eta = etaArr[i];
705 Int_t label = labelArr[i];
a9017e49 706
eb9356d5 707 Int_t particleWeight = 1;
2fa65f52 708
eb9356d5 709 // in case of systematic study, weight according to the change of the pt spectrum
710 if (fPtSpectrum)
711 {
712 TParticle* mother = 0;
2fa65f52 713
eb9356d5 714 // preserve label for later
715 Int_t labelCopy = label;
716 if (labelCopy >= 0)
717 labelCopy = AliPWG0Helper::FindPrimaryMotherLabel(stack, labelCopy);
718 if (labelCopy >= 0)
719 mother = stack->Particle(labelCopy);
2fa65f52 720
eb9356d5 721 // in case of pt study we do not count particles w/o label, because they cannot be scaled
722 if (!mother)
723 continue;
2fa65f52 724
eb9356d5 725 // it cannot be just multiplied because we cannot count "half of a particle"
726 // instead a random generator decides if the particle is counted twice (if value > 1)
727 // or not (if value < 0)
728 Int_t bin = fPtSpectrum->FindBin(mother->Pt());
729 if (bin > 0 && bin <= fPtSpectrum->GetNbinsX())
730 {
731 Float_t factor = fPtSpectrum->GetBinContent(bin);
732 if (factor > 0)
2fa65f52 733 {
eb9356d5 734 Float_t random = gRandom->Uniform();
735 if (factor > 1 && random < factor - 1)
2fa65f52 736 {
eb9356d5 737 particleWeight = 2;
2fa65f52 738 }
eb9356d5 739 else if (factor < 1 && random < 1 - factor)
740 particleWeight = 0;
2fa65f52 741 }
742 }
eb9356d5 743 }
2fa65f52 744
eb9356d5 745 //Printf("ESD weight is: %d", particleWeight);
2fa65f52 746
eb9356d5 747 if (TMath::Abs(eta) < 0.5)
748 nESDTracks05 += particleWeight;
a9017e49 749
eb9356d5 750 if (TMath::Abs(eta) < 1.0)
751 nESDTracks10 += particleWeight;
a9017e49 752
eb9356d5 753 if (TMath::Abs(eta) < 1.3)
754 nESDTracks14 += particleWeight;
a9017e49 755
eb9356d5 756 if (fParticleSpecies)
757 {
758 Int_t motherLabel = -1;
759 TParticle* mother = 0;
760
761 // find mother
762 if (label >= 0)
763 motherLabel = AliPWG0Helper::FindPrimaryMotherLabel(stack, label);
764 if (motherLabel >= 0)
765 mother = stack->Particle(motherLabel);
766
767 if (!mother)
768 {
769 // count tracks that did not have a label
770 if (TMath::Abs(eta) < etaRange)
771 nESDTracksSpecies[4]++;
772 }
773 else
2fa65f52 774 {
eb9356d5 775 // get particle type (pion, proton, kaon, other) of mother
776 Int_t idMother = -1;
777 switch (TMath::Abs(mother->GetPdgCode()))
778 {
779 case 211: idMother = 0; break;
780 case 321: idMother = 1; break;
781 case 2212: idMother = 2; break;
782 default: idMother = 3; break;
783 }
784
785 // double counting is ok for particle ratio study
786 if (TMath::Abs(eta) < etaRange)
787 nESDTracksSpecies[idMother]++;
2fa65f52 788
eb9356d5 789 // double counting is not ok for efficiency study
2fa65f52 790
eb9356d5 791 // 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)
792 if (foundTracks[label])
2fa65f52 793 {
2fa65f52 794 if (TMath::Abs(eta) < etaRange)
eb9356d5 795 nESDTracksSpecies[6]++;
2fa65f52 796 }
69b09e3b 797 else
2fa65f52 798 {
eb9356d5 799 foundTracks[label] = kTRUE;
69b09e3b 800
eb9356d5 801 // particle (primary) already counted?
802 if (foundPrimaries[motherLabel])
69b09e3b 803 {
804 if (TMath::Abs(eta) < etaRange)
eb9356d5 805 nESDTracksSpecies[5]++;
69b09e3b 806 }
807 else
eb9356d5 808 foundPrimaries[motherLabel] = kTRUE;
2fa65f52 809 }
69b09e3b 810 }
eb9356d5 811 }
69b09e3b 812
eb9356d5 813 if (fParticleCorrection[0])
814 {
815 if (label >= 0 && stack->IsPhysicalPrimary(label))
69b09e3b 816 {
eb9356d5 817 TParticle* particle = stack->Particle(label);
2fa65f52 818
eb9356d5 819 // get particle type (pion, proton, kaon, other)
820 Int_t id = -1;
821 switch (TMath::Abs(particle->GetPdgCode()))
822 {
823 case 211: id = 0; break;
824 case 321: id = 1; break;
825 case 2212: id = 2; break;
826 default: id = 3; break;
827 }
2fa65f52 828
eb9356d5 829 // todo check if values are not completely off??
2fa65f52 830
eb9356d5 831 // particle (primary) already counted?
832 if (!foundPrimaries2[label])
833 {
834 foundPrimaries2[label] = kTRUE;
835 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
69b09e3b 836 }
837 }
838 }
eb9356d5 839 }
840
841 if (fParticleCorrection[0])
842 {
843 // if the particle decays/stops before this radius we do not see it
844 // 8cm larger than SPD layer 2
845 // 123cm TPC radius where a track has about 50 clusters (cut limit)
846 const Float_t endRadius = (fAnalysisMode & AliPWG0Helper::kSPD) ? 8. : 123;
847
848 // loop over all primaries that have not been found
849 for (Int_t i=0; i<nPrim; i++)
69b09e3b 850 {
eb9356d5 851 // already found
852 if (foundPrimaries2[i])
853 continue;
69b09e3b 854
eb9356d5 855 TParticle* particle = 0;
856 TClonesArray* trackrefs = 0;
857 mcEvent->GetParticleAndTR(i, particle, trackrefs);
858
859 // true primary and charged
860 if (!AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
861 continue;
862
863 //skip particles with larger |eta| than 3, to keep the log clean, is anyway not included in correction map
864 if (TMath::Abs(particle->Eta()) > 3)
865 continue;
866
867 // skipping checking of process type of daughter: Neither kPBrem, kPDeltaRay nor kPCerenkov should appear in the event generation
868
869 // get particle type (pion, proton, kaon, other)
870 Int_t id = -1;
871 switch (TMath::Abs(particle->GetPdgCode()))
872 {
873 case 211: id = 4; break;
874 case 321: id = 5; break;
875 case 2212: id = 6; break;
876 default: id = 7; break;
877 }
878
879 if (!fParticleCorrection[id])
880 continue;
69b09e3b 881
eb9356d5 882 // get last track reference
883 AliTrackReference* trackref = dynamic_cast<AliTrackReference*> (trackrefs->Last());
884
885 if (!trackref)
886 {
887 Printf("ERROR: Could not get trackref of %d (count %d)", i, trackrefs->GetEntries());
888 particle->Print();
889 continue;
890 }
69b09e3b 891
eb9356d5 892 // particle in tracking volume long enough...
893 if (trackref->R() > endRadius)
894 continue;
895
896 if (particle->GetLastDaughter() >= 0)
897 {
898 Int_t uID = stack->Particle(particle->GetLastDaughter())->GetUniqueID();
899 //if (uID != kPBrem && uID != kPDeltaRay && uID < kPCerenkov)
900 if (uID == kPDecay)
69b09e3b 901 {
eb9356d5 902 // decayed
69b09e3b 903
eb9356d5 904 Printf("Particle %d (%s) decayed at %f, daugher uniqueID: %d:", i, particle->GetName(), trackref->R(), uID);
69b09e3b 905 particle->Print();
2a311c4d 906 Printf("Daugthers:");
eb9356d5 907 for (Int_t d = particle->GetFirstDaughter(); d <= particle->GetLastDaughter(); d++)
908 stack->Particle(d)->Print();
2a311c4d 909 Printf(" ");
69b09e3b 910
eb9356d5 911 fParticleCorrection[id]->GetTrackCorrection()->FillGene(vtxMC[2], particle->Eta(), particle->Pt());
2fa65f52 912 continue;
913 }
eb9356d5 914 }
915
916 if (trackref->DetectorId() == -1)
917 {
918 // stopped
919 Printf("Particle %d stopped at %f:", i, trackref->R());
69b09e3b 920 particle->Print();
2a311c4d 921 Printf(" ");
eb9356d5 922
923 fParticleCorrection[id]->GetTrackCorrection()->FillMeas(vtxMC[2], particle->Eta(), particle->Pt());
924 continue;
2fa65f52 925 }
eb9356d5 926
927 Printf("Particle %d simply not tracked", i);
928 particle->Print();
2a311c4d 929 Printf(" ");
2fa65f52 930 }
eb9356d5 931 }
932
933 delete[] foundTracks;
934 delete[] foundPrimaries;
935 delete[] foundPrimaries2;
a9017e49 936
eb9356d5 937// if ((Int_t) nMCTracks14 > 10 && nESDTracks14 <= 3)
938// {
939// TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
940// 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);
941// }
a9017e49 942
eb9356d5 943 if (eventTriggered)
944 {
945 fMultiplicity->FillTriggeredEvent(nESDTracks05, nESDTracks10, nESDTracks14);
946 fMultiplicity->FillNoVertexEvent(vtxMC[2], eventVertex, nMCTracks05, nMCTracks10, nMCTracks14, nESDTracks05, nESDTracks10, nESDTracks14);
947// if (!eventVertex)
948// {
949// if (nESDTracks05 == 0)
950// fTemp1->Fill(nMCTracks05, nESDTracks05);
951// }
952 }
953
954 if (eventTriggered && eventVertex)
955 {
2fa65f52 956 // fill response matrix using vtxMC (best guess)
eb9356d5 957 fMultiplicity->FillCorrection(vtxMC[2], nMCTracks05, nMCTracks10, nMCTracks14, nMCTracksAll, nESDTracks05, nESDTracks10, nESDTracks14);
a9017e49 958
b3b260d1 959 fMultiplicity->FillMeasured(vtx[2], nESDTracks05, nESDTracks10, nESDTracks14);
2fa65f52 960
961 if (fParticleSpecies)
962 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]);
963 }
964 }
965 }
a9017e49 966
69b09e3b 967 if (etaArr)
968 delete[] etaArr;
969 if (labelArr)
970 delete[] labelArr;
2701e5ae 971
a9017e49 972}
973
974void AliMultiplicityTask::Terminate(Option_t *)
975{
976 // The Terminate() function is the last function to be called during
977 // a query. It always runs on the client, it can be used to present
978 // the results graphically or save the results to file.
979
980 fOutput = dynamic_cast<TList*> (GetOutputData(0));
981 if (!fOutput) {
982 Printf("ERROR: fOutput not available");
983 return;
984 }
985
986 fMultiplicity = dynamic_cast<AliMultiplicityCorrection*> (fOutput->FindObject("Multiplicity"));
69b09e3b 987 for (Int_t i = 0; i < 8; ++i)
a9017e49 988 fParticleCorrection[i] = dynamic_cast<AliCorrection*> (fOutput->FindObject(Form("correction_%d", i)));
989 fParticleSpecies = dynamic_cast<TNtuple*> (fOutput->FindObject("fParticleSpecies"));
69b09e3b 990
991 fdNdpT = dynamic_cast<TH1*> (fOutput->FindObject("fdNdpT"));
a9017e49 992
993 if (!fMultiplicity)
994 {
995 AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fMultiplicity));
996 return;
997 }
998
eb9356d5 999 TString fileName("multiplicity");
1000 if (fSelectProcessType == 1)
1001 fileName += "ND";
1002 if (fSelectProcessType == 2)
1003 fileName += "SD";
1004 if (fSelectProcessType == 3)
1005 fileName += "DD";
1006 fileName += ".root";
1007 TFile* file = TFile::Open(fileName, "RECREATE");
a9017e49 1008
1009 fMultiplicity->SaveHistograms();
69b09e3b 1010 for (Int_t i = 0; i < 8; ++i)
a9017e49 1011 if (fParticleCorrection[i])
1012 fParticleCorrection[i]->SaveHistograms();
1013 if (fParticleSpecies)
1014 fParticleSpecies->Write();
69b09e3b 1015 if (fdNdpT)
1016 fdNdpT->Write();
a9017e49 1017
eb9356d5 1018 fTemp1 = dynamic_cast<TH1*> (fOutput->FindObject("fTemp1"));
1019 if (fTemp1)
1020 fTemp1->Write();
1021
f55b30d2 1022 fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
1023 if (fEtaPhi)
1024 fEtaPhi->Write();
1025
1026 fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
1027 if (fVertex)
1028 fVertex->Write();
1029
eb9356d5 1030 for (Int_t i=0; i<3; i++)
1031 {
1032 fEta[i] = dynamic_cast<TH1*> (fOutput->FindObject(Form("fEta_%d", i)));
1033 if (fEta[i])
1034 {
1035 fEta[i]->Sumw2();
1036 Float_t events = fMultiplicity->GetMultiplicityESD(i)->Integral(1, fMultiplicity->GetMultiplicityESD(i)->GetNbinsX());
1037 if (events > 0)
1038 fEta[i]->Scale(1.0 / events);
1039 fEta[i]->Scale(1.0 / fEta[i]->GetXaxis()->GetBinWidth(1));
1040 fEta[i]->Write();
1041 }
1042 }
1043
a9017e49 1044 TObjString option(fOption);
1045 option.Write();
1046
1047 file->Close();
1048
eb9356d5 1049 Printf("Written result to %s", fileName.Data());
a9017e49 1050}