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 // AliAnalysisTaskFlowEvent:
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 ////////////////////////////////////////////////////
28 #include "Riostream.h" //needed as include
31 #include "TFile.h" //needed as include
36 #include "TTimeStamp.h"
38 // ALICE Analysis Framework
39 #include "AliAnalysisManager.h"
40 #include "AliAnalysisTaskSE.h"
41 #include "AliESDtrackCuts.h"
44 #include "AliESDEvent.h"
45 #include "AliESDInputHandler.h"
48 #include "AliAODEvent.h"
49 #include "AliAODInputHandler.h"
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
55 // ALICE Correction Framework
56 #include "AliCFManager.h"
58 // Interface to Event generators to get Reaction Plane Angle
59 #include "AliGenCocktailEventHeader.h"
60 #include "AliGenHijingEventHeader.h"
61 #include "AliGenGeVSimEventHeader.h"
62 #include "AliGenEposEventHeader.h"
64 // Interface to Load short life particles
65 #include "TObjArray.h"
66 #include "AliFlowCandidateTrack.h"
68 // Interface to make the Flow Event Simple used in the flow analysis methods
69 #include "AliFlowEvent.h"
70 #include "AliFlowTrackCuts.h"
71 #include "AliFlowEventCuts.h"
72 #include "AliFlowCommonConstants.h"
73 #include "AliAnalysisTaskFlowEvent.h"
79 ClassImp(AliAnalysisTaskFlowEvent)
81 //________________________________________________________________________
82 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
85 fAnalysisType("AUTOMATIC"),
101 fLoadCandidates(kFALSE),
113 fPhiMax(TMath::TwoPi()),
120 fHistWeightvsPhiMin(0.),
121 fHistWeightvsPhiMax(3.),
126 fAfterburnerOn(kFALSE),
127 fNonFlowNumberOfTrackClones(0),
138 AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()");
141 //________________________________________________________________________
142 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
143 AliAnalysisTaskSE(name),
144 // fOutputFile(NULL),
145 fAnalysisType("AUTOMATIC"),
152 fCutContainer(new TList()),
161 fLoadCandidates(bCandidates),
173 fPhiMax(TMath::TwoPi()),
180 fHistWeightvsPhiMin(0.),
181 fHistWeightvsPhiMax(3.),
186 fAfterburnerOn(kFALSE),
187 fNonFlowNumberOfTrackClones(0),
198 AliDebug(2,"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)");
199 fMyTRandom3 = new TRandom3(iseed);
200 gRandom->SetSeed(fMyTRandom3->Integer(65539));
202 int availableINslot=1;
204 if (strcmp(RPtype,"FMD")==0) {
205 DefineInput(availableINslot++, TList::Class());
207 //Candidates input slot
208 if( fLoadCandidates )
209 DefineInput(availableINslot, TObjArray::Class());
211 // Define output slots here
212 // Define here the flow event output
213 DefineOutput(1, AliFlowEventSimple::Class());
214 DefineOutput(2, TList::Class());
216 // and for testing open an output file
217 // fOutputFile = new TFile("FlowEvents.root","RECREATE");
221 //________________________________________________________________________
222 AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
231 if (fCutContainer) fCutContainer->Delete(); delete fCutContainer;
232 // objects in the output list are deleted
233 // by the TSelector dtor (I hope)
237 //________________________________________________________________________
238 void AliAnalysisTaskFlowEvent::NotifyRun()
240 //at the beginning of each new run
241 if (fCutsRP) fCutsRP->SetRunsMuon(fInputHandler); // XZhang 20120604
242 if (fCutsPOI) fCutsPOI->SetRunsMuon(fInputHandler); // XZhang 20120604
243 AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
246 Int_t run = fESD->GetRunNumber();
247 AliInfo(Form("Stariting run #%i",run));
250 //________________________________________________________________________
251 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
253 // Called at every worker node to initialize
254 AliDebug(2,"AliAnalysisTaskFlowEvent::CreateOutputObjects()");
256 if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
258 AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
262 //set the common constants
263 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
264 cc->SetNbinsMult(fNbinsMult);
265 cc->SetNbinsPt(fNbinsPt);
266 cc->SetNbinsPhi(fNbinsPhi);
267 cc->SetNbinsEta(fNbinsEta);
268 cc->SetNbinsQ(fNbinsQ);
269 cc->SetNbinsMass(fNbinsMass);
270 cc->SetMultMin(fMultMin);
271 cc->SetMultMax(fMultMax);
272 cc->SetPtMin(fPtMin);
273 cc->SetPtMax(fPtMax);
274 cc->SetPhiMin(fPhiMin);
275 cc->SetPhiMax(fPhiMax);
276 cc->SetEtaMin(fEtaMin);
277 cc->SetEtaMax(fEtaMax);
280 cc->SetMassMin(fMassMin);
281 cc->SetMassMax(fMassMax);
282 cc->SetHistWeightvsPhiMax(fHistWeightvsPhiMax);
283 cc->SetHistWeightvsPhiMin(fHistWeightvsPhiMin);
285 fFlowEvent = new AliFlowEvent(3000);
291 fQAList->SetName(Form("%s QA",GetName()));
292 if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA()); //0
293 if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA()); //1
294 if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA()); //2
295 fQAList->Add(new TH1F("event plane angle","event plane angle;angle [rad];",100,0.,TMath::TwoPi())); //3
300 //________________________________________________________________________
301 void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
304 // Called for each event
306 AliMCEvent* mcEvent = MCEvent(); // from TaskSE
307 AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
308 AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
309 AliMultiplicity* myTracklets = NULL;
310 AliESDPmdTrack* pmdtracks = NULL;//pmd
312 int availableINslot=1;
314 if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
316 AliError("cuts not set");
320 //DEFAULT - automatically takes care of everything
321 if (fAnalysisType == "AUTOMATIC")
324 if (InputEvent() && !fCutsEvent->IsSelected(InputEvent())) return;
326 //first attach all possible information to the cuts
327 fCutsRP->SetEvent( InputEvent(), MCEvent() ); //attach event
328 fCutsPOI->SetEvent( InputEvent(), MCEvent() );
330 //then make the event
331 fFlowEvent->Fill( fCutsRP, fCutsPOI );
332 //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
335 fFlowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent()));
336 if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
339 // Make the FlowEvent for MC input
340 else if (fAnalysisType == "MC")
342 // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
343 // This handler can return the current MC event
344 if (!(fCFManager1&&fCFManager2))
346 AliError("ERROR: No pointer to correction framework cuts! ");
351 AliError("ERROR: Could not retrieve MC event");
355 fCFManager1->SetMCEventInfo(mcEvent);
356 fCFManager2->SetMCEventInfo(mcEvent);
358 AliInfo(Form("Number of MC particles: %d", mcEvent->GetNumberOfTracks()));
361 if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtGenCuts,mcEvent))
363 AliWarning("Event does not pass multiplicity cuts"); return;
366 fFlowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
369 // Make the FlowEvent for ESD input
370 else if (fAnalysisType == "ESD")
372 if (!(fCFManager1&&fCFManager2))
374 AliError("ERROR: No pointer to correction framework cuts!");
379 AliError("ERROR: ESD not available");
383 //check the offline trigger (check if the event has the correct trigger)
384 AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
387 if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
389 AliWarning("Event does not pass multiplicity cuts"); return;
393 if (fRPType == "Global") {
394 fFlowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
396 else if (fRPType == "TPCOnly") {
397 fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
399 else if (fRPType == "TPCHybrid") {
400 fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
402 else if (fRPType == "Tracklet"){
403 fFlowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
405 else if (fRPType == "FMD"){
406 TList* dataFMD = dynamic_cast<TList*>(GetInputData(availableINslot++));
408 cout<<" No dataFMD "<<endl;
411 TH2F* histFMD = dynamic_cast<TH2F*>(dataFMD->FindObject("dNdetadphiHistogramTrVtx"));
413 cout<< "No histFMD"<<endl;
416 fFlowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
418 else if (fRPType == "PMD"){
419 fFlowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
423 // if monte carlo event get reaction plane from monte carlo (depends on generator)
424 if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
425 //set reference multiplicity, TODO: maybe move it to the constructor?
426 fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
429 // Make the FlowEvent for ESD input combined with MC info
430 else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
432 if (!(fCFManager1&&fCFManager2))
434 AliError("ERROR: No pointer to correction framework cuts! ");
439 AliError("ERROR: ESD not available");
442 AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
446 AliError("ERROR: Could not retrieve MC event");
450 fCFManager1->SetMCEventInfo(mcEvent);
451 fCFManager2->SetMCEventInfo(mcEvent);
454 if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
456 AliWarning("Event does not pass multiplicity cuts"); return;
460 if (fAnalysisType == "ESDMCkineESD")
462 fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
464 else if (fAnalysisType == "ESDMCkineMC")
466 fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
468 if (!fFlowEvent) return;
469 // if monte carlo event get reaction plane from monte carlo (depends on generator)
470 if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
471 //set reference multiplicity, TODO: maybe move it to the constructor?
472 fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
474 // Make the FlowEventSimple for AOD input
475 else if (fAnalysisType == "AOD")
479 AliError("ERROR: AOD not available");
482 AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
483 fFlowEvent = new AliFlowEvent(myAOD);
487 if(fLoadCandidates) {
488 TObjArray* candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
489 //if(candidates->GetEntriesFast())
490 // printf("I received %d candidates\n",candidates->GetEntriesFast());
493 for(int iCand=0; iCand!=candidates->GetEntriesFast(); ++iCand ) {
494 AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(candidates->At(iCand));
496 //printf(" - Checking at candidate %d with %d daughters: mass %f\n",iCand,cand->GetNDaughters(),cand->Mass());
497 for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
498 //printf(" - Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
499 for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs ) {
500 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
502 if( !iRP->InRPSelection() )
504 if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
505 //printf(" was in RP set");
506 //cand->SetDaughter( iDau, iRP );
507 //temporarily untagging all daugters
508 iRP->SetForRPSelection(kFALSE);
509 fFlowEvent->SetNumberOfRPs( fFlowEvent->GetNumberOfRPs() -1 );
514 cand->SetForPOISelection(kTRUE);
515 fFlowEvent->InsertTrack( ((AliFlowTrack*) cand) );
520 if (!fFlowEvent) return; //shuts up coverity
522 //check final event cuts
523 Int_t mult = fFlowEvent->NumberOfTracks();
524 // AliInfo(Form("FlowEvent has %i tracks",mult));
525 if (mult<fMinMult || mult>fMaxMult)
527 AliWarning("FlowEvent cut on multiplicity"); return;
531 fFlowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
534 //////////////////////////////////////////////////////////////////////////////
535 ///////////////////////////AFTERBURNER
538 //if reaction plane not set from elsewhere randomize it before adding flow
539 if (!fFlowEvent->IsSetMCReactionPlaneAngle())
540 fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
543 fFlowEvent->AddV2(fDifferentialV2);
545 fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5); //add flow
546 fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
548 //////////////////////////////////////////////////////////////////////////////
551 fFlowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
556 TH1* h1 = static_cast<TH1*>(fQAList->FindObject("event plane angle"));
557 h1->Fill(fFlowEvent->GetMCReactionPlaneAngle());
560 //fListHistos->Print();
561 //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
562 PostData(1,fFlowEvent);
565 //________________________________________________________________________
566 void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
568 // Called once at the end of the query -- do not call in case of CAF