573ef5a4d0440c18e6180b5a0e66d7e85f8f795f
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaTask.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaTask.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>
11 #include <TH1F.h>
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
32 #include "esdTrackCuts/AliESDtrackCuts.h"
33 #include "AliPWG0Helper.h"
34 #include "AliCorrection.h"
35 #include "AliCorrectionMatrix3D.h"
36 #include "dNdEta/dNdEtaAnalysis.h"
37
38 ClassImp(AlidNdEtaTask)
39
40 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
41   AliAnalysisTask("AlidNdEtaTask", ""),
42   fESD(0),
43   fOutput(0),
44   fOption(opt),
45   fAnalysisMode(AliPWG0Helper::kTPC),
46   fReadMC(kFALSE),
47   fUseMCVertex(kFALSE),
48   fUseMCKine(kFALSE),
49   fEsdTrackCuts(0),
50   fdNdEtaAnalysisESD(0),
51   fMult(0),
52   fMultVtx(0),
53   fEvents(0),
54   fVertexResolution(0),
55   fdNdEtaAnalysis(0),
56   fdNdEtaAnalysisTr(0),
57   fdNdEtaAnalysisTrVtx(0),
58   fdNdEtaAnalysisTracks(0),
59   fVertex(0),
60   fPartPt(0),
61   fDeltaPhi(0)
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());
70 }
71
72 AlidNdEtaTask::~AlidNdEtaTask()
73 {
74   //
75   // Destructor
76   //
77
78   // histograms are in the output list and deleted when the output
79   // list is deleted by the TSelector dtor
80
81   if (fOutput) {
82     delete fOutput;
83     fOutput = 0;
84   }
85 }
86
87 //________________________________________________________________________
88 void AlidNdEtaTask::ConnectInputData(Option_t *)
89 {
90   // Connect ESD
91   // Called once
92
93   Printf("AlidNdEtaTask::ConnectInputData called");
94
95   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
96   if (!tree) {
97     Printf("ERROR: Could not read tree from input slot 0");
98   } else {
99     // Disable all branches and enable only the needed ones
100     //tree->SetBranchStatus("*", 0);
101
102     tree->SetBranchStatus("fTriggerMask", 1);
103     tree->SetBranchStatus("fSPDVertex*", 1);
104     // PrimaryVertex also needed
105
106     if (fAnalysisMode == AliPWG0Helper::kSPD) {
107       tree->SetBranchStatus("fSPDMult*", 1);
108       //AliPWG0Helper::SetBranchStatusRecursive(tree, "AliMultiplicity", kTRUE);
109     }
110     
111     if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
112       AliESDtrackCuts::EnableNeededBranches(tree);
113       tree->SetBranchStatus("fTracks.fLabel", 1);
114       //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
115     }
116
117     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
118
119     if (!esdH) {
120       Printf("ERROR: Could not get ESDInputHandler");
121     } else
122       fESD = esdH->GetEvent();
123   }
124
125   // disable info messages of AliMCEvent (per event)
126   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
127 }
128
129 void AlidNdEtaTask::CreateOutputObjects()
130 {
131   // create result objects and add to output list
132
133   Printf("AlidNdEtaTask::CreateOutputObjects");
134
135   fOutput = new TList;
136   fOutput->SetOwner();
137
138   fdNdEtaAnalysisESD = new dNdEtaAnalysis("fdNdEtaAnalysisESD", "fdNdEtaAnalysisESD", fAnalysisMode);
139   fOutput->Add(fdNdEtaAnalysisESD);
140
141   fMult = new TH1F("fMult", "fMult;Ntracks;Count", 201, -0.5, 200.5);
142   fOutput->Add(fMult);
143
144   fMultVtx = new TH1F("fMultVtx", "fMultVtx;Ntracks;Count", 201, -0.5, 200.5);
145   fOutput->Add(fMultVtx);
146
147   for (Int_t i=0; i<3; ++i)
148   {
149     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
150     fPartEta[i]->Sumw2();
151     fOutput->Add(fPartEta[i]);
152   }
153
154   fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
155   fOutput->Add(fEvents);
156
157   fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
158   fOutput->Add(fVertexResolution);
159
160   if (fReadMC)
161   {
162     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
163     fOutput->Add(fdNdEtaAnalysis);
164
165     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
166     fOutput->Add(fdNdEtaAnalysisTr);
167
168     fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx", fAnalysisMode);
169     fOutput->Add(fdNdEtaAnalysisTrVtx);
170
171     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
172     fOutput->Add(fdNdEtaAnalysisTracks);
173
174     fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
175     fOutput->Add(fVertex);
176
177     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
178     fPartPt->Sumw2();
179     fOutput->Add(fPartPt);
180   }
181
182   if (fAnalysisMode == AliPWG0Helper::kSPD)
183     fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
184 }
185
186 void AlidNdEtaTask::Exec(Option_t*)
187 {
188   // process the event
189
190   // Check prerequisites
191   if (!fESD)
192   {
193     AliDebug(AliLog::kError, "ESD branch not available");
194     return;
195   }
196
197   // trigger definition
198   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
199
200   // get the ESD vertex
201   const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
202   
203   Double_t vtx[3];
204   if (vtxESD) {
205     vtxESD->GetXYZ(vtx);
206   }
207   else
208     Printf("No vertex");
209   
210   // fill z vertex resolution
211   if (vtxESD)
212     fVertexResolution->Fill(vtxESD->GetZRes());
213
214   // post the data already here
215   PostData(0, fOutput);
216
217   // needed for syst. studies
218   AliStack* stack = 0;
219   
220   if (fUseMCVertex || fUseMCKine) {
221     if (!fReadMC) {
222       Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
223       return;
224     }
225
226     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
227     if (!eventHandler) {
228       Printf("ERROR: Could not retrieve MC event handler");
229       return;
230     }
231
232     AliMCEvent* mcEvent = eventHandler->MCEvent();
233     if (!mcEvent) {
234       Printf("ERROR: Could not retrieve MC event");
235       return;
236     }
237
238     AliHeader* header = mcEvent->Header();
239     if (!header)
240     {
241       AliDebug(AliLog::kError, "Header not available");
242       return;
243     }
244
245     if (fUseMCVertex)
246     {
247       Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
248       // get the MC vertex
249       AliGenEventHeader* genHeader = header->GenEventHeader();
250       TArrayF vtxMC(3);
251       genHeader->PrimaryVertex(vtxMC);
252
253       vtx[2] = vtxMC[2];
254     }
255
256     if (fUseMCKine)
257     {
258
259       stack = mcEvent->Stack();
260       if (!stack)
261       {
262         AliDebug(AliLog::kError, "Stack not available");
263         return;
264       }
265     }
266   }
267
268   // create list of (label, eta, pt) tuples
269   Int_t inputCount = 0;
270   Int_t* labelArr = 0;
271   Float_t* etaArr = 0;
272   Float_t* ptArr = 0;
273   if (fAnalysisMode == AliPWG0Helper::kSPD)
274   {
275     // get tracklets
276     const AliMultiplicity* mult = fESD->GetMultiplicity();
277     if (!mult)
278     {
279       AliDebug(AliLog::kError, "AliMultiplicity not available");
280       return;
281     }
282
283     labelArr = new Int_t[mult->GetNumberOfTracklets()];
284     etaArr = new Float_t[mult->GetNumberOfTracklets()];
285     ptArr = new Float_t[mult->GetNumberOfTracklets()];
286
287     if (fUseMCKine && stack)
288       Printf("Processing only primaries (MC information used). This is for systematical checks only.");
289
290     // get multiplicity from ITS tracklets
291     for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
292     {
293       //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
294
295       if (fUseMCKine && stack)
296         if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
297           continue;
298       
299       Float_t deltaPhi = mult->GetDeltaPhi(i);
300       // prevent values to be shifted by 2 Pi()
301       if (deltaPhi < -TMath::Pi())
302         deltaPhi += TMath::Pi() * 2;
303       if (deltaPhi > TMath::Pi())
304         deltaPhi -= TMath::Pi() * 2;
305
306       if (TMath::Abs(deltaPhi) > 1)
307         printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
308
309       fDeltaPhi->Fill(deltaPhi);
310
311       etaArr[inputCount] = mult->GetEta(i);
312       labelArr[inputCount] = mult->GetLabel(i, 0);
313       ptArr[inputCount] = 0; // no pt for tracklets
314       ++inputCount;
315     }
316
317     Printf("Accepted %d tracklets", inputCount);
318   }
319   else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
320   {
321     if (!fEsdTrackCuts)
322     {
323       AliDebug(AliLog::kError, "fESDTrackCuts not available");
324       return;
325     }
326
327     // get multiplicity from ESD tracks
328     TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
329     Int_t nGoodTracks = list->GetEntries();
330
331     labelArr = new Int_t[nGoodTracks];
332     etaArr = new Float_t[nGoodTracks];
333     ptArr = new Float_t[nGoodTracks];
334
335     // loop over esd tracks
336     for (Int_t i=0; i<nGoodTracks; i++)
337     {
338       AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
339       if (!esdTrack)
340       {
341         AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
342         continue;
343       }
344
345       etaArr[inputCount] = esdTrack->Eta();
346       labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
347       ptArr[inputCount] = esdTrack->Pt();
348       ++inputCount;
349     }
350
351     Printf("Accepted %d tracks", nGoodTracks);
352
353     delete list;
354   }
355   else
356     return;
357
358   // Processing of ESD information (always)
359   if (eventTriggered)
360   {
361     // control hist
362     fMult->Fill(inputCount);
363     fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
364
365     if (vtxESD)
366     {
367       // control hist
368       fMultVtx->Fill(inputCount);
369
370       for (Int_t i=0; i<inputCount; ++i)
371       {
372         Float_t eta = etaArr[i];
373         Float_t pt  = ptArr[i];
374
375         fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
376
377         if (TMath::Abs(vtx[2]) < 20)
378         {
379           fPartEta[0]->Fill(eta);
380
381           if (vtx[2] < 0)
382             fPartEta[1]->Fill(eta);
383           else
384             fPartEta[2]->Fill(eta);
385         }
386       }
387
388       // for event count per vertex
389       fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
390
391       // control hist
392       fEvents->Fill(vtx[2]);
393     }
394   }
395
396   if (fReadMC)   // Processing of MC information (optional)
397   {
398     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
399     if (!eventHandler) {
400       Printf("ERROR: Could not retrieve MC event handler");
401       return;
402     }
403
404     AliMCEvent* mcEvent = eventHandler->MCEvent();
405     if (!mcEvent) {
406       Printf("ERROR: Could not retrieve MC event");
407       return;
408     }
409
410     AliStack* stack = mcEvent->Stack();
411     if (!stack)
412     {
413       AliDebug(AliLog::kError, "Stack not available");
414       return;
415     }
416
417     AliHeader* header = mcEvent->Header();
418     if (!header)
419     {
420       AliDebug(AliLog::kError, "Header not available");
421       return;
422     }
423
424     // get the MC vertex
425     AliGenEventHeader* genHeader = header->GenEventHeader();
426     TArrayF vtxMC(3);
427     genHeader->PrimaryVertex(vtxMC);
428
429     // loop over mc particles
430     Int_t nPrim  = stack->GetNprimary();
431
432     Int_t nAcceptedParticles = 0;
433
434     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
435     {
436       //Printf("Getting particle %d", iMc);
437       TParticle* particle = stack->Particle(iMc);
438
439       if (!particle)
440         continue;
441
442       if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
443         continue;
444
445       AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
446       Float_t eta = particle->Eta();
447       Float_t pt = particle->Pt();
448
449        // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
450       if (TMath::Abs(eta) < 1.5 && pt > 0.3)
451         nAcceptedParticles++;
452
453       fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
454       fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
455
456       if (eventTriggered)
457       {
458         fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
459         if (vtxESD)
460           fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
461       }
462
463       if (TMath::Abs(eta) < 0.8)
464         fPartPt->Fill(pt);
465     }
466
467     fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
468     if (eventTriggered)
469     {
470       fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
471       if (vtxESD)
472         fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
473     }
474
475     if (eventTriggered && vtxESD)
476     {
477       // from tracks is only done for triggered and vertex reconstructed events
478
479       for (Int_t i=0; i<inputCount; ++i)
480       {
481         Int_t label = labelArr[i];
482
483         if (label < 0)
484           continue;
485
486         //Printf("Getting particle of track %d", label);
487         TParticle* particle = stack->Particle(label);
488
489         if (!particle)
490         {
491           AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
492           continue;
493         }
494
495         fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
496       } // end of track loop
497
498       // for event count per vertex
499       fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
500     }
501   }
502
503   delete[] etaArr;
504   delete[] labelArr;
505   delete[] ptArr;
506 }
507
508 void AlidNdEtaTask::Terminate(Option_t *)
509 {
510   // The Terminate() function is the last function to be called during
511   // a query. It always runs on the client, it can be used to present
512   // the results graphically or save the results to file.
513
514   fOutput = dynamic_cast<TList*> (GetOutputData(0));
515   if (!fOutput) {
516     Printf("ERROR: fOutput not available");
517     return;
518   }
519
520   fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
521   fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
522   fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
523   fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
524   fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
525   for (Int_t i=0; i<3; ++i)
526     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
527
528   if (!fdNdEtaAnalysisESD)
529   {
530     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
531     return;
532   }
533
534   if (fMult && fMultVtx)
535   {
536     new TCanvas;
537     fMult->Draw();
538     fMultVtx->SetLineColor(2);
539     fMultVtx->Draw("SAME");
540   }
541
542   if (fPartEta[0])
543   {
544     Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
545     Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
546
547     Printf("%d events with vertex used", events1 + events2);
548
549     if (events1 > 0 && events2 > 0)
550     {
551     fPartEta[0]->Scale(1.0 / (events1 + events2));
552     fPartEta[1]->Scale(1.0 / events1);
553     fPartEta[2]->Scale(1.0 / events2);
554
555     for (Int_t i=0; i<3; ++i)
556       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
557
558     new TCanvas("control", "control", 500, 500);
559     for (Int_t i=0; i<3; ++i)
560     {
561       fPartEta[i]->SetLineColor(i+1);
562       fPartEta[i]->Draw((i==0) ? "" : "SAME");
563     }
564     }
565   }
566
567   if (fEvents)
568   {
569     new TCanvas("control3", "control3", 500, 500);
570     fEvents->Draw();
571   }
572
573   TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
574
575   if (fdNdEtaAnalysisESD)
576     fdNdEtaAnalysisESD->SaveHistograms();
577
578   if (fEsdTrackCuts)
579     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
580
581   if (fMult)
582     fMult->Write();
583
584   if (fMultVtx)
585     fMultVtx->Write();
586
587   for (Int_t i=0; i<3; ++i)
588     if (fPartEta[i])
589       fPartEta[i]->Write();
590
591   if (fEvents)
592     fEvents->Write();
593
594   if (fVertexResolution)
595     fVertexResolution->Write();
596
597   if (fDeltaPhi)
598     fDeltaPhi->Write();
599
600   fout->Write();
601   fout->Close();
602
603   Printf("Writting result to analysis_esd_raw.root");
604
605   if (fReadMC)
606   {
607     fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
608     fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
609     fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
610     fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
611     fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
612
613     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
614     {
615       AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
616       return;
617     }
618
619     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
620     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
621     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
622     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
623
624     Int_t events = (Int_t) fdNdEtaAnalysis->GetData()->GetEventCorrection()->GetMeasuredHistogram()->Integral();
625     fPartPt->Scale(1.0/events);
626     fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
627
628     TFile* fout = new TFile("analysis_mc.root","RECREATE");
629
630     fdNdEtaAnalysis->SaveHistograms();
631     fdNdEtaAnalysisTr->SaveHistograms();
632     fdNdEtaAnalysisTrVtx->SaveHistograms();
633     fdNdEtaAnalysisTracks->SaveHistograms();
634     fPartPt->Write();
635
636     fout->Write();
637     fout->Close();
638
639     Printf("Writting result to analysis_mc.root");
640
641     if (fPartPt)
642     {
643       new TCanvas("control2", "control2", 500, 500);
644       fPartPt->Draw();
645     }
646   }
647 }