1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //*****************************************************
17 // Class AliEPSelectionTask
18 // Class to determine event plane
19 // author: Alberica Toia, Johanna Gramling
20 //*****************************************************
22 #include "AliEPSelectionTask.h"
30 #include <TObjString.h>
34 #include <TDirectory.h>
40 #include "AliAnalysisManager.h"
41 #include "AliVEvent.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"
60 #include "AliEventplane.h"
62 ClassImp(AliEPSelectionTask)
64 //________________________________________________________________________
65 AliEPSelectionTask::AliEPSelectionTask():
67 fAnalysisInput("ESD"),
70 fUsePhiWeight(kFALSE),
72 fSaveTrackContribution(kFALSE),
96 // Default constructor
97 AliInfo("Event Plane Selection enabled.");
100 //________________________________________________________________________
101 AliEPSelectionTask::AliEPSelectionTask(const char *name):
102 AliAnalysisTaskSE(name),
103 fAnalysisInput("ESD"),
106 fUsePhiWeight(kFALSE),
107 fUsePtWeight(kFALSE),
108 fSaveTrackContribution(kFALSE),
109 fUserphidist(kFALSE),
132 // Default constructor
133 AliInfo("Event Plane Selection enabled.");
134 DefineOutput(1, TList::Class());
137 //________________________________________________________________________
138 AliEPSelectionTask::~AliEPSelectionTask()
141 if (fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()){
146 delete fESDtrackCuts;
161 //________________________________________________________________________
162 void AliEPSelectionTask::UserCreateOutputObjects()
164 // Create the output containers
165 if (fDebug>1) printf("AliEPSelectionTask::UserCreateOutputObjects() \n");
166 AliLog::SetClassDebugLevel("AliEPSelectionTask", AliLog::kInfo);
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());
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);
188 PostData(1, fOutputList);
190 if(!fUserphidist) { // if it's already set and custom class is required, we use the one provided by the user
192 TString oadbfilename = (Form("%s/COMMON/EVENTPLANE/data/epphidist.root", AliAnalysisManager::GetOADBPath()));
194 TFile foadb(oadbfilename);
195 if(!foadb.IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
197 AliInfo("Using Standard OADB");
198 fEPContainer = (AliOADBContainer*) foadb.Get("epphidist");
199 if (!fEPContainer) AliFatal("Cannot fetch OADB container for EP selection");
204 //________________________________________________________________________
205 void AliEPSelectionTask::UserExec(Option_t */*option*/)
207 // Execute analysis for current event:
208 if (fDebug>1) printf(" **** AliEPSelectionTask::UserExec() \n");
212 AliEventplane* esdEP = 0;
215 Double_t fRP = 0.; // the monte carlo reaction plane angle
217 if (fAnalysisInput.CompareTo("ESD")==0){
219 AliVEvent* event = InputEvent();
220 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
222 if (!(fRunNumber == esd->GetRunNumber())) {
223 fRunNumber = esd->GetRunNumber();
229 AliMCEvent* mcEvent = MCEvent();
230 if (mcEvent && mcEvent->GenEventHeader()) {
231 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
232 if (headerH) fRP = headerH->ReactionPlaneAngle();
236 esdEP = esd->GetEventplane();
237 if (fSaveTrackContribution) {
238 esdEP->GetQContributionXArray()->Set(esd->GetNumberOfTracks());
239 esdEP->GetQContributionYArray()->Set(esd->GetNumberOfTracks());
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();
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);
255 esdEP->SetQVector(fQVector);
256 esdEP->SetEventplaneQ(fEventplaneQ);
257 esdEP->SetQsub(fQsub1,fQsub2);
258 esdEP->SetQsubRes(fQsubRes);
260 fHOutEventplaneQ->Fill(fEventplaneQ);
261 fHOutsub1sub2->Fill(fQsub1->Phi()/2,fQsub2->Phi()/2);
262 fHOutNTEPRes->Fill(nt,fQsubRes);
264 if (fUseMCRP) fHOutDiff->Fill(fEventplaneQ, fRP);
266 for (int iter = 0; iter<nt;iter++){
267 AliESDtrack* track = dynamic_cast<AliESDtrack*> (tracklist->At(iter));
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));
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;
283 fHOutleadPTPsi->Fill(trmax->Phi(),fEventplaneQ);
290 else if (fAnalysisInput.CompareTo("AOD")==0){
291 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
293 printf(" AOD analysis not yet implemented!!!\n\n");
297 printf(" Analysis Input not known!\n\n ");
300 PostData(1, fOutputList);
303 //________________________________________________________________________
304 void AliEPSelectionTask::Terminate(Option_t */*option*/)
306 // Terminate analysis
309 //__________________________________________________________________________
310 TVector2 AliEPSelectionTask::GetQ(AliEventplane* EP, TObjArray* tracklist)
318 int nt = tracklist->GetEntries();
320 for (int i=0; i<nt; i++){
322 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
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());
329 mQx += (weight*cos(2*track->Phi()));
330 mQy += (weight*sin(2*track->Phi()));
337 //________________________________________________________________________
338 void AliEPSelectionTask::GetQsub(TVector2 &Q1, TVector2 &Q2, TObjArray* tracklist)
342 float mQx1=0, mQy1=0, mQx2=0, mQy2=0;
348 int nt = tracklist->GetEntries();
349 int trackcounter1=0, trackcounter2=0;
351 for (Int_t i = 0; i < nt; i++) {
353 track = dynamic_cast<AliESDtrack*> (tracklist->At(i));
354 if (!track) continue;
355 weight = GetWeight(track);
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();
361 mQx1 += (weight*cos(2*track->Phi()));
362 mQy1 += (weight*sin(2*track->Phi()));
366 mQx2 += (weight*cos(2*track->Phi()));
367 mQy2 += (weight*sin(2*track->Phi()));
371 else if( trackcounter1 >= int(nt/2.)){
372 mQx2 += (weight*cos(2*track->Phi()));
373 mQy2 += (weight*sin(2*track->Phi()));
377 mQx1 += (weight*cos(2*track->Phi()));
378 mQy1 += (weight*sin(2*track->Phi()));
382 mQ[0].Set(mQx1,mQy1);
383 mQ[1].Set(mQx2,mQy2);
388 //________________________________________________________________________
389 void AliEPSelectionTask::SetPersonalESDtrackCuts(AliESDtrackCuts* trackcuts){
392 delete fESDtrackCuts;
397 fESDtrackCuts = trackcuts;
400 //_____________________________________________________________________________
401 void AliEPSelectionTask::SetTrackType(TString tracktype){
402 // Set the track type
403 fTrackType = tracktype;
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);
412 //________________________________________________________________________
413 Double_t AliEPSelectionTask::GetWeight(AliESDtrack* track)
415 // Get weight for track
419 if (track->Pt()<2) ptweight=track->Pt();
422 return ptweight*GetPhiWeight(track);
425 //________________________________________________________________________
426 Double_t AliEPSelectionTask::GetPhiWeight(AliESDtrack* track)
428 // Get phi weight for track
429 Double_t phiweight=1;
431 if (fUsePhiWeight && fPhiDist && track) {
432 Double_t nParticles = fPhiDist->Integral();
433 Double_t nPhibins = fPhiDist->GetNbinsX();
435 Double_t phi = track->Phi();
437 while (phi<0) phi += TMath::TwoPi();
438 while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
440 Double_t phiDistValue = fPhiDist->GetBinContent(1+TMath::FloorNint((track->Phi())*nPhibins/TMath::TwoPi()));
442 if (phiDistValue > 0) phiweight = nParticles/nPhibins/phiDistValue;
447 //__________________________________________________________________________
448 void AliEPSelectionTask::SetPhiDist()
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
453 fPhiDist = (TH1F*) fEPContainer->GetObject(fRunNumber, "Default");
454 if (!fPhiDist) AliFatal(Form("Cannot find OADB phi distribution for run %d", fRunNumber));
458 AliInfo("Using Custom Phi Distribution");
467 for (int i=1; i<fPhiDist->GetNbinsX(); i++){
468 if (!((fPhiDist->GetBinContent(i))>0)) {
473 cout << "empty bins - rebinning!" << endl;
481 AliError("After Maximum of rebinning still empty Phi-bins!!!");
485 //__________________________________________________________________________
486 void AliEPSelectionTask::SetPersonalPhiDistribution(const char* infilename, char* listname)
488 // Set a personal phi distribution
489 fUserphidist = kTRUE;
492 TObject* list = f.Get(listname);
493 fPhiDist = (TH1F*)list->FindObject("fHOutPhi");
494 if (!fPhiDist) AliFatal("Phi Distribution not found!!!");