]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/FLOW/AliFlowTasks/AliAnalysisTaskFlowEvent.cxx
added some QA to flow track cuts and event cuts, trackID fix, better redoFinish and...
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowTasks / AliAnalysisTaskFlowEvent.cxx
index 850c6b49f56c0dc291314859de5eea34e175465a..eb5d4d7a023c6b251361b1644de164b2f1d76078 100644 (file)
 #include "TTree.h"
 #include "TFile.h" //needed as include
 #include "TList.h"
+#include "TH2F.h"
 #include "TRandom3.h"
 #include "TTimeStamp.h"
 
 // ALICE Analysis Framework
 #include "AliAnalysisManager.h"
 #include "AliAnalysisTaskSE.h"
+#include "AliESDtrackCuts.h"
 
 // ESD interface
 #include "AliESDEvent.h"
 #include "AliGenGeVSimEventHeader.h"
 #include "AliGenEposEventHeader.h"
 
+// Interface to Load short life particles
+#include "TObjArray.h"
+#include "AliFlowCandidateTrack.h"
+
 // Interface to make the Flow Event Simple used in the flow analysis methods
 #include "AliFlowEvent.h"
-
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventCuts.h"
+#include "AliFlowCommonConstants.h"
 #include "AliAnalysisTaskFlowEvent.h"
 
 #include "AliLog.h"
@@ -71,72 +79,123 @@ ClassImp(AliAnalysisTaskFlowEvent)
 AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent() :
   AliAnalysisTaskSE(),
   //  fOutputFile(NULL),
-  fAnalysisType("ESD"),
+  fAnalysisType("AUTOMATIC"),
+  fRPType("Global"),
   fCFManager1(NULL),
   fCFManager2(NULL),
-  fQAInt(NULL),
-  fQADiff(NULL),
+  fCutsEvent(NULL),
+  fCutsRP(NULL),
+  fCutsPOI(NULL),
+  fQAList(NULL),
   fMinMult(0),
   fMaxMult(10000000),
   fMinA(-1.0),
   fMaxA(-0.01),
   fMinB(0.01),
   fMaxB(1.0),
-  fQA(kFALSE),
-  fMCReactionPlaneAngle(0.),
-  fCount(0),
-  fNoOfLoops(1),
-  fEllipticFlowValue(0.),
-  fSigmaEllipticFlowValue(0.),
-  fMultiplicityOfEvent(1000000000),
-  fSigmaMultiplicityOfEvent(0),
-  fMyTRandom3(NULL),
-  fbAfterburnerOn(kFALSE)
+  fQAon(kFALSE),
+  fLoadCandidates(kFALSE),
+  fNbinsMult(10000),
+  fNbinsPt(100),   
+  fNbinsPhi(100),
+  fNbinsEta(200),
+  fNbinsQ(500),
+  fMultMin(0.),            
+  fMultMax(10000.),
+  fPtMin(0.),       
+  fPtMax(10.),
+  fPhiMin(0.),      
+  fPhiMax(TMath::TwoPi()),
+  fEtaMin(-5.),             
+  fEtaMax(5.),      
+  fQMin(0.),        
+  fQMax(3.),
+  fExcludedEtaMin(0.), 
+  fExcludedEtaMax(0.), 
+  fExcludedPhiMin(0.),
+  fExcludedPhiMax(0.),
+  fAfterburnerOn(kFALSE),
+  fNonFlowNumberOfTrackClones(0),
+  fV1(0.),
+  fV2(0.),
+  fV3(0.),
+  fV4(0.),
+  fV5(0.),
+  fFlowEvent(NULL),
+  fMyTRandom3(NULL)
 {
   // Constructor
   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent()"<<endl;
 }
 
 //________________________________________________________________________
-AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed) :
+AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, TString RPtype, Bool_t on, UInt_t iseed, Bool_t bCandidates) :
   AliAnalysisTaskSE(name),
   //  fOutputFile(NULL),
-  fAnalysisType("ESD"),
+  fAnalysisType("AUTOMATIC"),
+  fRPType(RPtype),
   fCFManager1(NULL),
   fCFManager2(NULL),
-  fQAInt(NULL),
-  fQADiff(NULL),
+  fCutsEvent(NULL),
+  fCutsRP(NULL),
+  fCutsPOI(NULL),
+  fQAList(NULL),
   fMinMult(0),
   fMaxMult(10000000),
   fMinA(-1.0),
   fMaxA(-0.01),
   fMinB(0.01),
   fMaxB(1.0),
-  fQA(on),
-  fMCReactionPlaneAngle(0.),
-  fCount(0),
-  fNoOfLoops(1),
-  fEllipticFlowValue(0.),
-  fSigmaEllipticFlowValue(0.),
-  fMultiplicityOfEvent(1000000000),
-  fSigmaMultiplicityOfEvent(0),
-  fMyTRandom3(NULL),
-  fbAfterburnerOn(kFALSE)
+  fQAon(on),
+  fLoadCandidates(bCandidates),
+  fNbinsMult(10000),
+  fNbinsPt(100),   
+  fNbinsPhi(100),
+  fNbinsEta(200),
+  fNbinsQ(500),
+  fMultMin(0.),            
+  fMultMax(10000.),
+  fPtMin(0.),       
+  fPtMax(10.),
+  fPhiMin(0.),      
+  fPhiMax(TMath::TwoPi()),
+  fEtaMin(-5.),             
+  fEtaMax(5.),      
+  fQMin(0.),        
+  fQMax(3.),
+  fExcludedEtaMin(0.), 
+  fExcludedEtaMax(0.), 
+  fExcludedPhiMin(0.),
+  fExcludedPhiMax(0.),
+  fAfterburnerOn(kFALSE),
+  fNonFlowNumberOfTrackClones(0),
+  fV1(0.),
+  fV2(0.),
+  fV3(0.),
+  fV4(0.),
+  fV5(0.),
+  fFlowEvent(NULL),
+  fMyTRandom3(NULL)
 {
   // Constructor
   cout<<"AliAnalysisTaskFlowEvent::AliAnalysisTaskFlowEvent(const char *name, Bool_t on, UInt_t iseed)"<<endl;
   fMyTRandom3 = new TRandom3(iseed);
   gRandom->SetSeed(fMyTRandom3->Integer(65539));
 
+  int availableINslot=1;
+  //FMD input slot
+  if (strcmp(RPtype,"FMD")==0) {
+    DefineInput(availableINslot++, TList::Class());
+  }
+  //Candidates input slot
+  if( fLoadCandidates )
+    DefineInput(availableINslot, TObjArray::Class());
 
   // Define output slots here
   // Define here the flow event output
   DefineOutput(1, AliFlowEventSimple::Class());
-  if(on)
-  {
-    DefineOutput(2, TList::Class());
-    DefineOutput(3, TList::Class());
-  }
+  DefineOutput(2, TList::Class());
+
   // and for testing open an output file
   //  fOutputFile = new TFile("FlowEvents.root","RECREATE");
 
@@ -148,23 +207,69 @@ AliAnalysisTaskFlowEvent::~AliAnalysisTaskFlowEvent()
   //
   // Destructor
   //
-  if (fMyTRandom3) delete fMyTRandom3;
+  delete fMyTRandom3;
+  delete fFlowEvent;
+  delete fCutsEvent;
+  delete fCutsRP;
+  delete fCutsPOI;
+  delete fQAList;
   // objects in the output list are deleted
   // by the TSelector dtor (I hope)
 
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskFlowEvent::NotifyRun()
+{
+  //at the beginning of each new run
+       AliESDEvent* fESD = dynamic_cast<AliESDEvent*> (InputEvent());
+  if (!fESD) return;
+
+  Int_t run = fESD->GetRunNumber();  
+  AliInfo(Form("Stariting run #%i",run));
+}
+
 //________________________________________________________________________
 void AliAnalysisTaskFlowEvent::UserCreateOutputObjects()
 {
   // Called at every worker node to initialize
   cout<<"AliAnalysisTaskFlowEvent::CreateOutputObjects()"<<endl;
 
-  if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC"))
+  if (!(fAnalysisType == "AOD" || fAnalysisType == "ESD" || fAnalysisType == "ESDMCkineESD"  || fAnalysisType == "ESDMCkineMC" || fAnalysisType == "MC" || fAnalysisType == "AUTOMATIC"))
   {
-    AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD and MC are allowed.");
+    AliError("WRONG ANALYSIS TYPE! only ESD, ESDMCkineESD, ESDMCkineMC, AOD, MC and AUTOMATIC are allowed.");
     exit(1);
   }
+
+  //set the common constants
+  AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
+  cc->SetNbinsMult(fNbinsMult);
+  cc->SetNbinsPt(fNbinsPt);
+  cc->SetNbinsPhi(fNbinsPhi); 
+  cc->SetNbinsEta(fNbinsEta);
+  cc->SetNbinsQ(fNbinsQ);
+  cc->SetMultMin(fMultMin);
+  cc->SetMultMax(fMultMax);
+  cc->SetPtMin(fPtMin);
+  cc->SetPtMax(fPtMax);
+  cc->SetPhiMin(fPhiMin);
+  cc->SetPhiMax(fPhiMax);
+  cc->SetEtaMin(fEtaMin);
+  cc->SetEtaMax(fEtaMax);
+  cc->SetQMin(fQMin);
+  cc->SetQMax(fQMax);
+
+  fFlowEvent = new AliFlowEvent(3000);
+
+  if (fQAon)
+  {
+    fQAList=new TList();
+    fQAList->SetName(Form("%s QA",GetName()));
+    if (fCutsEvent->GetQA()) fQAList->Add(fCutsEvent->GetQA());
+    if (fCutsRP->GetQA()) fQAList->Add(fCutsRP->GetQA());
+    if (fCutsPOI->GetQA())fQAList->Add(fCutsPOI->GetQA());
+    PostData(2,fQAList);
+  }
 }
 
 //________________________________________________________________________
@@ -172,13 +277,42 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
 {
   // Main loop
   // Called for each event
-  AliFlowEvent* flowEvent = NULL;
+  //delete fFlowEvent;
   AliMCEvent*  mcEvent = MCEvent();                              // from TaskSE
   AliESDEvent* myESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
   AliAODEvent* myAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
+  AliMultiplicity* myTracklets = NULL;
+  AliESDPmdTrack* pmdtracks = NULL;//pmd      
+
+  int availableINslot=1;
+
+  if (!(fCutsRP&&fCutsPOI&&fCutsEvent))
+  {
+    AliError("cuts not set");
+    return;
+  }
+
+  //DEFAULT - automatically takes care of everything
+  if (fAnalysisType == "AUTOMATIC")
+  {
+    //check event cuts
+    if (!fCutsEvent->IsSelected(InputEvent())) return;
+
+    //first attach all possible information to the cuts
+    fCutsRP->SetEvent( InputEvent(), MCEvent() );  //attach event
+    fCutsPOI->SetEvent( InputEvent(), MCEvent() );
+
+    //then make the event
+    fFlowEvent->Fill( fCutsRP, fCutsPOI );
+    //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI );
+
+    if (myESD)
+      fFlowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent()));
+    if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+  }
 
   // Make the FlowEvent for MC input
-  if (fAnalysisType == "MC")
+  else if (fAnalysisType == "MC")
   {
     // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
     // This handler can return the current MC event
@@ -204,8 +338,9 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
       AliWarning("Event does not pass multiplicity cuts"); return;
     }
     //make event
-    flowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
+    fFlowEvent = new AliFlowEvent(mcEvent,fCFManager1,fCFManager2);
   }
+
   // Make the FlowEvent for ESD input
   else if (fAnalysisType == "ESD")
   {
@@ -221,7 +356,7 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
     }
 
     //check the offline trigger (check if the event has the correct trigger)
-    AliInfo(Form("There are %d tracks in this event", fInputEvent->GetNumberOfTracks()));
+    AliInfo(Form("ESD has %d tracks", fInputEvent->GetNumberOfTracks()));
 
     //check multiplicity
     if (!fCFManager1->CheckEventCuts(AliCFManager::kEvtRecCuts,myESD))
@@ -230,8 +365,42 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
     }
 
     //make event
-    flowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
+    if (fRPType == "Global") {
+      fFlowEvent = new AliFlowEvent(myESD,fCFManager1,fCFManager2);
+    }
+    else if (fRPType == "TPCOnly") {
+      fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kFALSE);
+    }
+    else if (fRPType == "TPCHybrid") {
+      fFlowEvent = new AliFlowEvent(myESD,fCFManager2,kTRUE);
+    }
+    else if (fRPType == "Tracklet"){
+      fFlowEvent = new AliFlowEvent(myESD,myTracklets,fCFManager2);
+    }
+    else if (fRPType == "FMD"){
+      TList* FMDdata = dynamic_cast<TList*>(GetInputData(availableINslot++));
+      if(!FMDdata) {
+        cout<<" No FMDdata "<<endl;
+        return;
+      }
+      TH2F* histFMD = dynamic_cast<TH2F*>(FMDdata->FindObject("dNdetadphiHistogramTrVtx"));
+      if (!histFMD) {
+        cout<< "No histFMD"<<endl;
+        return;
+      }
+      fFlowEvent = new AliFlowEvent(myESD,histFMD,fCFManager2);
+    }
+    else if (fRPType == "PMD"){
+      fFlowEvent = new AliFlowEvent(myESD,pmdtracks,fCFManager2);
+    }
+    else return;
+    
+    // if monte carlo event get reaction plane from monte carlo (depends on generator)
+    if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+    //set reference multiplicity, TODO: maybe move it to the constructor?
+    fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
   }
+
   // Make the FlowEvent for ESD input combined with MC info
   else if (fAnalysisType == "ESDMCkineESD" || fAnalysisType == "ESDMCkineMC" )
   {
@@ -265,12 +434,17 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
     //make event
     if (fAnalysisType == "ESDMCkineESD")
     {
-      flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
+      fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kESDkine, fCFManager1, fCFManager2 );
     }
     else if (fAnalysisType == "ESDMCkineMC")
     {
-      flowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
+      fFlowEvent = new AliFlowEvent(myESD, mcEvent, AliFlowEvent::kMCkine, fCFManager1, fCFManager2 );
     }
+    if (!fFlowEvent) return;
+    // if monte carlo event get reaction plane from monte carlo (depends on generator)
+    if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent);
+    //set reference multiplicity, TODO: maybe move it to the constructor?
+    fFlowEvent->SetReferenceMultiplicity(AliESDtrackCuts::GetReferenceMultiplicity(myESD,kTRUE));
   }
   // Make the FlowEventSimple for AOD input
   else if (fAnalysisType == "AOD")
@@ -280,57 +454,77 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *)
       AliError("ERROR: AOD not available");
       return;
     }
-    AliInfo(Form("There are %d tracks in this event", myAOD->GetNumberOfTracks()));
-    flowEvent = new AliFlowEvent(myAOD);
+    AliInfo(Form("AOD has %d tracks", myAOD->GetNumberOfTracks()));
+    fFlowEvent = new AliFlowEvent(myAOD);
   }
 
+  //inject candidates
+  if(fLoadCandidates) {
+    TObjArray* Candidates = dynamic_cast<TObjArray*>(GetInputData(availableINslot++));
+    //if(Candidates->GetEntriesFast()) 
+      //printf("I received %d candidates\n",Candidates->GetEntriesFast());
+    if (Candidates)
+    {
+      for(int iCand=0; iCand!=Candidates->GetEntriesFast(); ++iCand ) {
+        AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(Candidates->At(iCand));
+        if (!cand) continue;
+        cand->SetForPOISelection(kTRUE);
+        cand->SetForRPSelection(kFALSE);
+        //printf(" Ⱶ Checking at candidate %d with %d daughters\n",iCand,cand->GetNDaughters());
+        for(int iDau=0; iDau!=cand->GetNDaughters(); ++iDau) {
+          //printf("    Ⱶ Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau) );
+          for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs ) {
+            AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
+            if (!iRP) continue;
+            if( !iRP->InRPSelection() )
+              continue;
+            if( cand->GetIDDaughter(iDau) == iRP->GetID() ) {
+              //printf(" was in RP set");
+              cand->SetDaughter( iDau, iRP );
+              //temporarily untagging all daugters
+              iRP->SetForRPSelection(kFALSE);
+            }
+          }
+          //printf("\n");
+        }
+        fFlowEvent->AddTrack(cand);
+      }
+    }
+  }
+
+  if (!fFlowEvent) return; //shuts up coverity
+
   //check final event cuts
-  Int_t mult = flowEvent->NumberOfTracks();
-  if (mult<fMinMult && mult>fMaxMult)
+  Int_t mult = fFlowEvent->NumberOfTracks();
+  AliInfo(Form("FlowEvent has %i tracks",mult));
+  if (mult<fMinMult || mult>fMaxMult)
   {
     AliWarning("FlowEvent cut on multiplicity"); return;
   }
 
-  //tag subevents
-  flowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
-
-  ////TODO afterburner
-  //if (fbAfterburnerOn && fMyTRandom3) {
-  //  // set the new value of the values using a after burner
-  //  cout << "settings for afterburner in TaskFlowEvent.cxx:" << endl;
-  //  cout << "fCount = " << fCount << endl;
-  //  cout << "fNoOfLoops = " << fNoOfLoops << endl;
-  //  cout << "fEllipticFlowValue = " << fEllipticFlowValue << endl;
-  //  cout << "fSigmaEllipticFlowValue = " << fSigmaEllipticFlowValue << endl;
-  //  cout << "fMultiplicityOflowEvent = " << fMultiplicityOflowEvent << endl;
-  //  cout << "fSigmaMultiplicityOflowEvent = " << fSigmaMultiplicityOflowEvent << endl;
-  //  Double_t xRPangle;
-  //  if (!mcEvent) { xRPangle = TMath::TwoPi()*(fMyTRandom3->Rndm()); }
-  //  else { xRPangle = fRP; }
-  //  Double_t xNewFlowValue = fMyTRandom3->Gaus(fEllipticFlowValue,fSigmaEllipticFlowValue);
-  //  Int_t nNewMultOflowEvent =  TMath::Nint(fMyTRandom3->Gaus(fMultiplicityOflowEvent,fSigmaMultiplicityOflowEvent));
-  //  cout << "xRPangle = " << xRPangle << endl;
-  //  cout << "xNewFlowValue = " << xNewFlowValue << endl;
-  //  cout << "nNewMultOflowEvent = " << nNewMultOflowEvent << endl;
-  //  cout << "settings for after burner" << endl;
-  //  flowEventMaker->SetMCReactionPlaneAngle(xRPangle);
-  //  flowEventMaker->SetNoOfLoops(fNoOfLoops);
-  //  flowEventMaker->SetEllipticFlowValue(xNewFlowValue);
-  //  flowEventMaker->SetMultiplicityOflowEvent(nNewMultOflowEvent);
-  //  //end settings afterburner
-  //}
-
-  // if monte carlo event get reaction plane from monte carlo (depends on generator)
-  if (mcEvent && mcEvent->GenEventHeader()) flowEvent->SetMCReactionPlaneAngle(mcEvent);
+  //define dead zone
+  fFlowEvent->DefineDeadZone(fExcludedEtaMin, fExcludedEtaMax, fExcludedPhiMin, fExcludedPhiMax );
 
-  //fListHistos->Print();
-  //  fOutputFile->WriteObject(flowEvent,"myFlowEventSimple");
-  PostData(1,flowEvent);
-  if (fQA)
+
+  //////////////////////////////////////////////////////////////////////////////
+  ///////////////////////////AFTERBURNER
+  if (fAfterburnerOn)
   {
-    PostData(2,fQAInt);
-    PostData(3,fQADiff);
+    //if reaction plane not set from elsewhere randomize it before adding flow
+    if (!fFlowEvent->IsSetMCReactionPlaneAngle())
+      fFlowEvent->SetMCReactionPlaneAngle(gRandom->Uniform(0.0,TMath::TwoPi()));
+
+    fFlowEvent->AddFlow(fV1,fV2,fV3,fV4,fV5);     //add flow
+    fFlowEvent->CloneTracks(fNonFlowNumberOfTrackClones); //add nonflow by cloning tracks
   }
+  //////////////////////////////////////////////////////////////////////////////
+
+  //tag subEvents
+  fFlowEvent->TagSubeventsInEta(fMinA,fMaxA,fMinB,fMaxB);
+
+  //fListHistos->Print();
+  //fOutputFile->WriteObject(fFlowEvent,"myFlowEventSimple");
+  PostData(1,fFlowEvent);
 }
 
 //________________________________________________________________________
@@ -339,4 +533,3 @@ void AliAnalysisTaskFlowEvent::Terminate(Option_t *)
   // Called once at the end of the query -- do not call in case of CAF
 }
 
-