fixing warnings
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaCorrectionTask.cxx
CommitLineData
0f67a57c 1/* $Id$ */
2
3#include "AlidNdEtaCorrectionTask.h"
4
54b096ef 5#include <TROOT.h>
0f67a57c 6#include <TCanvas.h>
7#include <TChain.h>
8#include <TFile.h>
9#include <TH1F.h>
54b096ef 10#include <TH2F.h>
11#include <TProfile.h>
0f67a57c 12#include <TParticle.h>
13
14#include <AliLog.h>
15#include <AliESDVertex.h>
16#include <AliESDEvent.h>
17#include <AliStack.h>
18#include <AliHeader.h>
19#include <AliGenEventHeader.h>
20#include <AliMultiplicity.h>
21#include <AliAnalysisManager.h>
22#include <AliMCEventHandler.h>
23#include <AliMCEvent.h>
24#include <AliESDInputHandler.h>
25
745d6088 26#include "AliESDtrackCuts.h"
0f67a57c 27#include "AliPWG0Helper.h"
28//#include "AliCorrection.h"
29//#include "AliCorrectionMatrix3D.h"
30#include "dNdEta/dNdEtaAnalysis.h"
31#include "dNdEta/AlidNdEtaCorrection.h"
54b096ef 32#include "AliVertexerTracks.h"
0f67a57c 33
34ClassImp(AlidNdEtaCorrectionTask)
35
36AlidNdEtaCorrectionTask::AlidNdEtaCorrectionTask(const char* opt) :
37 AliAnalysisTask("AlidNdEtaCorrectionTask", ""),
38 fESD(0),
39 fOutput(0),
40 fOption(opt),
770a1f1d 41 fAnalysisMode(AliPWG0Helper::kTPC),
0fc41645 42 fTrigger(AliPWG0Helper::kMB1),
0f67a57c 43 fSignMode(0),
3d7758c1 44 fOnlyPrimaries(kFALSE),
0f67a57c 45 fEsdTrackCuts(0),
46 fdNdEtaCorrection(0),
47 fdNdEtaAnalysisMC(0),
48 fdNdEtaAnalysisESD(0),
49 fPIDParticles(0),
96fcc8a7 50 fPIDTracks(0),
54b096ef 51 fVertexCorrelation(0),
52 fVertexProfile(0),
53 fVertexShiftNorm(0),
3d7758c1 54 fEtaCorrelation(0),
55 fEtaProfile(0),
96fcc8a7 56 fSigmaVertexTracks(0),
4ef701a0 57 fSigmaVertexPrim(0),
58 fMultAll(0),
59 fMultTr(0),
50ec344d 60 fMultVtx(0),
61 fEventStats(0)
0f67a57c 62{
63 //
64 // Constructor. Initialization of pointers
65 //
66
67 // Define input and output slots here
68 DefineInput(0, TChain::Class());
69 DefineOutput(0, TList::Class());
96fcc8a7 70
71 for (Int_t i=0; i<2; i++)
72 fdNdEtaCorrectionProcessType[i] = 0;
4ef701a0 73
3d7758c1 74 for (Int_t i=0; i<8; i++)
4ef701a0 75 fDeltaPhi[i] = 0;
0f67a57c 76}
77
78AlidNdEtaCorrectionTask::~AlidNdEtaCorrectionTask()
79{
80 //
81 // Destructor
82 //
83
84 // histograms are in the output list and deleted when the output
85 // list is deleted by the TSelector dtor
86
87 if (fOutput) {
88 delete fOutput;
89 fOutput = 0;
90 }
91}
92
93//________________________________________________________________________
94void AlidNdEtaCorrectionTask::ConnectInputData(Option_t *)
95{
96 // Connect ESD
97 // Called once
98
99 Printf("AlidNdEtaCorrectionTask::ConnectInputData called");
100
101 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
102 if (!tree) {
103 Printf("ERROR: Could not read tree from input slot 0");
104 } else {
105 // Disable all branches and enable only the needed ones
96fcc8a7 106 //tree->SetBranchStatus("*", 0);
0f67a57c 107
108 tree->SetBranchStatus("fTriggerMask", 1);
109 tree->SetBranchStatus("fSPDVertex*", 1);
770a1f1d 110 // PrimaryVertex
0f67a57c 111
770a1f1d 112 if (fAnalysisMode == AliPWG0Helper::kSPD)
0f67a57c 113 tree->SetBranchStatus("fSPDMult*", 1);
114
770a1f1d 115 if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
0f67a57c 116 AliESDtrackCuts::EnableNeededBranches(tree);
117 tree->SetBranchStatus("fTracks.fLabel", 1);
118 }
119
0f67a57c 120 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
121
122 if (!esdH) {
123 Printf("ERROR: Could not get ESDInputHandler");
124 } else
125 fESD = esdH->GetEvent();
126 }
96fcc8a7 127
128 // disable info messages of AliMCEvent (per event)
129 AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
0f67a57c 130}
131
132void AlidNdEtaCorrectionTask::CreateOutputObjects()
133{
134 // create result objects and add to output list
135
136 if (fOption.Contains("only-positive"))
137 {
138 Printf("INFO: Processing only positive particles.");
139 fSignMode = 1;
140 }
141 else if (fOption.Contains("only-negative"))
142 {
143 Printf("INFO: Processing only negative particles.");
144 fSignMode = -1;
145 }
146
0f67a57c 147 fOutput = new TList;
148 fOutput->SetOwner();
149
770a1f1d 150 fdNdEtaCorrection = new AlidNdEtaCorrection("dndeta_correction", "dndeta_correction", fAnalysisMode);
0f67a57c 151 fOutput->Add(fdNdEtaCorrection);
152
153 fPIDParticles = new TH1F("pid_particles", "PID of generated primary particles", 10001, -5000.5, 5000.5);
154 fOutput->Add(fPIDParticles);
155
156 fPIDTracks = new TH1F("pid_tracks", "MC PID of reconstructed tracks", 10001, -5000.5, 5000.5);
157 fOutput->Add(fPIDTracks);
158
770a1f1d 159 fdNdEtaAnalysisMC = new dNdEtaAnalysis("dndetaMC", "dndetaMC", fAnalysisMode);
0f67a57c 160 fOutput->Add(fdNdEtaAnalysisMC);
161
770a1f1d 162 fdNdEtaAnalysisESD = new dNdEtaAnalysis("dndetaESD", "dndetaESD", fAnalysisMode);
0f67a57c 163 fOutput->Add(fdNdEtaAnalysisESD);
96fcc8a7 164
165 if (fOption.Contains("process-types")) {
54b096ef 166 fdNdEtaCorrectionProcessType[0] = new AlidNdEtaCorrection("dndeta_correction_ND", "dndeta_correction_ND", fAnalysisMode);
167 fdNdEtaCorrectionProcessType[1] = new AlidNdEtaCorrection("dndeta_correction_SD", "dndeta_correction_SD", fAnalysisMode);
168 fdNdEtaCorrectionProcessType[2] = new AlidNdEtaCorrection("dndeta_correction_DD", "dndeta_correction_DD", fAnalysisMode);
96fcc8a7 169
170 fOutput->Add(fdNdEtaCorrectionProcessType[0]);
171 fOutput->Add(fdNdEtaCorrectionProcessType[1]);
172 fOutput->Add(fdNdEtaCorrectionProcessType[2]);
173 }
174
175 if (fOption.Contains("sigma-vertex"))
176 {
177 fSigmaVertexTracks = new TH1F("fSigmaVertexTracks", "fSigmaVertexTracks;Nsigma2vertex;NacceptedTracks", 60, 0.05, 6.05);
178 fSigmaVertexPrim = new TH1F("fSigmaVertexPrim", "fSigmaVertexPrim;Nsigma2vertex;NacceptedPrimaries", 60, 0.05, 6.05);
179 fOutput->Add(fSigmaVertexTracks);
180 fOutput->Add(fSigmaVertexPrim);
181 Printf("WARNING: sigma-vertex analysis enabled. This will produce weird results in the AliESDtrackCuts histograms");
182 }
54b096ef 183
0fc41645 184 fVertexCorrelation = new TH2F("fVertexCorrelation", "fVertexCorrelation;MC z-vtx;ESD z-vtx", 120, -30, 30, 120, -30, 30);
185 fOutput->Add(fVertexCorrelation);
54b096ef 186 fVertexProfile = new TProfile("fVertexProfile", "fVertexProfile;MC z-vtx;MC z-vtx - ESD z-vtx", 40, -20, 20);
0fc41645 187 fOutput->Add(fVertexProfile);
54b096ef 188 fVertexShiftNorm = new TH1F("fVertexShiftNorm", "fVertexShiftNorm;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries", 200, -100, 100);
0fc41645 189 fOutput->Add(fVertexShiftNorm);
190
3d7758c1 191 fEtaCorrelation = new TH2F("fEtaCorrelation", "fEtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);
0fc41645 192 fOutput->Add(fEtaCorrelation);
3d7758c1 193 fEtaProfile = new TProfile("fEtaProfile", "fEtaProfile;MC #eta;MC #eta - ESD #eta", 120, -3, 3);
0fc41645 194 fOutput->Add(fEtaProfile);
3d7758c1 195
4ef701a0 196 fMultAll = new TH1F("fMultAll", "fMultAll", 500, -0.5, 499.5);
0fc41645 197 fOutput->Add(fMultAll);
4ef701a0 198 fMultVtx = new TH1F("fMultVtx", "fMultVtx", 500, -0.5, 499.5);
0fc41645 199 fOutput->Add(fMultVtx);
4ef701a0 200 fMultTr = new TH1F("fMultTr", "fMultTr", 500, -0.5, 499.5);
0fc41645 201 fOutput->Add(fMultTr);
4ef701a0 202
3d7758c1 203 for (Int_t i=0; i<8; i++)
0fc41645 204 {
3d7758c1 205 fDeltaPhi[i] = new TH1F(Form("fDeltaPhi_%d", i), ";#Delta phi;Entries", 2000, -0.1, 0.1);
0fc41645 206 fOutput->Add(fDeltaPhi[i]);
207 }
50ec344d 208
209 fEventStats = new TH2F("fEventStats", "fEventStats;event type;status;count", 2, -0.5, 1.5, 4, -0.5, 3.5);
0fc41645 210 fOutput->Add(fEventStats);
50ec344d 211 fEventStats->GetXaxis()->SetBinLabel(1, "INEL");
212 fEventStats->GetXaxis()->SetBinLabel(2, "NSD");
213
214 fEventStats->GetYaxis()->SetBinLabel(1, "nothing");
215 fEventStats->GetYaxis()->SetBinLabel(2, "trg");
216 fEventStats->GetYaxis()->SetBinLabel(3, "vtx");
217 fEventStats->GetYaxis()->SetBinLabel(4, "trgvtx");
0f67a57c 218}
219
220void AlidNdEtaCorrectionTask::Exec(Option_t*)
221{
222 // process the event
223
224 // Check prerequisites
225 if (!fESD)
226 {
227 AliDebug(AliLog::kError, "ESD branch not available");
228 return;
229 }
230
3d7758c1 231 if (fOnlyPrimaries)
232 Printf("WARNING: Processing only primaries. For systematical studies only!");
233
0f67a57c 234 // trigger definition
0fc41645 235 Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), fTrigger);
54b096ef 236
3d7758c1 237 if (!eventTriggered)
238 Printf("No trigger");
239
54b096ef 240 // post the data already here
241 PostData(0, fOutput);
242
243 // MC info
244 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
245 if (!eventHandler) {
246 Printf("ERROR: Could not retrieve MC event handler");
247 return;
248 }
249
250 AliMCEvent* mcEvent = eventHandler->MCEvent();
251 if (!mcEvent) {
252 Printf("ERROR: Could not retrieve MC event");
253 return;
254 }
255
256 AliStack* stack = mcEvent->Stack();
257 if (!stack)
258 {
259 AliDebug(AliLog::kError, "Stack not available");
260 return;
261 }
262
263 AliHeader* header = mcEvent->Header();
264 if (!header)
265 {
266 AliDebug(AliLog::kError, "Header not available");
267 return;
268 }
269
270 // get process type; NB: this only works for Pythia
271 Int_t processType = AliPWG0Helper::GetPythiaEventProcessType(header);
272 AliDebug(AliLog::kDebug+1, Form("Found pythia process type %d", processType));
273
274 if (processType<0)
275 AliDebug(AliLog::kError, Form("Unknown Pythia process type %d.", processType));
276
277 // get the MC vertex
278 AliGenEventHeader* genHeader = header->GenEventHeader();
279 TArrayF vtxMC(3);
280 genHeader->PrimaryVertex(vtxMC);
281
282 // get the ESD vertex
283 const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
0f67a57c 284
770a1f1d 285 Bool_t eventVertex = kFALSE;
54b096ef 286 Double_t vtx[3];
287 if (vtxESD)
288 {
289 vtxESD->GetXYZ(vtx);
770a1f1d 290 eventVertex = kTRUE;
54b096ef 291
292 Double_t diff = vtxMC[2] - vtx[2];
293 if (vtxESD->GetZRes() > 0)
294 fVertexShiftNorm->Fill(diff / vtxESD->GetZRes());
295 }
296 else
297 Printf("No vertex found");
0f67a57c 298
50ec344d 299 // fill process type
300 Int_t biny = (Int_t) eventTriggered + 2 * (Int_t) eventVertex;
301 // INEL
302 fEventStats->Fill(0.0, biny);
303 // NDS
304 if (processType != 92 && processType != 93)
305 fEventStats->Fill(1.0, biny);
306
0f67a57c 307 // create list of (label, eta, pt) tuples
308 Int_t inputCount = 0;
309 Int_t* labelArr = 0;
3d7758c1 310 Int_t* labelArr2 = 0; // only for case of SPD
0f67a57c 311 Float_t* etaArr = 0;
312 Float_t* ptArr = 0;
4ef701a0 313 Float_t* deltaPhiArr = 0;
770a1f1d 314 if (fAnalysisMode == AliPWG0Helper::kSPD)
0f67a57c 315 {
316 // get tracklets
317 const AliMultiplicity* mult = fESD->GetMultiplicity();
318 if (!mult)
319 {
320 AliDebug(AliLog::kError, "AliMultiplicity not available");
321 return;
322 }
323
324 labelArr = new Int_t[mult->GetNumberOfTracklets()];
3d7758c1 325 labelArr2 = new Int_t[mult->GetNumberOfTracklets()];
0f67a57c 326 etaArr = new Float_t[mult->GetNumberOfTracklets()];
327 ptArr = new Float_t[mult->GetNumberOfTracklets()];
4ef701a0 328 deltaPhiArr = new Float_t[mult->GetNumberOfTracklets()];
0f67a57c 329
330 // get multiplicity from ITS tracklets
331 for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
332 {
333 //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
334
4ef701a0 335 Float_t deltaPhi = mult->GetDeltaPhi(i);
336 // prevent values to be shifted by 2 Pi()
337 if (deltaPhi < -TMath::Pi())
338 deltaPhi += TMath::Pi() * 2;
339 if (deltaPhi > TMath::Pi())
340 deltaPhi -= TMath::Pi() * 2;
341
342 if (TMath::Abs(deltaPhi) > 1)
343 printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
0f67a57c 344
3d7758c1 345 if (fOnlyPrimaries)
346 if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
347 continue;
348
0f67a57c 349 etaArr[inputCount] = mult->GetEta(i);
3d7758c1 350 labelArr[inputCount] = mult->GetLabel(i, 0);
351 labelArr2[inputCount] = mult->GetLabel(i, 1);
0f67a57c 352 ptArr[inputCount] = 0; // no pt for tracklets
4ef701a0 353 deltaPhiArr[inputCount] = deltaPhi;
0f67a57c 354 ++inputCount;
355 }
356 }
770a1f1d 357 else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
0f67a57c 358 {
359 if (!fEsdTrackCuts)
360 {
361 AliDebug(AliLog::kError, "fESDTrackCuts not available");
362 return;
363 }
364
365 // get multiplicity from ESD tracks
366 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
367 Int_t nGoodTracks = list->GetEntries();
368
54b096ef 369 Printf("Accepted %d tracks", nGoodTracks);
370
0f67a57c 371 labelArr = new Int_t[nGoodTracks];
3d7758c1 372 labelArr2 = new Int_t[nGoodTracks];
0f67a57c 373 etaArr = new Float_t[nGoodTracks];
374 ptArr = new Float_t[nGoodTracks];
4ef701a0 375 deltaPhiArr = new Float_t[nGoodTracks];
0f67a57c 376
377 // loop over esd tracks
378 for (Int_t i=0; i<nGoodTracks; i++)
379 {
380 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
381 if (!esdTrack)
382 {
383 AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
384 continue;
385 }
386
387 etaArr[inputCount] = esdTrack->Eta();
388 labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
3d7758c1 389 labelArr2[inputCount] = labelArr[inputCount]; // no second label for tracks
0f67a57c 390 ptArr[inputCount] = esdTrack->Pt();
4ef701a0 391 deltaPhiArr[inputCount] = 0; // no delta phi for tracks
0f67a57c 392 ++inputCount;
393 }
0f67a57c 394
54b096ef 395 delete list;
0f67a57c 396 }
54b096ef 397 else
0f67a57c 398 return;
0f67a57c 399
54b096ef 400 if (eventTriggered && eventVertex)
0f67a57c 401 {
54b096ef 402 fVertexCorrelation->Fill(vtxMC[2], vtx[2]);
403 fVertexProfile->Fill(vtxMC[2], vtxMC[2] - vtx[2]);
0f67a57c 404 }
405
0f67a57c 406
407 // loop over mc particles
408 Int_t nPrim = stack->GetNprimary();
4ef701a0 409 Int_t nAccepted = 0;
0f67a57c 410
411 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
412 {
413 //Printf("Getting particle %d", iMc);
414 TParticle* particle = stack->Particle(iMc);
415
416 if (!particle)
417 continue;
418
419 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
420 continue;
421
422 if (SignOK(particle->GetPDG()) == kFALSE)
423 continue;
424
425 Float_t eta = particle->Eta();
426 Float_t pt = particle->Pt();
427
428 fdNdEtaCorrection->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
429
96fcc8a7 430 if (fdNdEtaCorrectionProcessType[0])
431 {
432 // non diffractive
433 if (processType!=92 && processType!=93 && processType!=94)
434 fdNdEtaCorrectionProcessType[0]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
435
436 // single diffractive
437 if (processType==92 || processType==93)
438 fdNdEtaCorrectionProcessType[1]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
439
440 // double diffractive
441 if (processType==94)
442 fdNdEtaCorrectionProcessType[2]->FillMCParticle(vtxMC[2], eta, pt, eventTriggered, eventVertex, processType);
443 }
444
0f67a57c 445 if (eventTriggered)
446 if (eventVertex)
447 fdNdEtaAnalysisMC->FillTrack(vtxMC[2], eta, pt);
4ef701a0 448
0fc41645 449 // TODO this value might be needed lower for the SPD study (only used for control histograms anyway)
4ef701a0 450 if (TMath::Abs(eta) < 1 && pt > 0.2)
451 nAccepted++;
452 }
453
454 fMultAll->Fill(nAccepted);
455 if (eventTriggered) {
456 fMultTr->Fill(nAccepted);
457 if (eventVertex)
458 fMultVtx->Fill(nAccepted);
0f67a57c 459 }
460
3d7758c1 461 Int_t processed = 0;
462
0f67a57c 463 for (Int_t i=0; i<inputCount; ++i)
464 {
465 Int_t label = labelArr[i];
3d7758c1 466 Int_t label2 = labelArr2[i];
0f67a57c 467
468 if (label < 0)
469 {
0fc41645 470 Printf("WARNING: cannot find corresponding mc particle for track(let) %d with label %d.", i, label);
0f67a57c 471 continue;
472 }
473
474 TParticle* particle = stack->Particle(label);
475 if (!particle)
476 {
477 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
478 continue;
479 }
3d7758c1 480
481 // find particle that is filled in the correction map
482 // this should be the particle that has been reconstructed
483 // for TPC this is clear and always identified by the track's label
484 // for SPD the labels might not be identical. In this case the particle closest to the beam line that is a primary is taken:
485 // L1 L2 (P = primary, S = secondary)
486 // P P' : --> P
487 // P S : --> P
488 // S P : --> P
489 // S S' : --> S
490
491 if (label != label2 && label2 >= 0)
492 {
493 if (stack->IsPhysicalPrimary(label) == kFALSE && stack->IsPhysicalPrimary(label2))
494 {
495 particle = stack->Particle(label2);
496 if (!particle)
497 {
498 AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label2));
499 continue;
500 }
501 }
502 }
0f67a57c 503
504 if (eventTriggered && eventVertex)
505 {
3d7758c1 506 if (SignOK(particle->GetPDG()) == kFALSE)
507 continue;
508
509 processed++;
510
511 Bool_t firstIsPrim = stack->IsPhysicalPrimary(label);
512 // in case of primary the MC values are filled, otherwise (background) the reconstructed values
513 if (label == label2 && firstIsPrim)
514 {
515 fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
516 }
517 else
518 {
519 fdNdEtaCorrection->FillTrackedParticle(vtxMC[2], etaArr[i], ptArr[i]);
3d7758c1 520 }
745d6088 521
522 fEtaProfile->Fill(particle->Eta(), particle->Eta() - etaArr[i]);
0f67a57c 523 fdNdEtaAnalysisESD->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
3d7758c1 524
525 fEtaCorrelation->Fill(etaArr[i], particle->Eta());
526
0f67a57c 527 if (particle->Pt() > 0.1 && particle->Pt() < 0.2)
528 {
529 fPIDTracks->Fill(particle->GetPdgCode());
530 }
96fcc8a7 531
532 if (fdNdEtaCorrectionProcessType[0])
533 {
534 // non diffractive
535 if (processType!=92 && processType!=93 && processType!=94)
536 fdNdEtaCorrectionProcessType[0]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
537
538 // single diffractive
539 if (processType==92 || processType==93)
540 fdNdEtaCorrectionProcessType[1]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
541
542 // double diffractive
543 if (processType==94)
544 fdNdEtaCorrectionProcessType[2]->FillTrackedParticle(vtxMC[2], particle->Eta(), particle->Pt());
545 }
4ef701a0 546
745d6088 547 // control histograms
3d7758c1 548 Int_t hist = -1;
549 if (label == label2)
4ef701a0 550 {
3d7758c1 551 if (firstIsPrim)
552 {
553 hist = 0;
554 }
555 else
556 hist = 1;
557 }
558 else if (label2 >= 0)
559 {
560 Bool_t secondIsPrim = stack->IsPhysicalPrimary(label2);
561 if (firstIsPrim && secondIsPrim)
562 {
563 hist = 2;
564 }
565 else if (firstIsPrim && !secondIsPrim)
566 {
567 hist = 3;
568
569 // check if secondary is caused by the primary or it is a fake combination
570 //Printf("PS case --> Label 1 is %d, label 2 is %d", label, label2);
571 Int_t mother = label2;
572 while (stack->Particle(mother) && stack->Particle(mother)->GetMother(0) >= 0)
573 {
574 //Printf(" %d created by %d, %d", mother, stack->Particle(mother)->GetMother(0), stack->Particle(mother)->GetMother(1));
575 if (stack->Particle(mother)->GetMother(0) == label)
576 {
577 /*Printf("The untraceable decay was:");
578 Printf(" from:");
579 particle->Print();
580 Printf(" to (status code %d):", stack->Particle(mother)->GetStatusCode());
581 stack->Particle(mother)->Print();*/
582 hist = 4;
583 }
584 mother = stack->Particle(mother)->GetMother(0);
585 }
586 }
587 else if (!firstIsPrim && secondIsPrim)
588 {
589 hist = 5;
590 }
591 else if (!firstIsPrim && !secondIsPrim)
592 {
593 hist = 6;
594 }
595
4ef701a0 596 }
597 else
3d7758c1 598 hist = 7;
599
600 fDeltaPhi[hist]->Fill(deltaPhiArr[i]);
0f67a57c 601 }
602 }
603
3d7758c1 604 if (processed < inputCount)
605 Printf("Only %d out of %d track(let)s were used", processed, inputCount);
606
0f67a57c 607 if (eventTriggered && eventVertex)
96fcc8a7 608 {
0f67a57c 609 fdNdEtaAnalysisMC->FillEvent(vtxMC[2], inputCount);
96fcc8a7 610 fdNdEtaAnalysisESD->FillEvent(vtxMC[2], inputCount);
611 }
0f67a57c 612
96fcc8a7 613 // stuff regarding the vertex reco correction and trigger bias correction
0f67a57c 614 fdNdEtaCorrection->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
96fcc8a7 615
616 if (fdNdEtaCorrectionProcessType[0])
617 {
618 // non diffractive
619 if (processType!=92 && processType!=93 && processType!=94)
620 fdNdEtaCorrectionProcessType[0]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
621
622 // single diffractive
623 if (processType==92 || processType==93)
624 fdNdEtaCorrectionProcessType[1]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
625
626 // double diffractive
627 if (processType==94)
628 fdNdEtaCorrectionProcessType[2]->FillEvent(vtxMC[2], inputCount, eventTriggered, eventVertex, processType);
0f67a57c 629 }
630
631 delete[] etaArr;
632 delete[] labelArr;
3d7758c1 633 delete[] labelArr2;
0f67a57c 634 delete[] ptArr;
4ef701a0 635 delete[] deltaPhiArr;
96fcc8a7 636
637 // fills the fSigmaVertex histogram (systematic study)
638 if (fSigmaVertexTracks)
639 {
640 // save the old value
641 Float_t oldSigmaVertex = fEsdTrackCuts->GetMinNsigmaToVertex();
642
643 // set to maximum
644 fEsdTrackCuts->SetMinNsigmaToVertex(6);
645
646 TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
647 Int_t nGoodTracks = list->GetEntries();
648
649 // loop over esd tracks
650 for (Int_t i=0; i<nGoodTracks; i++)
651 {
652 AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
653 if (!esdTrack)
654 {
655 AliError(Form("ERROR: Could not retrieve track %d.", i));
656 continue;
657 }
658
659 Float_t sigma2Vertex = fEsdTrackCuts->GetSigmaToVertex(esdTrack);
660
661 for (Double_t nSigma = 0.1; nSigma < 6.05; nSigma += 0.1)
662 {
663 if (sigma2Vertex < nSigma)
664 {
665 fSigmaVertexTracks->Fill(nSigma);
666
667 Int_t label = TMath::Abs(esdTrack->GetLabel());
668 TParticle* particle = stack->Particle(label);
669 if (!particle || label >= nPrim)
670 continue;
671
672 if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim))
673 fSigmaVertexPrim->Fill(nSigma);
674 }
675 }
676 }
677
678 delete list;
679 list = 0;
680
681 // set back the old value
682 fEsdTrackCuts->SetMinNsigmaToVertex(oldSigmaVertex);
683 }
0f67a57c 684}
685
686void AlidNdEtaCorrectionTask::Terminate(Option_t *)
687{
688 // The Terminate() function is the last function to be called during
689 // a query. It always runs on the client, it can be used to present
690 // the results graphically or save the results to file.
691
692 fOutput = dynamic_cast<TList*> (GetOutputData(0));
693 if (!fOutput) {
694 Printf("ERROR: fOutput not available");
695 return;
696 }
697
698 fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
699 fdNdEtaAnalysisMC = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaMC"));
700 fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaESD"));
701 if (!fdNdEtaCorrection || !fdNdEtaAnalysisMC || !fdNdEtaAnalysisESD)
702 {
703 AliDebug(AliLog::kError, "Could not read object from output list");
704 return;
705 }
706
707 fdNdEtaCorrection->Finish();
708
709 TString fileName;
710 fileName.Form("correction_map%s.root", fOption.Data());
711 TFile* fout = new TFile(fileName, "RECREATE");
712
713 if (fEsdTrackCuts)
714 fEsdTrackCuts->SaveHistograms("esd_track_cuts");
715 fdNdEtaCorrection->SaveHistograms();
716 fdNdEtaAnalysisMC->SaveHistograms();
717 fdNdEtaAnalysisESD->SaveHistograms();
718
96fcc8a7 719 fdNdEtaCorrectionProcessType[0] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_ND"));
720 fdNdEtaCorrectionProcessType[1] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_SD"));
721 fdNdEtaCorrectionProcessType[2] = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction_DD"));
722 for (Int_t i=0; i<3; ++i)
723 if (fdNdEtaCorrectionProcessType[i])
724 fdNdEtaCorrectionProcessType[i]->SaveHistograms();
725
726 fSigmaVertexTracks = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexTracks"));
727 if (fSigmaVertexTracks)
728 fSigmaVertexTracks->Write();
729 fSigmaVertexPrim = dynamic_cast<TH1F*> (fOutput->FindObject("fSigmaVertexPrim"));
730 if (fSigmaVertexPrim)
731 fSigmaVertexPrim->Write();
732
0fc41645 733 fVertexCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fVertexCorrelation"));
54b096ef 734 if (fVertexCorrelation)
735 fVertexCorrelation->Write();
0fc41645 736 fVertexProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fVertexProfile"));
54b096ef 737 if (fVertexProfile)
738 fVertexProfile->Write();
0fc41645 739 fVertexShiftNorm = dynamic_cast<TH1F*> (fOutput->FindObject("fVertexShiftNorm"));
54b096ef 740 if (fVertexShiftNorm)
741 fVertexShiftNorm->Write();
3d7758c1 742
0fc41645 743 fEtaCorrelation = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaCorrelation"));
3d7758c1 744 if (fEtaCorrelation)
745 fEtaCorrelation->Write();
0fc41645 746 fEtaProfile = dynamic_cast<TProfile*> (fOutput->FindObject("fEtaProfile"));
3d7758c1 747 if (fEtaProfile)
748 fEtaProfile->Write();
749
0fc41645 750 fMultAll = dynamic_cast<TH1F*> (fOutput->FindObject("fMultAll"));
4ef701a0 751 if (fMultAll)
752 fMultAll->Write();
753
0fc41645 754 fMultTr = dynamic_cast<TH1F*> (fOutput->FindObject("fMultTr"));
4ef701a0 755 if (fMultTr)
756 fMultTr->Write();
0fc41645 757
758 fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
4ef701a0 759 if (fMultVtx)
760 fMultVtx->Write();
761
3d7758c1 762 for (Int_t i=0; i<8; ++i)
0fc41645 763 {
764 fDeltaPhi[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("fDeltaPhi_%d", i)));
4ef701a0 765 if (fDeltaPhi[i])
766 fDeltaPhi[i]->Write();
0fc41645 767 }
54b096ef 768
0fc41645 769 fEventStats = dynamic_cast<TH2F*> (fOutput->FindObject("fEventStats"));
50ec344d 770 if (fEventStats)
771 fEventStats->Write();
772
0f67a57c 773 fout->Write();
774 fout->Close();
775
54b096ef 776 //fdNdEtaCorrection->DrawHistograms();
0f67a57c 777
96fcc8a7 778 Printf("Writting result to %s", fileName.Data());
779
0f67a57c 780 if (fPIDParticles && fPIDTracks)
781 {
782 new TCanvas("pidcanvas", "pidcanvas", 500, 500);
783
784 fPIDParticles->Draw();
785 fPIDTracks->SetLineColor(2);
786 fPIDTracks->Draw("SAME");
787
788 TDatabasePDG* pdgDB = new TDatabasePDG;
789
790 for (Int_t i=0; i <= fPIDParticles->GetNbinsX()+1; ++i)
791 if (fPIDParticles->GetBinContent(i) > 0)
792 printf("PDG = %d (%s): generated: %d, reconstructed: %d, ratio: %f\n", (Int_t) fPIDParticles->GetBinCenter(i), pdgDB->GetParticle((Int_t) fPIDParticles->GetBinCenter(i))->GetName(), (Int_t) fPIDParticles->GetBinContent(i), (Int_t) fPIDTracks->GetBinContent(i), ((fPIDTracks->GetBinContent(i) > 0) ? fPIDParticles->GetBinContent(i) / fPIDTracks->GetBinContent(i) : -1));
793
794 delete pdgDB;
795 pdgDB = 0;
796 }
0f67a57c 797}
798
799Bool_t AlidNdEtaCorrectionTask::SignOK(TParticlePDG* particle)
800{
801 // returns if a particle with this sign should be counted
802 // this is determined by the value of fSignMode, which should have the same sign
803 // as the charge
804 // if fSignMode is 0 all particles are counted
805
806 if (fSignMode == 0)
807 return kTRUE;
808
809 if (!particle)
810 {
811 printf("WARNING: not counting a particle that does not have a pdg particle\n");
812 return kFALSE;
813 }
814
815 Double_t charge = particle->Charge();
816
817 if (fSignMode > 0)
818 if (charge < 0)
819 return kFALSE;
820
821 if (fSignMode < 0)
822 if (charge > 0)
823 return kFALSE;
824
825 return kTRUE;
826}
96fcc8a7 827