fd11abfcb2f3c166514d5b6547fa7c9253925f22
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisTaskNetParticle.cxx
1 //-*- Mode: C++ -*-
2
3 #include "TFile.h"
4 #include "TChain.h"
5 #include "TTree.h"
6 #include "TH1F.h"
7 #include "TH2F.h"
8 #include "TH3F.h"
9 #include "TCanvas.h"
10 #include "TStopwatch.h"
11 #include "TMath.h"
12 #include "THashList.h"
13
14 #include "AliAnalysisTask.h"
15 #include "AliAnalysisManager.h"
16 #include "AliTracker.h" 
17 #include "AliESDEvent.h"
18 #include "AliESDInputHandler.h"
19 #include "AliESDpid.h"
20 #include "AliStack.h"
21 #include "AliMCEvent.h"
22 #include "AliMCEventHandler.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliKineTrackCuts.h"
25 #include "AliMCParticle.h"
26 #include "AliESDVZERO.h"
27 #include "AliAnalysisTaskNetParticle.h"
28 #include "AliGenEventHeader.h"
29 #include "AliCentrality.h"
30 #include "AliAODEvent.h"
31 #include "AliAODInputHandler.h"
32
33 using namespace std;
34
35 // Task for NetParticle checks
36 // Author: Jochen Thaeder <jochen@thaeder.de>
37
38 ClassImp(AliAnalysisTaskNetParticle)
39
40 /*
41  * ---------------------------------------------------------------------------------
42  *                            Constructor / Destructor
43  * ---------------------------------------------------------------------------------
44  */
45
46 //________________________________________________________________________
47 AliAnalysisTaskNetParticle::AliAnalysisTaskNetParticle(const char *name) :
48   AliAnalysisTaskSE(name),
49   fHelper(new AliAnalysisNetParticleHelper),
50   fEffCont(NULL),
51   fDCA(NULL),
52   fDist(NULL),
53
54   fOutList(NULL),
55   fOutListEff(NULL),
56   fOutListCont(NULL),
57   fOutListDCA(NULL),
58   fOutListQA(NULL),
59
60   fESD(NULL), 
61   fESDHandler(NULL),
62
63   fESDTrackCutsBase(new AliESDtrackCuts),
64   fESDTrackCuts(NULL),
65   fESDTrackCutsBkg(NULL),
66   fESDTrackCutsEff(NULL),
67
68   fAOD(NULL), 
69   fAODHandler(NULL),
70
71   fIsMC(kFALSE),
72   fIsAOD(kFALSE),
73   fESDTrackCutMode(0),
74   fModeEffCreation(0),
75   fModeDCACreation(0),
76   fModeDistCreation(0),
77
78   fMCEvent(NULL),
79   fMCStack(NULL),
80
81   fHnQA(NULL),
82   fUseQATHnSparse(kFALSE),
83
84   fEtaMax(0.9),
85   fPtRange(new Float_t[2]),
86   fPtRangeEff(new Float_t[2]),
87
88   fAODtrackCutBit(1024){
89   // Constructor   
90
91   AliLog::SetClassDebugLevel("AliAnalysisTaskNetParticle",10);
92
93   fPtRange[0] = 0.4;
94   fPtRange[1] = 0.8;
95   fPtRangeEff[0] = 0.2;
96   fPtRangeEff[1] = 1.6;
97
98   // -- Output slots
99   // -------------------------------------------------
100   DefineOutput(1, TList::Class());
101   DefineOutput(2, TList::Class());
102   DefineOutput(3, TList::Class());
103   DefineOutput(4, TList::Class());
104   DefineOutput(5, TList::Class());
105 }
106
107 //________________________________________________________________________
108 AliAnalysisTaskNetParticle::~AliAnalysisTaskNetParticle() {
109   // Destructor
110
111   if (fPtRange)          delete[] fPtRange;
112   if (fPtRangeEff)       delete[] fPtRangeEff;
113
114   if (fESDTrackCutsBase) delete fESDTrackCutsBase;
115   if (fESDTrackCuts)     delete fESDTrackCuts;
116   if (fESDTrackCutsBkg)  delete fESDTrackCutsBkg;
117   if (fESDTrackCutsEff)  delete fESDTrackCutsEff;
118
119   if (fHelper)           delete fHelper;
120   if (fEffCont)          delete fEffCont;
121   if (fDCA)              delete fDCA;
122   if (fDist)             delete fDist;
123 }
124
125 /*
126  * ---------------------------------------------------------------------------------
127  *                                 Public Methods
128  * ---------------------------------------------------------------------------------
129  */
130
131 //________________________________________________________________________
132 void AliAnalysisTaskNetParticle::UserCreateOutputObjects() {
133   // Create histograms
134
135   Bool_t oldStatus = TH1::AddDirectoryStatus();
136   TH1::AddDirectory(kFALSE);
137
138   fOutList = new TList;
139   fOutList->SetName(GetName()) ;
140   fOutList->SetOwner(kTRUE);
141  
142   fOutListEff = new TList;
143   fOutListEff->SetName(Form("%s_eff",GetName()));
144   fOutListEff->SetOwner(kTRUE) ;
145
146   fOutListCont = new TList;
147   fOutListCont->SetName(Form("%s_cont",GetName()));
148   fOutListCont->SetOwner(kTRUE) ;
149
150   fOutListDCA = new TList;
151   fOutListDCA->SetName(Form("%s_dca",GetName()));
152   fOutListDCA->SetOwner(kTRUE) ;
153  
154   fOutListQA = new TList;
155   fOutListQA->SetName(Form("%s_qa",GetName()));
156   fOutListQA->SetOwner(kTRUE) ;
157  
158   // ------------------------------------------------------------------
159   // -- Get event / trigger statistics histograms
160   // ------------------------------------------------------------------
161
162   fOutList->Add(fHelper->GetHEventStat0());
163   fOutList->Add(fHelper->GetHEventStat1());
164   fOutList->Add(fHelper->GetHTriggerStat());
165   fOutList->Add(fHelper->GetHCentralityStat());
166
167   // ------------------------------------------------------------------
168   // -- Add histograms from distribution class
169   // ------------------------------------------------------------------
170   if (fModeDistCreation == 1)
171     fDist->CreateHistograms(fOutList);
172
173   // ------------------------------------------------------------------
174   // -- Add histograms from efficiency/contamination class
175   // ------------------------------------------------------------------
176   if ((fIsAOD||fIsMC) && fModeEffCreation == 1) {
177     fOutListEff->Add(fEffCont->GetHnEff());
178     fOutListCont->Add(fEffCont->GetHnCont());
179   }
180
181   // ------------------------------------------------------------------
182   // -- Add histograms from DCA class
183   // ------------------------------------------------------------------
184   if (fModeDCACreation == 1)
185     fOutListDCA->Add(fDCA->GetHnDCA());
186
187   // ------------------------------------------------------------------
188   // -- Add Multiplicity Correlations
189   // ------------------------------------------------------------------
190   /*
191   fOutList->Add(new TH2F("fHMultCorrTPCref0", "Centrality vs Multiplicity (TPCref);Centrality;Multiplicity (TPCref)",          
192                          24, -0.5, 11.49, 10001, 0.5, 10000.49));
193   fOutList->Add(new TH2F("fHMultCorrTracks0", "Centrality vs Multiplicity (Tracks);Centrality;Multiplicity (Tracks)",          
194                          24, -0.5, 11.49, 10001, 0.5, 10000.49));
195   fOutList->Add(new TH2F("fHMultCorrESDTracks0", "Centrality vs Multiplicity (ESDTracks);Centrality;Multiplicity (ESDTracks)", 
196                          24, -0.5, 11.49, 10001, 0.5, 10000.49));
197   
198   fOutList->Add(new TH2F("fHMultCorrTPCref1", "Centrality vs Multiplicity (TPCref);Centrality;Multiplicity (TPCref)",          
199                          101, -0.5, 100.49, 10001, 0.5, 10000.49));
200   fOutList->Add(new TH2F("fHMultCorrTracks1", "Centrality vs Multiplicity (Tracks);Centrality;Multiplicity (Tracks)",          
201                          101, -0.5, 100.49, 10001, 0.5, 10000.49));
202   fOutList->Add(new TH2F("fHMultCorrESDTracks1", "Centrality vs Multiplicity (ESDTracks);Centrality;Multiplicity (ESDTracks)", 
203                          101, -0.5, 100.49, 10001, 0.5, 10000.49));
204   */
205   // ------------------------------------------------------------------
206   // -- Create THnSparseF - QA
207   // ------------------------------------------------------------------
208   /*
209   Double_t dCent[2] = {-0.5, 8.5};
210   Int_t iCent       = 9;
211   
212   Double_t dEta[2]  = {-0.9, 0.9}; // -> 37 bins
213   Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
214
215   Double_t dRap[2]  = {-0.5, 0.5}; 
216   Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
217
218   Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
219   Int_t iPhi        = 42;
220
221   Double_t dPt[2]   = {0.1, 3.0}; 
222   Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
223
224   Double_t dSign[2] = {-1.5, 1.5};
225   Int_t iSign       = 3;
226
227   //                      cent:isAccepted: pInner:     pt:     sign:tpcSignal:nSigmaTPC:nSigmaTOF:      eta:     phi:       y: dcar: dcaz: nSigmaCdd: nSigmaCzz  
228   Int_t    bin[15] = {   iCent,       2  ,    iPt,    iPt,    iSign,      500,     50  ,     50  ,     iEta,    iPhi,    iRap,  50 ,  50 ,       50 ,       50 };
229   Double_t min[15] = {dCent[0],      -0.5, dPt[0], dPt[0], dSign[0],       30,     -5.0,     -5.0,  dEta[0], dPhi[0], dRap[0], -10., -10.,      -10.,      -10.};
230   Double_t max[15] = {dCent[1],       1.5, dPt[1], dPt[1], dSign[1],      500,      5.0,      5.0,  dEta[1], dPhi[1], dRap[1],  10.,  10.,       10.,       10.};
231   
232   if (fUseQATHnSparse) {
233     fOutListQA->Add(new THnSparseF("fHnQA", "cent:isAccepted:pInner:pt:sign:tpcSignal:nSigmaTPC:nSigmaTOF:eta:phi:u:dcar:dcaz:nSigmaCdd:nSigmaCzz",
234                                    15, bin, min, max));
235     
236     fHnQA = static_cast<THnSparseF*>(fOutListQA->Last());
237     fHnQA->Sumw2();
238
239     fHnQA->GetAxis(0)->SetTitle("centrality");
240     fHnQA->GetAxis(1)->SetTitle("isAccepted");
241     fHnQA->GetAxis(2)->SetTitle("#it{p}_{Inner} (GeV/#it{c})");
242     fHnQA->GetAxis(3)->SetTitle("#it{p}_{T} (GeV/#it{c})");
243     fHnQA->GetAxis(4)->SetTitle("sign");
244     
245     fHnQA->GetAxis(5)->SetTitle("TPC signal");
246     fHnQA->GetAxis(6)->SetTitle("n #sigma TPC");
247     fHnQA->GetAxis(7)->SetTitle("n #sigma TOF");
248     
249     fHnQA->GetAxis(8)->SetTitle("#eta");
250     fHnQA->GetAxis(9)->SetTitle("#varphi");
251     fHnQA->GetAxis(10)->SetTitle("#it{y}");    
252
253     fHnQA->GetAxis(11)->SetTitle("DCAr");
254     fHnQA->GetAxis(12)->SetTitle("DCAz");
255
256     fHnQA->GetAxis(13)->SetTitle("n #sigma #sqrt(Cdd)/DCAr");
257     fHnQA->GetAxis(14)->SetTitle("n #sigma #sqrt(Czz)/DCAz");
258
259     fHelper->BinLogAxis(fHnQA, 2);
260     fHelper->BinLogAxis(fHnQA, 3);
261   }
262   */
263   // ------------------------------------------------------------------
264
265   TH1::AddDirectory(oldStatus);
266
267   return;
268 }
269
270 //________________________________________________________________________
271 void AliAnalysisTaskNetParticle::UserExec(Option_t *) {
272   // Called for each event
273
274   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
275   // -- Setup Event
276   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
277   if (SetupEvent() < 0) {
278     PostData(1,fOutList);
279     PostData(2,fOutListEff);
280     PostData(3,fOutListCont);
281     PostData(4,fOutListDCA);
282     PostData(5,fOutListQA);
283     return;
284   }
285
286   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
287   // -- Process Efficiency / Contamination Determination
288   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
289   if ((fIsMC||fIsAOD) && fModeEffCreation == 1)
290     fEffCont->Process();
291
292   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
293   // -- Process DCA Determination
294   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
295   if (fModeDCACreation == 1)
296     fDCA->Process();
297
298   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
299   // -- Process Distributions 
300   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
301   if (fModeDistCreation == 1)
302     fDist->Process();
303
304   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
305   // -- Post output data
306   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
307
308   PostData(1,fOutList);
309   PostData(2,fOutListEff);
310   PostData(3,fOutListCont);
311   PostData(4,fOutListDCA);
312   PostData(5,fOutListQA);
313
314   return;
315 }      
316
317 //________________________________________________________________________
318 void AliAnalysisTaskNetParticle::Terminate(Option_t *){
319   // Terminate
320 }
321
322 /*
323  * ---------------------------------------------------------------------------------
324  *                                 Public Methods
325  * ---------------------------------------------------------------------------------
326  */
327
328 //________________________________________________________________________
329 Int_t AliAnalysisTaskNetParticle::Initialize() {
330   // Initialize event
331
332   // ------------------------------------------------------------------
333   // -- ESD TrackCuts
334   // ------------------------------------------------------------------
335   TString sModeName("");
336   
337   // -- Create ESD track cuts
338   // --------------------------
339   fESDTrackCutsBase->SetMinNCrossedRowsTPC(70);                                             // TPC
340   fESDTrackCutsBase->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8);                    // TPC
341
342   fESDTrackCutsBase->SetMaxChi2PerClusterTPC(4);                                            // TPC
343   fESDTrackCutsBase->SetAcceptKinkDaughters(kFALSE);                                        // TPC
344   fESDTrackCutsBase->SetRequireTPCRefit(kTRUE);                                             // TPC
345
346   fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kOff); // ITS
347   fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSDD,AliESDtrackCuts::kOff); // ITS
348   fESDTrackCutsBase->SetClusterRequirementITS(AliESDtrackCuts::kSSD,AliESDtrackCuts::kOff); // ITS
349
350   fESDTrackCutsBase->SetDCAToVertex2D(kFALSE);                                              // VertexConstrained 
351   fESDTrackCutsBase->SetRequireSigmaToVertex(kFALSE);                                       // VertexConstrained 
352
353   fESDTrackCutsBase->SetEtaRange(-1.*fEtaMax, fEtaMax);                                     // Acceptance
354   fESDTrackCutsBase->SetPtRange(fPtRange[0],fPtRange[1]);                                   // Acceptance
355
356   // -- Mode : clean cuts -> small contamination
357   if (fESDTrackCutMode == 0) {
358     sModeName = "Clean";
359     fESDTrackCutsBase->SetRequireITSRefit(kTRUE);                                           // ITS
360     fESDTrackCutsBase->SetMaxChi2PerClusterITS(36);                                         // ITS
361   }  
362   // -- Mode : dirty cuts -> high efficiency
363   else if (fESDTrackCutMode == 1) {
364     sModeName = "Dirty";
365     fESDTrackCutsBase->SetRequireITSRefit(kFALSE);                                          // ITS
366   }
367   // -- Mode : Default
368   else {
369     sModeName = "Base";
370   }
371
372   fESDTrackCutsBase->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
373
374   // -- Create ESD BKG track cuts
375   // ------------------------------
376   fESDTrackCutsBkg = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
377   fESDTrackCutsBkg->SetName(Form("NetParticleCuts2010_%s_Bkg",sModeName.Data()));
378   fESDTrackCutsBkg->SetMaxDCAToVertexZ(10.);                                                // VertexConstrained 
379   
380   // -- Create ESD track cuts
381   // ------------------------------
382   fESDTrackCuts = static_cast<AliESDtrackCuts*>(fESDTrackCutsBase->Clone());
383   fESDTrackCuts->SetName(Form("NetParticleCuts2010_%s",sModeName.Data()));
384   fESDTrackCuts->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");                         // VertexConstrained  ->  7*(0.0026+0.0050/pt^1.01)
385   fESDTrackCuts->SetMaxChi2TPCConstrainedGlobal(36);                                        // VertexConstrained
386   fESDTrackCuts->SetMaxDCAToVertexZ(2);                                                     // VertexConstrained 
387
388   // -- Create ESD Eff track cuts
389   // ------------------------------
390   fESDTrackCutsEff = static_cast<AliESDtrackCuts*>(fESDTrackCuts->Clone());
391   fESDTrackCutsEff->SetName(Form("NetParticleCuts2010_%s_Eff",sModeName.Data()));
392   fESDTrackCutsEff->SetPtRange(fPtRangeEff[0],fPtRangeEff[1]);                              // Acceptance
393
394   // ------------------------------------------------------------------
395   // -- Initialize Helper
396   // ------------------------------------------------------------------
397   if (fHelper->Initialize(fIsMC))
398     return -1;
399
400   // ------------------------------------------------------------------
401   // -- Create / Initialize Efficiency/Contamination
402   // ------------------------------------------------------------------
403   if ((fIsMC||fIsAOD) && fModeEffCreation == 1) {
404     fEffCont = new AliAnalysisNetParticleEffCont;
405     fEffCont->Initialize(fESDTrackCutsEff, fHelper,fAODtrackCutBit);
406   }
407
408   // ------------------------------------------------------------------
409   // -- Create / Initialize DCA Determination
410   // ------------------------------------------------------------------
411   if (fModeDCACreation == 1) {
412     fDCA = new AliAnalysisNetParticleDCA;
413     fDCA->Initialize(fESDTrackCuts, fESDTrackCutsBkg, fHelper);
414   }
415
416   // ------------------------------------------------------------------
417   // -- Create / Initialize DCA Determination
418   // ------------------------------------------------------------------
419   if (fModeDistCreation == 1) {
420     fDist = new AliAnalysisNetParticleDistribution;
421     fDist->Initialize(fHelper, fESDTrackCuts, fIsMC, fPtRange, fEtaMax,fAODtrackCutBit);
422   }
423
424   // ------------------------------------------------------------------
425   // -- Reset Event
426   // ------------------------------------------------------------------
427   ResetEvent();
428
429   return 0;
430 }
431
432 /*
433  * ---------------------------------------------------------------------------------
434  *                            Setup/Reset Methods - private
435  * ---------------------------------------------------------------------------------
436  */
437
438 //________________________________________________________________________
439 Int_t AliAnalysisTaskNetParticle::SetupEvent() {
440   // Setup Reading of event
441   // > return 0 for success / accepted event
442   // > return -1 for failed setup
443   // > return -2 for rejected event
444
445   ResetEvent();
446
447   // -- ESD Event
448   // ------------------------------------------------------------------
449   if(!fIsAOD){
450     if (SetupESDEvent() < 0) {
451       AliError("Setup ESD Event failed");
452       return -1;
453     }
454   }
455
456   // -- AOD Event
457   // ------------------------------------------------------------------
458   else {
459     if (SetupAODEvent() < 0) {
460       AliError("Setup AOD Event failed");
461       return -1;
462     }
463   }
464   
465   // -- Setup MC Event
466   // ------------------------------------------------------------------
467   if (fIsMC && SetupMCEvent() < 0) {
468     AliError("Setup MC Event failed");
469     return -1;
470   }
471
472   // -- Setup Event for Helper / EffCont  / DCA / Dist classes
473   // ------------------------------------------------------------------
474   fHelper->SetupEvent(fESDHandler, fAODHandler, fMCEvent);
475
476   if (fIsMC && fModeEffCreation)
477     fEffCont->SetupEvent(fESDHandler, fMCEvent);
478
479   if (fIsAOD && fModeEffCreation)
480     fEffCont->SetupEvent(fAODHandler);
481
482   if (fModeDCACreation == 1)
483     fDCA->SetupEvent(fESDHandler, fMCEvent);
484
485   if (fModeDistCreation == 1)
486     fDist->SetupEvent(fESDHandler, fAODHandler, fMCEvent);
487
488   // -- Evaluate Event cuts
489   // ------------------------------------------------------------------
490   return fHelper->IsEventRejected() ? -2 : 0;
491 }
492
493 //________________________________________________________________________
494 Int_t AliAnalysisTaskNetParticle::SetupESDEvent() {
495   // -- Setup ESD Event
496   // > return 0 for success 
497   // > return -1 for failed setup
498
499   fESDHandler= dynamic_cast<AliESDInputHandler*> 
500     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
501   if (!fESDHandler) {
502     AliError("Could not get ESD input handler");
503     return -1;
504   } 
505
506   fESD = fESDHandler->GetEvent();
507   if (!fESD) {
508     AliError("Could not get ESD event");
509     return -1;
510   }
511
512   // -- Check PID response
513   // ------------------------------------------------------------------
514   if (!fESDHandler->GetPIDResponse()) {
515     AliError("Could not get PID response");
516     return -1;
517   } 
518
519   // -- Check Vertex
520   // ------------------------------------------------------------------
521   if (!fESD->GetPrimaryVertexTracks()) {
522     AliError("Could not get vertex from tracks");
523     return -1;
524   }
525
526   // -- Check Centrality
527   // ------------------------------------------------------------------
528   if (!fESD->GetCentrality()) {
529     AliError("Could not get centrality");
530     return -1;
531   }
532
533   return 0;
534 }
535
536 //________________________________________________________________________
537 Int_t AliAnalysisTaskNetParticle::SetupAODEvent() {
538   // -- Setup AOD Event
539   // > return 0 for success 
540   // > return -1 for failed setup
541
542   fAODHandler= dynamic_cast<AliAODInputHandler*> 
543     (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
544   if (!fAODHandler) {
545     AliError("Could not get AOD input handler");
546     return -1;
547   } 
548
549   fAOD = fAODHandler->GetEvent();
550   if (!fAOD) {
551     AliError("Could not get AOD event");
552     return -1;
553   }
554
555   // -- Check PID response
556   // ------------------------------------------------------------------
557   if (!fAODHandler->GetPIDResponse()) {
558     AliError("Could not get PID response");
559     return -1;
560   } 
561
562   // -- Check Vertex
563   // ------------------------------------------------------------------
564   if (!fAOD->GetPrimaryVertex()) {
565     AliError("Could not get primary vertex");
566     return -1;
567   }
568
569   // -- Check Centrality
570   // ------------------------------------------------------------------
571   if (!fAOD->GetHeader()->GetCentralityP()) {
572     AliError("Could not get centrality");
573     return -1;
574   }
575
576   return 0;
577 }
578
579 //________________________________________________________________________
580 Int_t AliAnalysisTaskNetParticle::SetupMCEvent() {
581   // -- Setup MC Event
582   // > return 0 for success 
583   // > return -1 for failed setup
584
585   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> 
586     (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
587   
588   fMCEvent = mcH->MCEvent();
589   if (!fMCEvent) {
590     AliError("MC event not available");
591     return -1;
592   }
593
594   // -- Get MC header
595   // ------------------------------------------------------------------
596   AliHeader* header = fMCEvent->Header();
597   if (!header) {
598     AliError("MC header not available");
599     return -1;
600   }
601
602   // -- Check Stack
603   // ------------------------------------------------------------------
604   fMCStack = fMCEvent->Stack(); 
605   if (!fMCStack) {
606     AliError("MC stack not available");
607     return -1;
608   }
609     
610   // -- Check GenHeader
611   // ------------------------------------------------------------------
612   if (!header->GenEventHeader()) {
613     AliError("Could not retrieve genHeader from header");
614     return -1;
615   }
616
617   // -- Check primary vertex
618   // ------------------------------------------------------------------
619   if (!fMCEvent->GetPrimaryVertex()){
620     AliError("Could not get MC vertex");
621     return -1;
622   }
623
624   return 0;
625 }
626
627 //________________________________________________________________________
628 void AliAnalysisTaskNetParticle::ResetEvent() {
629   // -- Reset event
630   
631   // -- Reset ESD Event
632   fESD       = NULL;
633
634   // -- Reset ESD Event
635   fAOD       = NULL;
636
637   // -- Reset MC Event
638   if (fIsMC)
639     fMCEvent = NULL;
640
641   // -- Reset Dist Creation 
642   if (fModeDistCreation == 1)
643     fDist->ResetEvent();
644
645   return;
646 }
647
648 /*
649  * ---------------------------------------------------------------------------------
650  *                           Process Methods - private
651  * ---------------------------------------------------------------------------------
652  */
653
654 #if 0
655
656 //________________________________________________________________________
657 void AliAnalysisTaskNetParticle::ProcessESDTrackLoop() {
658   // -- Process ESD tracks and fill histograms
659   
660   Int_t nTracks = 0;
661
662   for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
663     AliESDtrack *track = fESD->GetTrack(idxTrack); 
664
665     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
666     // -- Check track cuts
667     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
668     
669     // -- Check if charged track is accepted for basic parameters
670     if (!fHelper->IsTrackAcceptedBasicCharged(track))
671       continue;
672     
673     // -- Check if accepted - BKG
674     if (!fESDTrackCutsBkg->AcceptTrack(track)) 
675       continue;
676
677     // -- Check if accepted in rapidity window
678     Double_t yP;
679     if (!fHelper->IsTrackAcceptedRapidity(track, yP))
680       continue;
681
682     // -- Check if accepted bt PID from TPC or TPC+TOF
683     Double_t pid[2];
684     Bool_t isAcceptedPID = fHelper->IsTrackAcceptedPID(track, pid);
685
686     // -- Check if accepted with thighter DCA cuts
687     Bool_t isAcceptedDCA = fHelper->IsTrackAcceptedDCA(track);
688
689     // -- Check track cuts
690     Bool_t isAcceptedVertex = fESDTrackCuts->AcceptTrack(track);
691
692     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
693     // -- Fill tracks in QA THnSparseF
694     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
695
696     if (fUseQATHnSparse && fHnQA->GetEntries() < 5e5) {
697
698       // -- Get dca r/z
699       Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
700       track->GetImpactParameters(dca,cov);
701       
702       Float_t dcaRoverCdd = ( TMath::Sqrt(cov[0]) != 0. )  ? dca[0]/TMath::Sqrt(cov[0]) : -9.99;
703       Float_t dcaZoverCzz = ( TMath::Sqrt(cov[2]) != 0. )  ? dca[1]/TMath::Sqrt(cov[2]) : -9.99;
704       
705       Double_t aQA[15] = {
706         Double_t(fHelper->GetCentralityBin()),                     //  0 centrality 
707         Double_t(isAcceptedVertex || isAcceptedDCA),  //  1 isAccepted -> Vertex || isAcceptedDCA
708         track->GetInnerParam()->GetP(),               //  2 p at InnerParam
709         track->Pt(),                                  //  3 pt
710         track->GetSign(),                             //  4 sign
711         track->GetTPCsignal(),                        //  5 TPC dE/dx
712         pid[0],                                       //  6 n Sigma TPC
713         pid[1],                                       //  7 n Sigma TOF
714         track->Eta(),                                 //  8 eta
715         track->Phi(),                                 //  9 phi
716         yP,                                           // 10 rapidity   
717         dca[0],                                       // 11 dca r
718         dca[1],                                       // 12 dca z
719         dcaRoverCdd,                                  // 13 sqrt(cov[dd])
720         dcaZoverCzz,                                  // 14 sqrt(cov[zz])
721       };
722     
723       fHnQA->Fill(aQA);
724     }
725
726     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
727     // -- Apply track cuts
728     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
729
730     if (!isAcceptedVertex || !isAcceptedDCA)
731       continue;
732
733     ++nTracks;
734
735     if (!isAcceptedPID)
736       continue;
737     
738   } // for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
739
740   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
741   // -- Fill Particle Fluctuation Histograms
742   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
743
744   /*
745   Int_t tpcRef = fESDTrackCuts->GetReferenceMultiplicity(fESD,kTRUE);
746
747   (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTPCref0")))->Fill(fHelper->GetCentralityBin(),  tpcRef);
748   (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTracks0")))->Fill(fHelper->GetCentralityBin(),  nTracks);
749
750   Float_t centPercentile = fHelper->GetCentralityPercentile();
751   (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTPCref1")))->Fill(centPercentile,  tpcRef);
752   (static_cast<TH2F*>(fOutList->FindObject("fHMultCorrTracks1")))->Fill(centPercentile,  nTracks);
753 */
754   return;
755 }
756
757 #endif