Fix in the last caall to CleanOwnPrimaryVertex
[u/mrichter/AliRoot.git] / ANALYSIS / AliEPSelectionTask.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //*****************************************************
17 //   Class AliEPSelectionTask
18 //   Class to determine event plane            
19 //   author: Alberica Toia, Johanna Gramling
20 //*****************************************************
21
22 #include "AliEPSelectionTask.h"
23
24 #include <TTree.h>
25 #include <TList.h>
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TProfile.h>
29 #include <TFile.h>
30 #include <TObjString.h>
31 #include <TString.h>
32 #include <TCanvas.h>
33 #include <TROOT.h>
34 #include <TDirectory.h>
35 #include <TSystem.h>
36 #include <iostream>
37 #include <TRandom2.h>
38 #include <TArrayF.h>
39
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.h"
42 #include "AliESD.h"
43 #include "AliESDEvent.h"
44 #include "AliMCEvent.h"
45 #include "AliESDtrack.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDHeader.h"
48 #include "AliESDInputHandler.h"
49 #include "AliAODHandler.h"
50 #include "AliAODEvent.h"
51 #include "AliHeader.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliPhysicsSelectionTask.h"
55 #include "AliPhysicsSelection.h"
56 #include "AliBackgroundSelection.h"
57 #include "AliESDUtils.h"
58 #include "AliOADBContainer.h"
59
60 #include "AliEventplane.h"
61
62 ClassImp(AliEPSelectionTask)
63
64 //________________________________________________________________________
65 AliEPSelectionTask::AliEPSelectionTask():
66 AliAnalysisTaskSE(),
67   fAnalysisInput("ESD"),
68   fTrackType("TPC"),
69   fUseMCRP(kFALSE),
70   fUsePhiWeight(kFALSE),
71   fUsePtWeight(kFALSE),
72   fSaveTrackContribution(kFALSE),
73   fUserphidist(kFALSE),
74   fUsercuts(kFALSE),
75   fRunNumber(-15),
76   fESDtrackCuts(0),
77   fEPContainer(0),
78   fPhiDist(0),
79   fQVector(0),
80   fQContributionX(0),
81   fQContributionY(0),
82   fEventplaneQ(0),
83   fQsub1(0),
84   fQsub2(0),
85   fQsubRes(0),  
86   fOutputList(0),
87   fHOutEventplaneQ(0),
88   fHOutPhi(0),
89   fHOutPhiCorr(0),
90   fHOutsub1sub2(0),
91   fHOutNTEPRes(0),
92   fHOutPTPsi(0),
93   fHOutDiff(0),
94   fHOutleadPTPsi(0)
95 {   
96   // Default constructor
97   AliInfo("Event Plane Selection enabled.");
98 }   
99
100 //________________________________________________________________________
101 AliEPSelectionTask::AliEPSelectionTask(const char *name):
102   AliAnalysisTaskSE(name),
103   fAnalysisInput("ESD"),
104   fTrackType("TPC"),
105   fUseMCRP(kFALSE),
106   fUsePhiWeight(kFALSE),
107   fUsePtWeight(kFALSE),  
108   fSaveTrackContribution(kFALSE),
109   fUserphidist(kFALSE),
110   fUsercuts(kFALSE),
111   fRunNumber(-15),
112   fESDtrackCuts(0),
113   fEPContainer(0),
114   fPhiDist(0),
115   fQVector(0),
116   fQContributionX(0),
117   fQContributionY(0),
118   fEventplaneQ(0),
119   fQsub1(0),
120   fQsub2(0),
121   fQsubRes(0),
122   fOutputList(0),
123   fHOutEventplaneQ(0),
124   fHOutPhi(0),
125   fHOutPhiCorr(0),
126   fHOutsub1sub2(0),
127   fHOutNTEPRes(0),
128   fHOutPTPsi(0),
129   fHOutDiff(0),
130   fHOutleadPTPsi(0)
131 {
132   // Default constructor
133   AliInfo("Event Plane Selection enabled.");
134   DefineOutput(1, TList::Class());
135 }
136  
137 //________________________________________________________________________
138 AliEPSelectionTask::~AliEPSelectionTask()
139 {
140   // Destructor  
141   if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
142       delete fOutputList;
143       fOutputList = 0;
144   }
145   if (fESDtrackCuts){
146       delete fESDtrackCuts;
147       fESDtrackCuts = 0;
148   }
149   if (fUserphidist) {
150     if (fPhiDist) {
151       delete fPhiDist;
152       fPhiDist = 0;
153     }
154   }
155   if (fEPContainer){
156       delete fEPContainer;
157       fEPContainer = 0;
158   }
159 }  
160
161 //________________________________________________________________________
162 void AliEPSelectionTask::UserCreateOutputObjects()
163 {  
164   // Create the output containers
165   if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
166   AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
167
168   fOutputList = new TList();
169   fOutputList->SetOwner();
170   fHOutEventplaneQ = new TH1F("fHOutEventplaneQ","fHOutEventplaneQ; Event Plane Q",100,0,TMath::Pi());
171   fHOutPhi = new TH1F("fHOutPhi","fHOutPhi; Phi Distribution",100,0,TMath::TwoPi());
172   fHOutPhiCorr = new TH1F("fHOutPhiCorr","fHOutPhiCorr; Corrected Phi Distribution",100,0,TMath::TwoPi());
173   fHOutsub1sub2 = new TH2F("fHOutsub1sub2","fHOutsub1sub2; EP1; EP2",100,0,TMath::Pi(),100,0,TMath::Pi());
174   fHOutNTEPRes = new TH2F("fHOutNTEPRes","fHOutNTEPRes; Number of Tracks; Event Plane Resolution",100,0,5000,100,-TMath::Pi(),TMath::Pi());
175   fHOutPTPsi = new TH2F("fHOutPTPsi","fHOutPTPsi; PT; Phi-EP",100,0,20,100,0,TMath::Pi());
176   fHOutDiff = new TH2F("fHOutDiff","fHOutDiff; EP; MCEP",100,0,TMath::Pi(),100,0,TMath::Pi());
177   fHOutleadPTPsi = new TH2F("fHOutleadPTPsi","fHOutleadPTPsi; leadPT; EP",100,0,TMath::Pi(),100,0,TMath::Pi());
178
179   fOutputList->Add(fHOutEventplaneQ);
180   fOutputList->Add(fHOutPhi);
181   fOutputList->Add(fHOutPhiCorr);
182   fOutputList->Add(fHOutsub1sub2);
183   fOutputList->Add(fHOutNTEPRes);
184   fOutputList->Add(fHOutPTPsi);
185   fOutputList->Add(fHOutleadPTPsi);
186   fOutputList->Add(fHOutDiff);
187   
188   PostData(1, fOutputList); 
189   
190     if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
191
192     TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
193
194     TFile foadb(oadbfilename); 
195     if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
196
197     AliInfo("Using Standard OADB");
198     fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");    
199     if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
200     foadb.Close();
201     }
202 }
203
204 //________________________________________________________________________
205 void AliEPSelectionTask::UserExec(Option_t */*option*/)
206
207   // Execute analysis for current event:
208   if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
209   
210 //   fRunNumber = -15;
211  
212   AliEventplane* esdEP = 0;
213   TVector2 qq1;
214   TVector2 qq2;
215   Double_t fRP = 0.; // the monte carlo reaction plane angle
216     
217   if (fAnalysisInput.CompareTo("ESD")==0){
218
219     AliVEvent* event = InputEvent();
220     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
221     if (esd){    
222       if (!(fRunNumber == esd->GetRunNumber())) {
223           fRunNumber = esd->GetRunNumber();
224             SetPhiDist();      
225       }
226       
227       
228       if (fUseMCRP) {
229           AliMCEvent* mcEvent  = MCEvent();      
230           if (mcEvent && mcEvent->GenEventHeader()) {
231               AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
232               if (headerH) fRP = headerH->ReactionPlaneAngle();
233           }
234       }
235       
236       esdEP = esd->GetEventplane();
237       if (fSaveTrackContribution) {
238         esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
239         esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
240       }
241       
242       TObjArray* tracklist = new TObjArray;
243       if (fTrackType.CompareTo("GLOBAL")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kFALSE);
244       if (fTrackType.CompareTo("TPC")==0) tracklist = fESDtrackCuts->GetAcceptedTracks(esd,kTRUE);
245       const int nt = tracklist->GetEntries();
246       
247       if (nt>4){
248         fQVector = new TVector2(GetQ(esdEP,tracklist));
249         fEventplaneQ = fQVector->Phi()/2; 
250         GetQsub(qq1, qq2, tracklist);
251         fQsub1 = new TVector2(qq1);
252         fQsub2 = new TVector2(qq2);
253         fQsubRes = (fQsub1->Phi()/2 - fQsub2->Phi()/2);
254         
255         esdEP->SetQVector(fQVector);
256         esdEP->SetEventplaneQ(fEventplaneQ);
257         esdEP->SetQsub(fQsub1,fQsub2);
258         esdEP->SetQsubRes(fQsubRes);
259         
260         fHOutEventplaneQ->Fill(fEventplaneQ);
261         fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
262         fHOutNTEPRes->Fill(nt,fQsubRes);
263
264         if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
265         
266         for (int iter = 0; iter<nt;iter++){
267           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
268           if (track) {
269             float delta = track->Phi()-fEventplaneQ;
270             while (delta < 0) delta += TMath::Pi();
271             while (delta > TMath::Pi()) delta -= TMath::Pi();
272             fHOutPTPsi->Fill(track->Pt(),delta);
273             fHOutPhi->Fill(track->Phi());
274             fHOutPhiCorr->Fill(track->Phi(),GetPhiWeight(track));
275           }
276         }
277         
278         AliESDtrack* trmax = esd->GetTrack(0);
279         for (int iter = 1; iter<nt;iter++){
280           AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
281           if (track && (track->Pt() > trmax->Pt())) trmax = track;
282         }
283         fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);      
284       }
285       delete tracklist;
286       tracklist = 0;
287     }
288   }
289   
290   else if (fAnalysisInput.CompareTo("AOD")==0){
291     //AliAODEvent *aod =  dynamic_cast<AliAODEvent*> (InputEvent());
292     // to be implemented
293     printf("  AOD analysis not yet implemented!!!\n\n");
294     return;
295   }  
296   else {
297     printf(" Analysis Input not known!\n\n ");
298     return;
299   }
300   PostData(1, fOutputList); 
301 }
302
303 //________________________________________________________________________
304 void AliEPSelectionTask::Terminate(Option_t */*option*/)
305 {
306   // Terminate analysis
307 }
308
309 //__________________________________________________________________________
310 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
311 {
312 // Get the Q vector
313   TVector2 mQ;
314   float mQx=0, mQy=0;
315   AliESDtrack* track;
316   Double_t weight;
317   
318   int nt = tracklist->GetEntries();
319
320   for (int i=0; i<nt; i++){
321     weight = 1;
322     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
323     if (track) {
324       weight = GetWeight(track);
325       if (fSaveTrackContribution){
326         EP->GetQContributionXArray()->AddAt(weight*cos(2*track->Phi()),track->GetID());
327         EP->GetQContributionYArray()->AddAt(weight*sin(2*track->Phi()),track->GetID());
328       }
329       mQx += (weight*cos(2*track->Phi()));
330       mQy += (weight*sin(2*track->Phi()));
331     }
332   }
333   mQ.Set(mQx,mQy);
334   return mQ;
335 }
336   
337   //________________________________________________________________________
338 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
339 {
340 // Get Qsub
341   TVector2 mQ[2];
342   float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
343   Double_t weight;
344
345   AliESDtrack* track;
346   TRandom2 rn = 0;
347   
348   int nt = tracklist->GetEntries();
349   int trackcounter1=0, trackcounter2=0;
350   
351   for (Int_t i = 0; i < nt; i++) {
352     weight = 1;
353     track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
354     if (!track) continue;
355     weight = GetWeight(track);
356     
357     // This loop splits the track set into 2 random subsets
358     if( trackcounter1 < int(nt/2.) && trackcounter2 < int(nt/2.)){
359       float random = rn.Rndm();
360       if(random < .5){
361         mQx1 += (weight*cos(2*track->Phi()));
362         mQy1 += (weight*sin(2*track->Phi()));
363         trackcounter1++;
364       }
365       else {
366         mQx2 += (weight*cos(2*track->Phi()));
367         mQy2 += (weight*sin(2*track->Phi()));
368         trackcounter2++;
369       }
370     }
371     else if( trackcounter1 >= int(nt/2.)){
372       mQx2 += (weight*cos(2*track->Phi()));
373       mQy2 += (weight*sin(2*track->Phi()));
374       trackcounter2++;
375     }
376     else {
377       mQx1 += (weight*cos(2*track->Phi()));
378       mQy1 += (weight*sin(2*track->Phi()));
379       trackcounter1++;
380     }
381   }
382   mQ[0].Set(mQx1,mQy1);
383   mQ[1].Set(mQx2,mQy2);
384   Q1 = mQ[0];
385   Q2 = mQ[1];
386 }
387
388 //________________________________________________________________________
389 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
390   
391   if(fESDtrackCuts){ 
392     delete fESDtrackCuts;
393     fESDtrackCuts = 0;
394   }
395     
396   fUsercuts = kTRUE;
397   fESDtrackCuts = trackcuts;
398 }
399
400 //_____________________________________________________________________________
401 void AliEPSelectionTask::SetTrackType(TString tracktype){
402 // Set the track type
403   fTrackType = tracktype;
404   if (!fUsercuts) {
405   if (fTrackType.CompareTo("GLOBAL")==0) fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
406   if (fTrackType.CompareTo("TPC")==0)    fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
407   fESDtrackCuts->SetPtRange(0.15,20);
408   fESDtrackCuts->SetEtaRange(-0.8,0.8);
409   }
410 }
411
412 //________________________________________________________________________
413 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
414 {
415 // Get weight for track
416   Double_t ptweight=1;
417
418   if (fUsePtWeight) {      
419     if (track->Pt()<2) ptweight=track->Pt();
420     else ptweight=2;
421   }
422   return ptweight*GetPhiWeight(track);
423 }
424
425 //________________________________________________________________________
426 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
427 {
428 // Get phi weight for track
429   Double_t phiweight=1;
430   
431   if (fUsePhiWeight && fPhiDist && track) {
432     Double_t nParticles = fPhiDist->Integral();
433     Double_t nPhibins = fPhiDist->GetNbinsX();
434   
435     Double_t phi = track->Phi();
436     
437     while (phi<0) phi += TMath::TwoPi();
438     while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
439       
440     Double_t phiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
441     
442     if (phiDistValue > 0) phiweight = nParticles/nPhibins/phiDistValue;
443   }
444   return phiweight;
445 }
446
447 //__________________________________________________________________________
448 void AliEPSelectionTask::SetPhiDist() 
449 {
450 // Set the phi distribution
451   if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
452
453     fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
454     if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
455
456   } 
457   else {
458     AliInfo("Using Custom Phi Distribution");
459   }
460     
461   Bool_t emptybins;
462
463   int iter = 0;  
464   while (iter<3){
465       emptybins = kFALSE;
466    
467       for (int i=1; i<fPhiDist->GetNbinsX(); i++){
468         if (!((fPhiDist->GetBinContent(i))>0)) {
469           emptybins = kTRUE;
470         }
471       }  
472       if (emptybins) {
473         cout << "empty bins - rebinning!" << endl;
474         fPhiDist->Rebin();
475         iter++;
476       }      
477       else iter = 3;
478   }
479   
480   if (emptybins) {
481     AliError("After Maximum of rebinning still empty Phi-bins!!!");
482   }
483 }
484
485 //__________________________________________________________________________
486 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
487 {
488     // Set a personal phi distribution
489   fUserphidist = kTRUE;
490   
491   TFile f(infilename);
492   TObject* list = f.Get(listname);
493   fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
494   if (!fPhiDist) AliFatal("Phi Distribution not found!!!");
495
496   f.Close();
497