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