AliPhysics  ff1d528 (ff1d528)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  fNbinsMult(10000),
107  fNbinsPt(100),
108  fNbinsPhi(100),
109  fNbinsEta(200),
110  fNbinsQ(500),
111  fNbinsMass(1),
112  fMultMin(0.),
113  fMultMax(10000.),
114  fPtMin(0.),
115  fPtMax(10.),
116  fPhiMin(0.),
117  fPhiMax(TMath::TwoPi()),
118  fEtaMin(-5.),
119  fEtaMax(5.),
120  fQMin(0.),
121  fQMax(3.),
122  fMassMin(-1.),
123  fMassMax(0.),
124  fHistWeightvsPhiMin(0.),
125  fHistWeightvsPhiMax(3.),
126  fExcludedEtaMin(0.),
127  fExcludedEtaMax(0.),
128  fExcludedPhiMin(0.),
129  fExcludedPhiMax(0.),
130  fAfterburnerOn(kFALSE),
131  fNonFlowNumberOfTrackClones(0),
132  fV1(0.),
133  fV2(0.),
134  fV3(0.),
135  fV4(0.),
136  fV5(0.),
137  fDifferentialV2(0),
138  fFlowEvent(NULL),
139  fShuffleTracks(kFALSE),
140  fMyTRandom3(NULL)
141 {
142  // Constructor
143  AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()");
144 }
145 
146 //________________________________________________________________________
147 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
148  AliAnalysisTaskSE(name),
149  // fOutputFile(NULL),
150  fAnalysisType("AUTOMATIC"),
151  fRPType(RPtype),
152  fCFManager1(NULL),
153  fCFManager2(NULL),
154  fCutsEvent(NULL),
155  fCutsRP(NULL),
156  fCutsPOI(NULL),
157  fCutContainer(new TList()),
158  fQAList(NULL),
159  fMinMult(0),
160  fMaxMult(10000000),
161  fMinA(-1.0),
162  fMaxA(-0.01),
163  fMinB(0.01),
164  fMaxB(1.0),
165  fQAon(on),
166  fLoadCandidates(bCandidates),
167  fNbinsMult(10000),
168  fNbinsPt(100),
169  fNbinsPhi(100),
170  fNbinsEta(200),
171  fNbinsQ(500),
172  fNbinsMass(1),
173  fMultMin(0.),
174  fMultMax(10000.),
175  fPtMin(0.),
176  fPtMax(10.),
177  fPhiMin(0.),
178  fPhiMax(TMath::TwoPi()),
179  fEtaMin(-5.),
180  fEtaMax(5.),
181  fQMin(0.),
182  fQMax(3.),
183  fMassMin(-1.),
184  fMassMax(0.),
185  fHistWeightvsPhiMin(0.),
186  fHistWeightvsPhiMax(3.),
187  fExcludedEtaMin(0.),
188  fExcludedEtaMax(0.),
189  fExcludedPhiMin(0.),
190  fExcludedPhiMax(0.),
191  fAfterburnerOn(kFALSE),
192  fNonFlowNumberOfTrackClones(0),
193  fV1(0.),
194  fV2(0.),
195  fV3(0.),
196  fV4(0.),
197  fV5(0.),
198  fDifferentialV2(0),
199  fFlowEvent(NULL),
200  fShuffleTracks(kFALSE),
201  fMyTRandom3(NULL)
202 {
203  // Constructor
204  AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)");
205  fMyTRandom3 = new TRandom3(iseed);
206  gRandom->SetSeed(fMyTRandom3->Integer(65539));
207 
208  int availableINslot=1;
209  //FMD input slot
210  if (strcmp(RPtype,"FMD")==0) {
211  DefineInput(availableINslot++, TList::Class());
212  }
213  //Candidates input slot
214  if( fLoadCandidates )
215  DefineInput(availableINslot, TObjArray::Class());
216 
217  // Define output slots here
218  // Define here the flow event output
219  DefineOutput(1, AliFlowEventSimple::Class());
220  DefineOutput(2, TList::Class());
221 
222  // and for testing open an output file
223  // fOutputFile = new TFile("FlowEvents.root","RECREATE");
224 
225 }
226 
227 //________________________________________________________________________
229 {
230  //
231  // Destructor
232  //
233  delete fMyTRandom3;
234  delete fFlowEvent;
235  delete fCutsEvent;
236  delete fQAList;
237  if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
238  // objects in the output list are deleted
239  // by the TSelector dtor (I hope)
240 
241 }
242 
243 //________________________________________________________________________
245 {
246  //at the beginning of each new run
247  if (fCutsRP) fCutsRP->SetRunsMuon(fInputHandler); // XZhang 20120604
248  if (fCutsPOI) fCutsPOI->SetRunsMuon(fInputHandler); // XZhang 20120604
249  AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
250  if (!fESD) return;
251 
252  Int_t run = fESD->GetRunNumber();
253  AliInfo(Form("Stariting run #%i",run));
254 }
255 
256 //________________________________________________________________________
258 {
259  // Called at every worker node to initialize
260  AliDebug(2,"AliAnalysisTaskFlowEvent::CreateOutputObjects()");
261 
262  if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
263  {
264  AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
265  exit(1);
266  }
267 
268  //set the common constants
271  cc->SetNbinsPt(fNbinsPt);
272  cc->SetNbinsPhi(fNbinsPhi);
273  cc->SetNbinsEta(fNbinsEta);
274  cc->SetNbinsQ(fNbinsQ);
276  cc->SetMultMin(fMultMin);
277  cc->SetMultMax(fMultMax);
278  cc->SetPtMin(fPtMin);
279  cc->SetPtMax(fPtMax);
280  cc->SetPhiMin(fPhiMin);
281  cc->SetPhiMax(fPhiMax);
282  cc->SetEtaMin(fEtaMin);
283  cc->SetEtaMax(fEtaMax);
284  cc->SetQMin(fQMin);
285  cc->SetQMax(fQMax);
286  cc->SetMassMin(fMassMin);
287  cc->SetMassMax(fMassMax);
290 
291  fFlowEvent = new AliFlowEvent(10000);
292 
293  if (fQAon)
294  {
295  fQAList=new TList();
296  fQAList->SetOwner();
297  fQAList->SetName(Form("%s QA",GetName()));
298  if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
299  if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
300  if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
301  fQAList->Add(new TH1F("event plane angle","event plane angle;angle [rad];",100,0.,TMath::TwoPi())); //3
302  PostData(2,fQAList);
303  }
304 }
305 
306 //________________________________________________________________________
307 
309 {
310  // Main loop
311  // Called for each event
312  //delete fFlowEvent;
313  AliMCEvent* mcEvent = MCEvent(); // from TaskSE
314  AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
315  AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
316  AliMultiplicity* myTracklets = NULL;
317  AliESDPmdTrack* pmdtracks = NULL;//pmd
318 
319  int availableINslot=1;
320 
321  if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
322  {
323  AliError("cuts not set");
324  return;
325  }
326 
327  //DEFAULT - automatically takes care of everything
328  if (fAnalysisType == "AUTOMATIC")
329  {
330  //check event cuts
331  if (InputEvent() && !fCutsEvent->IsSelected(InputEvent(),MCEvent()))
332  return;
333 
334  //first attach all possible information to the cuts
335  fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
336  fCutsPOI->SetEvent( InputEvent(), MCEvent() );
337 
338  //then make the event
340  //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
341 
342  // if (myESD)
344  fFlowEvent->SetCentrality(fCutsEvent->GetCentrality(InputEvent(),mcEvent));
345  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
346  }
347 
348  // Make the FlowEvent for MC input
349  else if (fAnalysisType == "MC")
350  {
351  // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
352  // This handler can return the current MC event
353  if (!(fCFManager1&&fCFManager2))
354  {
355  AliError("ERROR: No pointer to correction framework cuts! ");
356  return;
357  }
358  if (!mcEvent)
359  {
360  AliError("ERROR: Could not retrieve MC event");
361  return;
362  }
363 
364  fCFManager1->SetMCEventInfo(mcEvent);
365  fCFManager2->SetMCEventInfo(mcEvent);
366 
367  AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
368 
369  //check multiplicity
370  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
371  {
372  AliWarning("Event does not pass multiplicity cuts"); return;
373  }
374  //make event
376  }
377 
378  // Make the FlowEvent for ESD input
379  else if (fAnalysisType == "ESD")
380  {
381  if (!(fCFManager1&&fCFManager2))
382  {
383  AliError("ERROR: No pointer to correction framework cuts!");
384  return;
385  }
386  if (!myESD)
387  {
388  AliError("ERROR: ESD not available");
389  return;
390  }
391 
392  //check the offline trigger (check if the event has the correct trigger)
393  AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
394 
395  //check multiplicity
396  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
397  {
398  AliWarning("Event does not pass multiplicity cuts"); return;
399  }
400 
401  //make event
402  if (fRPType == "Global") {
404  }
405  else if (fRPType == "TPCOnly") {
406  fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
407  }
408  else if (fRPType == "TPCHybrid") {
409  fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
410  }
411  else if (fRPType == "Tracklet"){
412  fFlowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
413  }
414  else if (fRPType == "FMD"){
415  TList* dataFMD = dynamic_cast<TList*>(GetInputData(availableINslot++));
416  if(!dataFMD) {
417  cout<<" No dataFMD "<<endl;
418  return;
419  }
420  TH2F* histFMD = dynamic_cast<TH2F*>(dataFMD->FindObject("dNdetadphiHistogramTrVtx"));
421  if (!histFMD) {
422  cout<< "No histFMD"<<endl;
423  return;
424  }
425  fFlowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
426  }
427  else if (fRPType == "PMD"){
428  fFlowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
429  }
430  else return;
431 
432  // if monte carlo event get reaction plane from monte carlo (depends on generator)
433  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
434  //set reference multiplicity, TODO: maybe move it to the constructor?
435  fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
436  }
437 
438  // Make the FlowEvent for ESD input combined with MC info
439  else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
440  {
441  if (!(fCFManager1&&fCFManager2))
442  {
443  AliError("ERROR: No pointer to correction framework cuts! ");
444  return;
445  }
446  if (!myESD)
447  {
448  AliError("ERROR: ESD not available");
449  return;
450  }
451  AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
452 
453  if (!mcEvent)
454  {
455  AliError("ERROR: Could not retrieve MC event");
456  return;
457  }
458 
459  fCFManager1->SetMCEventInfo(mcEvent);
460  fCFManager2->SetMCEventInfo(mcEvent);
461 
462  //check multiplicity
463  if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
464  {
465  AliWarning("Event does not pass multiplicity cuts"); return;
466  }
467 
468  //make event
469  if (fAnalysisType == "ESDMCkineESD")
470  {
472  }
473  else if (fAnalysisType == "ESDMCkineMC")
474  {
476  }
477  if (!fFlowEvent) return;
478  // if monte carlo event get reaction plane from monte carlo (depends on generator)
479  if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
480  //set reference multiplicity, TODO: maybe move it to the constructor?
481  fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
482  }
483  // Make the FlowEventSimple for AOD input
484  else if (fAnalysisType == "AOD")
485  {
486  if (!myAOD)
487  {
488  AliError("ERROR: AOD not available");
489  return;
490  }
491  AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
492  fFlowEvent = new AliFlowEvent(myAOD);
493  }
494 
495  //inject candidates
496  if(fLoadCandidates) {
497  TObjArray* candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
498  //if(candidates->GetEntriesFast())
499  // printf("I received %d candidates\n",candidates->GetEntriesFast());
500  if (candidates)
501  {
502  for(int iCand=0; iCand!=candidates->GetEntriesFast(); ++iCand ) {
503  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(candidates->At(iCand));
504  if (!cand) continue;
505  //printf(" - Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
506  for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
507  //printf(" - Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
508  for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs ) {
509  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
510  if (!iRP) continue;
511  if( !iRP->InRPSelection() )
512  continue;
513  if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
514  //printf(" was in RP set");
515  //cand->SetDaughter( iDau, iRP );
516  //temporarily untagging all daugters
517  iRP->SetForRPSelection(kFALSE);
519  }
520  }
521  //printf("\n");
522  }
523  cand->SetForPOISelection(kTRUE);
524  fFlowEvent->InsertTrack( ((AliFlowTrack*) cand) );
525  }
526  }
527  }
528 
529  if (!fFlowEvent) return; //shuts up coverity
530 
531  //check final event cuts
532  Int_t mult = fFlowEvent->NumberOfTracks();
533  // AliInfo(Form("FlowEvent has %i tracks",mult));
534  if (mult<fMinMult || mult>fMaxMult)
535  {
536  AliWarning("FlowEvent cut on multiplicity"); return;
537  }
538 
539  //define dead zone
541 
542 
545  if (fAfterburnerOn)
546  {
547  //if reaction plane not set from elsewhere randomize it before adding flow
549  fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
550 
551  if(fDifferentialV2)
553  else
554  fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
555  fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
556  }
558 
559  //tag subEvents
561 
562  //QA
563  if (fQAon)
564  {
565  TH1* h1 = static_cast<TH1*>(fQAList->FindObject("event plane angle"));
566  h1->Fill(fFlowEvent->GetMCReactionPlaneAngle());
567  }
568 
569  //do we want to serve shullfed tracks to everybody?
571 
572  // associate the mother particles to their daughters in the flow event (if any)
574 
575  //fListHistos->Print();
576  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
577  PostData(1,fFlowEvent);
578 }
579 
580 //________________________________________________________________________
582 {
583  // Called once at the end of the query -- do not call in case of CAF
584 }
585 
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
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 SetRunsMuon(const AliInputEventHandler *eventHandler)
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)
ClassImp(AliAnalysisTaskFlowEvent) AliAnalysisTaskFlowEvent
void InsertTrack(AliFlowTrack *)
Definition: External.C:196
Int_t NumberOfTracks() const