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