AliPhysics  b555aef (b555aef)
AliAnalysisTaskFlowEvent.cxx
Go to the documentation of this file.
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 
17 // AliAnalysisTaskFlowEvent:
18 //
19 // analysis task for filling the flow event
20 // from MCEvent, ESD, AOD ....
21 // and put it in an output stream so it can
22 // be used by the various flow analysis methods
23 // for cuts the correction framework is used
24 // which also outputs QA histograms to view
25 // the effects of the cuts
27 
28 #include "Riostream.h" //needed as include
29 #include "TChain.h"
30 #include "TTree.h"
31 #include "TFile.h" //needed as include
32 #include "TList.h"
33 #include "TF1.h"
34 #include "TH2F.h"
35 #include "TRandom3.h"
36 #include "TTimeStamp.h"
37 #include "AliAODHandler.h"
38 
39 // ALICE Analysis Framework
40 #include "AliAnalysisManager.h"
41 #include "AliAnalysisTaskSE.h"
42 #include "AliESDtrackCuts.h"
43 
44 // ESD interface
45 #include "AliESDEvent.h"
46 #include "AliESDInputHandler.h"
47 
48 // AOD interface
49 #include "AliAODEvent.h"
50 #include "AliAODInputHandler.h"
51 
52 // Monte Carlo Event
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 
56 // ALICE Correction Framework
57 #include "AliCFManager.h"
58 
59 // Interface to Event generators to get Reaction Plane Angle
60 #include "AliGenCocktailEventHeader.h"
61 #include "AliGenHijingEventHeader.h"
62 #include "AliGenGeVSimEventHeader.h"
63 #include "AliGenEposEventHeader.h"
64 
65 // Interface to Load short life particles
66 #include "TObjArray.h"
67 #include "AliFlowCandidateTrack.h"
68 
69 // Interface to make the Flow Event Simple used in the flow analysis methods
70 #include "AliFlowEvent.h"
71 #include "AliFlowTrackCuts.h"
72 #include "AliFlowEventCuts.h"
73 #include "AliFlowCommonConstants.h"
75 
76 #include "AliLog.h"
77 #include "AliAODHeader.h"
78 #include "AliAODVertex.h"
79 #include "AliHeader.h"
80 
81 using std::cout;
82 using std::endl;
84 
85 //________________________________________________________________________
88  // fOutputFile(NULL),
89  fAnalysisType("AUTOMATIC"),
90  fRPType("Global"),
91  fCFManager1(NULL),
92  fCFManager2(NULL),
93  fCutsEvent(NULL),
94  fCutsRP(NULL),
95  fCutsPOI(NULL),
96  fCutContainer(NULL),
97  fQAList(NULL),
98  fMinMult(0),
99  fMaxMult(10000000),
100  fMinA(-1.0),
101  fMaxA(-0.01),
102  fMinB(0.01),
103  fMaxB(1.0),
104  fQAon(kFALSE),
105  fLoadCandidates(kFALSE),
106  fPassMCeventToCutsObject(kTRUE),
107  fNbinsMult(10000),
108  fNbinsPt(100),
109  fNbinsPhi(100),
110  fNbinsEta(200),
111  fNbinsQ(500),
112  fNbinsMass(1),
113  fMultMin(0.),
114  fMultMax(10000.),
115  fPtMin(0.),
116  fPtMax(10.),
117  fPhiMin(0.),
118  fPhiMax(TMath::TwoPi()),
119  fEtaMin(-5.),
120  fEtaMax(5.),
121  fQMin(0.),
122  fQMax(3.),
123  fMassMin(-1.),
124  fMassMax(0.),
125  fHistWeightvsPhiMin(0.),
126  fHistWeightvsPhiMax(3.),
127  fExcludedEtaMin(0.),
128  fExcludedEtaMax(0.),
129  fExcludedPhiMin(0.),
130  fExcludedPhiMax(0.),
131  fAfterburnerOn(kFALSE),
132  fNonFlowNumberOfTrackClones(0),
133  fV1(0.),
134  fV2(0.),
135  fV3(0.),
136  fV4(0.),
137  fV5(0.),
138  fDifferentialV2(0),
139  fFlowEvent(NULL),
140  fShuffleTracks(kFALSE),
141  fMyTRandom3(NULL)
142 {
143  // Constructor
144  AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()");
145 }
146 
147 //________________________________________________________________________
148 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
149  AliAnalysisTaskSE(name),
150  // fOutputFile(NULL),
151  fAnalysisType("AUTOMATIC"),
152  fRPType(RPtype),
153  fCFManager1(NULL),
154  fCFManager2(NULL),
155  fCutsEvent(NULL),
156  fCutsRP(NULL),
157  fCutsPOI(NULL),
158  fCutContainer(new TList()),
159  fQAList(NULL),
160  fMinMult(0),
161  fMaxMult(10000000),
162  fMinA(-1.0),
163  fMaxA(-0.01),
164  fMinB(0.01),
165  fMaxB(1.0),
166  fQAon(on),
167  fLoadCandidates(bCandidates),
168  fPassMCeventToCutsObject(kTRUE),
169  fNbinsMult(10000),
170  fNbinsPt(100),
171  fNbinsPhi(100),
172  fNbinsEta(200),
173  fNbinsQ(500),
174  fNbinsMass(1),
175  fMultMin(0.),
176  fMultMax(10000.),
177  fPtMin(0.),
178  fPtMax(10.),
179  fPhiMin(0.),
180  fPhiMax(TMath::TwoPi()),
181  fEtaMin(-5.),
182  fEtaMax(5.),
183  fQMin(0.),
184  fQMax(3.),
185  fMassMin(-1.),
186  fMassMax(0.),
187  fHistWeightvsPhiMin(0.),
188  fHistWeightvsPhiMax(3.),
189  fExcludedEtaMin(0.),
190  fExcludedEtaMax(0.),
191  fExcludedPhiMin(0.),
192  fExcludedPhiMax(0.),
193  fAfterburnerOn(kFALSE),
194  fNonFlowNumberOfTrackClones(0),
195  fV1(0.),
196  fV2(0.),
197  fV3(0.),
198  fV4(0.),
199  fV5(0.),
200  fDifferentialV2(0),
201  fFlowEvent(NULL),
202  fShuffleTracks(kFALSE),
203  fMyTRandom3(NULL)
204 {
205  // Constructor
206  AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)");
207  fMyTRandom3 = new TRandom3(iseed);
208  gRandom->SetSeed(fMyTRandom3->Integer(65539));
209 
210  int availableINslot=1;
211  //FMD input slot
212  if (strcmp(RPtype,"FMD")==0) {
213  DefineInput(availableINslot++, TList::Class());
214  }
215  //Candidates input slot
216  if( fLoadCandidates )
217  DefineInput(availableINslot, TObjArray::Class());
218 
219  // Define output slots here
220  // Define here the flow event output
221  DefineOutput(1, AliFlowEventSimple::Class());
222  DefineOutput(2, TList::Class());
223 
224  // and for testing open an output file
225  // fOutputFile = new TFile("FlowEvents.root","RECREATE");
226 
227 }
228 
229 //________________________________________________________________________
231 {
232  //
233  // Destructor
234  //
235  delete fMyTRandom3;
236  delete fFlowEvent;
237  delete fCutsEvent;
238  delete fQAList;
239  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
240  // objects in the output list are deleted
241  // by the TSelector dtor (I hope)
242 
243 }
244 
245 //________________________________________________________________________
247 {
248  //at the beginning of each new run
249  if (fCutsRP) fCutsRP->SetRunsMuon(fInputHandler); // XZhang 20120604
250  if (fCutsPOI) fCutsPOI->SetRunsMuon(fInputHandler); // XZhang 20120604
251  AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
252  if (!fESD) return;
253 
254  Int_t run = fESD->GetRunNumber();
255  AliInfo(Form("Stariting run #%i",run));
256 }
257 
258 //________________________________________________________________________
260 {
261  // Called at every worker node to initialize
262  AliDebug(2,"AliAnalysisTaskFlowEvent::CreateOutputObjects()");
263 
264  if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
265  {
266  AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
267  exit(1);
268  }
269 
270  //set the common constants
273  cc->SetNbinsPt(fNbinsPt);
274  cc->SetNbinsPhi(fNbinsPhi);
275  cc->SetNbinsEta(fNbinsEta);
276  cc->SetNbinsQ(fNbinsQ);
278  cc->SetMultMin(fMultMin);
279  cc->SetMultMax(fMultMax);
280  cc->SetPtMin(fPtMin);
281  cc->SetPtMax(fPtMax);
282  cc->SetPhiMin(fPhiMin);
283  cc->SetPhiMax(fPhiMax);
284  cc->SetEtaMin(fEtaMin);
285  cc->SetEtaMax(fEtaMax);
286  cc->SetQMin(fQMin);
287  cc->SetQMax(fQMax);
288  cc->SetMassMin(fMassMin);
289  cc->SetMassMax(fMassMax);
292 
293  fFlowEvent = new AliFlowEvent(10000);
294 
295  if (fQAon)
296  {
297  fQAList=new TList();
298  fQAList->SetOwner();
299  fQAList->SetName(Form("%s QA",GetName()));
300  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
301  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
302  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
303  fQAList->Add(new TH1F("event plane angle","event plane angle;angle [rad];",100,0.,TMath::TwoPi())); //3
304  PostData(2,fQAList);
305  }
306 }
307 
308 //________________________________________________________________________
309 
311 {
312  // Main loop
313  // Called for each event
314  //delete fFlowEvent;
315  AliMCEvent* mcEvent = MCEvent(); // from TaskSE
316  AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
317  AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
318  AliMultiplicity* myTracklets = NULL;
319  AliESDPmdTrack* pmdtracks = NULL;//pmd
320 
321  int availableINslot=1;
322 
323  if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
324  {
325  AliError("cuts not set");
326  return;
327  }
328 
329  //DEFAULT - automatically takes care of everything
330  if (fAnalysisType == "AUTOMATIC")
331  {
332  //check event cuts
333  if (InputEvent() && !fCutsEvent->IsSelected(InputEvent(),MCEvent()))
334  return;
335 
336 
338  //first attach all possible information to the cuts
339  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
340  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
341  }
342  else{
343  AliMCEvent* McEventFake = NULL;
344  fCutsRP->SetEvent( InputEvent(), McEventFake ); //attach event
345  fCutsPOI->SetEvent( InputEvent(), McEventFake);
346  }
347 
348 
349 
350  //then make the event
352  //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
353 
354  // if (myESD)
356  fFlowEvent->SetCentrality(fCutsEvent->GetCentrality(InputEvent(),mcEvent));
357  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
358  }
359 
360  // Make the FlowEvent for MC input
361  else if (fAnalysisType == "MC")
362  {
363  // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
364  // This handler can return the current MC event
365  if (!(fCFManager1&&fCFManager2))
366  {
367  AliError("ERROR: No pointer to correction framework cuts! ");
368  return;
369  }
370  if (!mcEvent)
371  {
372  AliError("ERROR: Could not retrieve MC event");
373  return;
374  }
375 
376  fCFManager1->SetMCEventInfo(mcEvent);
377  fCFManager2->SetMCEventInfo(mcEvent);
378 
379  AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
380 
381  //check multiplicity
382  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
383  {
384  AliWarning("Event does not pass multiplicity cuts"); return;
385  }
386  //make event
388  }
389 
390  // Make the FlowEvent for ESD input
391  else if (fAnalysisType == "ESD")
392  {
393  if (!(fCFManager1&&fCFManager2))
394  {
395  AliError("ERROR: No pointer to correction framework cuts!");
396  return;
397  }
398  if (!myESD)
399  {
400  AliError("ERROR: ESD not available");
401  return;
402  }
403 
404  //check the offline trigger (check if the event has the correct trigger)
405  AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
406 
407  //check multiplicity
408  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
409  {
410  AliWarning("Event does not pass multiplicity cuts"); return;
411  }
412 
413  //make event
414  if (fRPType == "Global") {
416  }
417  else if (fRPType == "TPCOnly") {
418  fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
419  }
420  else if (fRPType == "TPCHybrid") {
421  fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
422  }
423  else if (fRPType == "Tracklet"){
424  fFlowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
425  }
426  else if (fRPType == "FMD"){
427  TList* dataFMD = dynamic_cast<TList*>(GetInputData(availableINslot++));
428  if(!dataFMD) {
429  cout<<" No dataFMD "<<endl;
430  return;
431  }
432  TH2F* histFMD = dynamic_cast<TH2F*>(dataFMD->FindObject("dNdetadphiHistogramTrVtx"));
433  if (!histFMD) {
434  cout<< "No histFMD"<<endl;
435  return;
436  }
437  fFlowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
438  }
439  else if (fRPType == "PMD"){
440  fFlowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
441  }
442  else return;
443 
444  // if monte carlo event get reaction plane from monte carlo (depends on generator)
445  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
446  //set reference multiplicity, TODO: maybe move it to the constructor?
447  fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
448  }
449 
450  // Make the FlowEvent for ESD input combined with MC info
451  else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
452  {
453  if (!(fCFManager1&&fCFManager2))
454  {
455  AliError("ERROR: No pointer to correction framework cuts! ");
456  return;
457  }
458  if (!myESD)
459  {
460  AliError("ERROR: ESD not available");
461  return;
462  }
463  AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
464 
465  if (!mcEvent)
466  {
467  AliError("ERROR: Could not retrieve MC event");
468  return;
469  }
470 
471  fCFManager1->SetMCEventInfo(mcEvent);
472  fCFManager2->SetMCEventInfo(mcEvent);
473 
474  //check multiplicity
475  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
476  {
477  AliWarning("Event does not pass multiplicity cuts"); return;
478  }
479 
480  //make event
481  if (fAnalysisType == "ESDMCkineESD")
482  {
484  }
485  else if (fAnalysisType == "ESDMCkineMC")
486  {
488  }
489  if (!fFlowEvent) return;
490  // if monte carlo event get reaction plane from monte carlo (depends on generator)
491  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
492  //set reference multiplicity, TODO: maybe move it to the constructor?
493  fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
494  }
495  // Make the FlowEventSimple for AOD input
496  else if (fAnalysisType == "AOD")
497  {
498  if (!myAOD)
499  {
500  AliError("ERROR: AOD not available");
501  return;
502  }
503  AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
504  fFlowEvent = new AliFlowEvent(myAOD);
505  }
506 
507  //inject candidates
508  if(fLoadCandidates) {
509  TObjArray* candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
510  //if(candidates->GetEntriesFast())
511  // printf("I received %d candidates\n",candidates->GetEntriesFast());
512  if (candidates)
513  {
514  for(int iCand=0; iCand!=candidates->GetEntriesFast(); ++iCand ) {
515  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(candidates->At(iCand));
516  if (!cand) continue;
517  //printf(" - Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
518  for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
519  //printf(" - Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
520  for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs ) {
521  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
522  if (!iRP) continue;
523  if( !iRP->InRPSelection() )
524  continue;
525  if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
526  //printf(" was in RP set");
527  //cand->SetDaughter( iDau, iRP );
528  //temporarily untagging all daugters
529  iRP->SetForRPSelection(kFALSE);
531  }
532  }
533  //printf("\n");
534  }
535  cand->SetForPOISelection(kTRUE);
536  fFlowEvent->InsertTrack( ((AliFlowTrack*) cand) );
537  }
538  }
539  }
540 
541  if (!fFlowEvent) return; //shuts up coverity
542 
543  //check final event cuts
544  Int_t mult = fFlowEvent->NumberOfTracks();
545  // AliInfo(Form("FlowEvent has %i tracks",mult));
546  if (mult<fMinMult || mult>fMaxMult)
547  {
548  AliWarning("FlowEvent cut on multiplicity"); return;
549  }
550 
551  //define dead zone
553 
554 
557  if (fAfterburnerOn)
558  {
559  //if reaction plane not set from elsewhere randomize it before adding flow
561  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
562 
563  if(fDifferentialV2)
565  else
566  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
567  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
568  }
570 
571  //tag subEvents
573 
574  //QA
575  if (fQAon)
576  {
577  TH1* h1 = static_cast<TH1*>(fQAList->FindObject("event plane angle"));
578  h1->Fill(fFlowEvent->GetMCReactionPlaneAngle());
579  }
580 
581  //do we want to serve shullfed tracks to everybody?
583 
584  // associate the mother particles to their daughters in the flow event (if any)
586 
587  //fListHistos->Print();
588  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
589  PostData(1,fFlowEvent);
590 }
591 
592 //________________________________________________________________________
594 {
595  // Called once at the end of the query -- do not call in case of CAF
596 }
597 
virtual void UserExec(Option_t *option)
void FindDaughters(Bool_t keepDaughtersInRPselection=kFALSE)
const Color_t cc[]
Definition: DrawKs.C:1
Float_t GetCentrality(AliVEvent *event, AliMCEvent *mcEvent)
Definition: External.C:236
Int_t GetNumberOfRPs() const
void AddFlow(Double_t v1, Double_t v2, Double_t v3, Double_t v4, Double_t v5)
Int_t GetReferenceMultiplicity(AliVEvent *event, AliMCEvent *mcEvent)
void SetMCReactionPlaneAngle(const AliMCEvent *mcEvent)
void SetHistWeightvsPhiMax(Double_t d)
void SetReferenceMultiplicity(Int_t m)
TRandom * gRandom
void SetRunsMuon(const AliVEventHandler *eventHandler)
Bool_t InRPSelection() const
AliFlowTrack * GetTrack(Int_t i)
TList * GetQA() const
void SetForRPSelection(Bool_t b=kTRUE)
Int_t GetNDaughters(void) const
void SetNumberOfRPs(Int_t nr)
virtual void Terminate(Option_t *)
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
void SetShuffleTracks(Bool_t b)
static AliFlowCommonConstants * GetMaster()
void SetCentrality(Double_t c)
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
Int_t GetIDDaughter(Int_t value) const
void SetEvent(AliVEvent *event, AliMCEvent *mcEvent=NULL)
Bool_t IsSetMCReactionPlaneAngle() const
TList * GetQA() const
void SetForPOISelection(Bool_t b=kTRUE)
void SetHistWeightvsPhiMin(Double_t d)
void AddV2(Double_t v2)
void TagSubeventsInEta(Double_t etaMinA, Double_t etaMaxA, Double_t etaMinB, Double_t etaMaxB)
Double_t GetMCReactionPlaneAngle() const
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
virtual Bool_t IsSelected(TObject *obj, TObject *objmc)
void InsertTrack(AliFlowTrack *)
Definition: External.C:196
Int_t NumberOfTracks() const