53a1666019b3e656f9ebac57debe4ba4d3df74f5
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliAnalysisTaskFlowEvent.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 // 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
26 ////////////////////////////////////////////////////
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 "TH2F.h"
34 #include "TRandom3.h"
35 #include "TTimeStamp.h"
36
37 // ALICE Analysis Framework
38 #include "AliAnalysisManager.h"
39 #include "AliAnalysisTaskSE.h"
40 #include "AliESDtrackCuts.h"
41
42 // ESD interface
43 #include "AliESDEvent.h"
44 #include "AliESDInputHandler.h"
45
46 // AOD interface
47 #include "AliAODEvent.h"
48 #include "AliAODInputHandler.h"
49
50 // Monte Carlo Event
51 #include "AliMCEventHandler.h"
52 #include "AliMCEvent.h"
53
54 // ALICE Correction Framework
55 #include "AliCFManager.h"
56
57 // Interface to Event generators to get Reaction Plane Angle
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliGenHijingEventHeader.h"
60 #include "AliGenGeVSimEventHeader.h"
61 #include "AliGenEposEventHeader.h"
62
63 // Interface to Load short life particles
64 #include "TObjArray.h"
65 #include "AliFlowCandidateTrack.h"
66
67 // Interface to make the Flow Event Simple used in the flow analysis methods
68 #include "AliFlowEvent.h"
69 #include "AliFlowTrackCuts.h"
70 #include "AliFlowEventCuts.h"
71 #include "AliFlowCommonConstants.h"
72 #include "AliAnalysisTaskFlowEvent.h"
73
74 #include "AliLog.h"
75
76 ClassImp(AliAnalysisTaskFlowEvent)
77
78 //________________________________________________________________________
79 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
80   AliAnalysisTaskSE(),
81   //  fOutputFile(NULL),
82   fAnalysisType("AUTOMATIC"),
83   fRPType("Global"),
84   fCFManager1(NULL),
85   fCFManager2(NULL),
86   fCutsEvent(NULL),
87   fCutsRP(NULL),
88   fCutsPOI(NULL),
89   fQAInt(NULL),
90   fQADiff(NULL),
91   fMinMult(0),
92   fMaxMult(10000000),
93   fMinA(-1.0),
94   fMaxA(-0.01),
95   fMinB(0.01),
96   fMaxB(1.0),
97   fQA(kFALSE),
98   fLoadCandidates(kFALSE),
99   fNbinsMult(10000),
100   fNbinsPt(100),   
101   fNbinsPhi(100),
102   fNbinsEta(200),
103   fNbinsQ(500),
104   fMultMin(0.),            
105   fMultMax(10000.),
106   fPtMin(0.),        
107   fPtMax(10.),
108   fPhiMin(0.),       
109   fPhiMax(TMath::TwoPi()),
110   fEtaMin(-5.),      
111   fEtaMax(5.),       
112   fQMin(0.),         
113   fQMax(3.),
114   fExcludedEtaMin(0.), 
115   fExcludedEtaMax(0.), 
116   fExcludedPhiMin(0.),
117   fExcludedPhiMax(0.),
118   fAfterburnerOn(kFALSE),
119   fNonFlowNumberOfTrackClones(0),
120   fV1(0.),
121   fV2(0.),
122   fV3(0.),
123   fV4(0.),
124   fMyTRandom3(NULL)
125 {
126   // Constructor
127   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
128 }
129
130 //________________________________________________________________________
131 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
132   AliAnalysisTaskSE(name),
133   //  fOutputFile(NULL),
134   fAnalysisType("AUTOMATIC"),
135   fRPType(RPtype),
136   fCFManager1(NULL),
137   fCFManager2(NULL),
138   fCutsEvent(NULL),
139   fCutsRP(NULL),
140   fCutsPOI(NULL),
141   fQAInt(NULL),
142   fQADiff(NULL),
143   fMinMult(0),
144   fMaxMult(10000000),
145   fMinA(-1.0),
146   fMaxA(-0.01),
147   fMinB(0.01),
148   fMaxB(1.0),
149   fQA(on),
150   fLoadCandidates(bCandidates),
151   fNbinsMult(10000),
152   fNbinsPt(100),   
153   fNbinsPhi(100),
154   fNbinsEta(200),
155   fNbinsQ(500),
156   fMultMin(0.),            
157   fMultMax(10000.),
158   fPtMin(0.),        
159   fPtMax(10.),
160   fPhiMin(0.),       
161   fPhiMax(TMath::TwoPi()),
162   fEtaMin(-5.),      
163   fEtaMax(5.),       
164   fQMin(0.),         
165   fQMax(3.),
166   fExcludedEtaMin(0.), 
167   fExcludedEtaMax(0.), 
168   fExcludedPhiMin(0.),
169   fExcludedPhiMax(0.),
170   fAfterburnerOn(kFALSE),
171   fNonFlowNumberOfTrackClones(0),
172   fV1(0.),
173   fV2(0.),
174   fV3(0.),
175   fV4(0.),
176   fMyTRandom3(NULL)
177 {
178   // Constructor
179   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
180   fMyTRandom3 = new TRandom3(iseed);
181   gRandom->SetSeed(fMyTRandom3->Integer(65539));
182
183   int availableINslot=1;
184   //FMD input slot
185   if (strcmp(RPtype,"FMD")==0) {
186     DefineInput(availableINslot++, TList::Class());
187   }
188   //Candidates input slot
189   if( fLoadCandidates )
190     DefineInput(availableINslot, TObjArray::Class());
191
192   // Define output slots here
193   // Define here the flow event output
194   DefineOutput(1, AliFlowEventSimple::Class());
195   if(on)
196   {
197     DefineOutput(2, TList::Class());
198     DefineOutput(3, TList::Class());
199   }
200
201   // and for testing open an output file
202   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
203
204 }
205
206 //________________________________________________________________________
207 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
208 {
209   //
210   // Destructor
211   //
212   delete fMyTRandom3;
213   // objects in the output list are deleted
214   // by the TSelector dtor (I hope)
215
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskFlowEvent::NotifyRun()
220 {
221   //at the beginning of each new run
222         AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
223   if (!fESD) return;
224
225   Int_t run = fESD->GetRunNumber();  
226   AliInfo(Form("Stariting run #%i",run));
227 }
228
229 //________________________________________________________________________
230 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
231 {
232   // Called at every worker node to initialize
233   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
234
235   if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
236   {
237     AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
238     exit(1);
239   }
240
241   //set the common constants
242   AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
243   cc->SetNbinsMult(fNbinsMult);
244   cc->SetNbinsPt(fNbinsPt);
245   cc->SetNbinsPhi(fNbinsPhi); 
246   cc->SetNbinsEta(fNbinsEta);
247   cc->SetNbinsQ(fNbinsQ);
248   cc->SetMultMin(fMultMin);
249   cc->SetMultMax(fMultMax);
250   cc->SetPtMin(fPtMin);
251   cc->SetPtMax(fPtMax);
252   cc->SetPhiMin(fPhiMin);
253   cc->SetPhiMax(fPhiMax);
254   cc->SetEtaMin(fEtaMin);
255   cc->SetEtaMax(fEtaMax);
256   cc->SetQMin(fQMin);
257   cc->SetQMax(fQMax);
258
259 }
260
261 //________________________________________________________________________
262 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
263 {
264   // Main loop
265   // Called for each event
266   AliFlowEvent* flowEvent = NULL;
267   AliMCEvent*  mcEvent = MCEvent();                              // from TaskSE
268   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
269   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
270   AliMultiplicity* myTracklets = NULL;
271   AliESDPmdTrack* pmdtracks = NULL;//pmd      
272   TH2F* histFMD = NULL;
273
274   int availableINslot=1;
275   if(strcmp(fRPType,"FMD")==0) {
276     TList* FMDdata = dynamic_cast<TList*>(GetInputData(availableINslot++));
277     if(!FMDdata) {
278       cout<<" No FMDdata "<<endl;
279       exit(2);
280     }
281     histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
282     if (!histFMD) {
283       cout<< "No histFMD"<<endl;
284       exit(2);
285     }
286   }
287
288   if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
289   {
290     AliError("cuts not set");
291     return;
292   }
293
294   //use the new and temporarily incomplete way of doing things
295   if (fAnalysisType == "AUTOMATIC")
296   {
297     //check event cuts
298     if (!fCutsEvent->IsSelected(InputEvent())) return;
299
300     //first attach all possible information to the cuts
301     fCutsRP->SetEvent( InputEvent(), MCEvent() );  //attach event
302     fCutsPOI->SetEvent( InputEvent(), MCEvent() );
303
304     //then make the event
305     flowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
306     if (myESD)
307       flowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent()));
308     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
309   }
310
311   // Make the FlowEvent for MC input
312   if (fAnalysisType == "MC")
313   {
314     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
315     // This handler can return the current MC event
316     if (!(fCFManager1&&fCFManager2))
317     {
318       AliError("ERROR: No pointer to correction framework cuts! ");
319       return;
320     }
321     if (!mcEvent)
322     {
323       AliError("ERROR: Could not retrieve MC event");
324       return;
325     }
326
327     fCFManager1->SetMCEventInfo(mcEvent);
328     fCFManager2->SetMCEventInfo(mcEvent);
329     
330     AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
331
332     //check multiplicity 
333     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
334     {
335       AliWarning("Event does not pass multiplicity cuts"); return;
336     }
337     //make event
338     flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
339   }
340
341   // Make the FlowEvent for ESD input
342   else if (fAnalysisType == "ESD")
343   {
344     if (!(fCFManager1&&fCFManager2))
345     {
346       AliError("ERROR: No pointer to correction framework cuts!");
347       return;
348     }
349     if (!myESD)
350     {
351       AliError("ERROR: ESD not available");
352       return;
353     }
354
355     //check the offline trigger (check if the event has the correct trigger)
356     AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
357
358     //check multiplicity
359     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
360     {
361       AliWarning("Event does not pass multiplicity cuts"); return;
362     }
363
364     //make event
365     if (fRPType == "Global") {
366       flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
367     }
368     if (fRPType == "TPCOnly") {
369       flowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
370     }
371     if (fRPType == "TPCHybrid") {
372       flowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
373     }
374     else if (fRPType == "Tracklet"){
375       flowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
376     }
377     else if (fRPType == "FMD"){
378       flowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
379     }
380     //pmd
381     else if (fRPType == "PMD"){
382       flowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
383     }
384     //pmd
385     
386     // if monte carlo event get reaction plane from monte carlo (depends on generator)
387     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
388     //set reference multiplicity, TODO: maybe move it to the constructor?
389     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
390   }
391
392   // Make the FlowEvent for ESD input combined with MC info
393   else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
394   {
395     if (!(fCFManager1&&fCFManager2))
396     {
397       AliError("ERROR: No pointer to correction framework cuts! ");
398       return;
399     }
400     if (!myESD)
401     {
402       AliError("ERROR: ESD not available");
403       return;
404     }
405     AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
406
407     if (!mcEvent)
408     {
409       AliError("ERROR: Could not retrieve MC event");
410       return;
411     }
412
413     fCFManager1->SetMCEventInfo(mcEvent);
414     fCFManager2->SetMCEventInfo(mcEvent);
415
416     //check multiplicity
417     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
418     {
419       AliWarning("Event does not pass multiplicity cuts"); return;
420     }
421
422     //make event
423     if (fAnalysisType == "ESDMCkineESD")
424     {
425       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
426     }
427     else if (fAnalysisType == "ESDMCkineMC")
428     {
429       flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
430     }
431     // if monte carlo event get reaction plane from monte carlo (depends on generator)
432     if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
433     //set reference multiplicity, TODO: maybe move it to the constructor?
434     flowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
435   }
436   // Make the FlowEventSimple for AOD input
437   else if (fAnalysisType == "AOD")
438   {
439     if (!myAOD)
440     {
441       AliError("ERROR: AOD not available");
442       return;
443     }
444     AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
445     flowEvent = new AliFlowEvent(myAOD);
446   }
447
448   //inject candidates
449   if(fLoadCandidates) {
450     TObjArray* Candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
451     AliFlowCandidateTrack *cand;
452     for(int iCand=0; iCand!=Candidates->GetEntriesFast(); ++iCand ) {
453       cand = dynamic_cast<AliFlowCandidateTrack*>(Candidates->At(iCand));
454       cand->SetForPOISelection();
455       flowEvent->AddTrack(cand);
456     }
457   }
458
459   //check final event cuts
460   Int_t mult = flowEvent->NumberOfTracks();
461   AliInfo(Form("FlowEvent has %i tracks",mult));
462   if (mult<fMinMult || mult>fMaxMult)
463   {
464     AliWarning("FlowEvent cut on multiplicity"); return;
465   }
466
467   //define dead zone
468   flowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
469
470
471   //////////////////////////////////////////////////////////////////////////////
472   ///////////////////////////AFTERBURNER
473   if (fAfterburnerOn)
474   {
475     //if reaction plane not set from elsewhere randomize it before adding flow
476     if (!flowEvent->IsSetMCReactionPlaneAngle())
477       flowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
478
479     flowEvent->AddFlow(fV1,fV2,fV3,fV4);     //add flow
480     flowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
481   }
482   //////////////////////////////////////////////////////////////////////////////
483
484   //tag subEvents
485   flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
486
487   //fListHistos->Print();
488   //fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
489   PostData(1,flowEvent);
490   if (fQA)
491   {
492     PostData(2,fQAInt);
493     PostData(3,fQADiff);
494   }
495 }
496
497 //________________________________________________________________________
498 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
499 {
500   // Called once at the end of the query -- do not call in case of CAF
501 }
502