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