]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/FLOW/Tasks/AliAnalysisTaskJetFlow.cxx
from Redmer Bertens:
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskJetFlow.cxx
index 8c1d7e011ea696eeefb972324eeac03943bdae78..dfe5a7531052e680278315b4664d10dd57a019dd 100644 (file)
  **************************************************************************/
 
 /* 
- * AliAnaysisTaskJet
+ * AliAnaysisTaskJetFlow
  * author: Redmer Alexander Bertens
  * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl 
  *
- * jet correlation task - correlates jets to the vzero reaction plane using
- * the scalar product method 
+ * Interface task between EMCal jet framework and flow package
  *
- * this task expects POI's in a TClonesArray, e.g. from running it after 
- * $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation()
- * if subtracted jet pt is stored by using jet->PtSub() run in default mode
- * to run on unsubtracted jets - or any other class inheriting from AliVParticle - call SetDoVParticleAnalysis(kTRUE)
+ * This task expects POI's in a TClonesArray  (e.g. from running it after 
+ * $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation() )
+ * and connects them into an AliFlowEvent which is filled with either VZERO tracks or 
+ * TPC trakcs. 
+ * 
+ * POI's can be supplied as AliEmcalJets or as AliVParticles, note the different behavior
+ * with respect to the Pt value: for AliEmcalJets subtracted Pt is used by default, for VParticles
+ * Pt is taken directly.
  *
  * */
 
+// root includes
 #include <TChain.h>
 #include <TH1F.h>
+#include <TArrayD.h>
+#include <TProfile.h>
+#include <TString.h>
+#include <TList.h>
+// aliroot includes
+#include <AliAODEvent.h>
+#include <AliESDEvent.h>
+#include <AliEmcalJet.h>
 #include <AliAnalysisTask.h>
 #include <AliAnalysisManager.h>
 #include <AliFlowEventSimple.h>
+#include <AliFlowEvent.h>
 #include <AliFlowTrack.h>
 #include <AliFlowTrackCuts.h>
 #include <AliFlowEventCuts.h>
 #include <AliFlowCommonConstants.h>
-#include <TString.h>
-#include <AliAODEvent.h>
-#include <AliESDEvent.h>
-#include <AliEmcalJet.h>
-
+// local includes
 #include "AliAnalysisTaskJetFlow.h"
 
 class AliAnalysisTaskJetFlow;
@@ -51,24 +60,26 @@ using namespace std;
 ClassImp(AliAnalysisTaskJetFlow)
 
 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow() : AliAnalysisTaskSE(), 
-    fDebug(-1), fJetsName(0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(kFALSE), fInitialized(kFALSE), fPtBump(0), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPOIPtMin(0.15), fPOIPtMax(150), fCutsRP(0), fCutsPOI(0), fCutsNull(0), fCutsEvent(0), fFlowEvent(0), fHistAnalysisSummary(0), fCentralitySelection(0)
-{ /* default constructor */ }
+    fDebug(-1), fJetsName(0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(kFALSE), fDoTestFlowAnalysis(kFALSE), fInitialized(kFALSE), fPtBump(0), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPOIPtMin(0.15), fPOIPtMax(150), fPtBins(0), fCutsRP_TPC(0), fCutsRP_VZERO(0), fCutsPOI(0), fCutsNull(0), fCutsEvent(0), fFlowEvent_TPC(0), fFlowEvent_VZERO(0), fHistAnalysisSummary(0), fCentralitySelection(0), fv2VZEROA(0), fv2VZEROC(0)
+{ /* default constructor for ROOT IO */ }
 //_____________________________________________________________________________
 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow(const char* name) : AliAnalysisTaskSE(name),
-    fDebug(-1), fJetsName(0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(kFALSE), fInitialized(kFALSE), fPtBump(0), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPOIPtMin(0.15), fPOIPtMax(150), fCutsRP(0), fCutsPOI(0), fCutsNull(0), fCutsEvent(0), fFlowEvent(0), fHistAnalysisSummary(0), fCentralitySelection(0)
+    fDebug(-1), fJetsName(0), fOutputList(0), fDataType(kESD), fVParticleAnalysis(kFALSE), fDoTestFlowAnalysis(kFALSE), fInitialized(kFALSE), fPtBump(0), fCCMaxPt(150), fCCBinsInPt(50), fCentralityMin(-1), fCentralityMax(-1), fPOIPtMin(0.15), fPOIPtMax(150), fPtBins(0), fCutsRP_TPC(0), fCutsRP_VZERO(0), fCutsPOI(0), fCutsNull(0), fCutsEvent(0), fFlowEvent_TPC(0), fFlowEvent_VZERO(0), fHistAnalysisSummary(0), fCentralitySelection(0), fv2VZEROA(0), fv2VZEROC(0)
 {
     // constructor
     DefineInput(0, TChain::Class());
     DefineOutput(1, TList::Class());
     DefineOutput(2, AliFlowEventSimple::Class());
+    DefineOutput(3, AliFlowEventSimple::Class());
 }
 //_____________________________________________________________________________
 AliAnalysisTaskJetFlow::~AliAnalysisTaskJetFlow()
 {
     // destructor
-    if(fOutputList)     delete fOutputList;
-    if(fFlowEvent)      delete fFlowEvent;
-    if(fCutsEvent)      delete fCutsEvent;
+    if(fOutputList)             delete fOutputList;
+    if(fFlowEvent_TPC)          delete fFlowEvent_TPC;
+    if(fFlowEvent_VZERO)        delete fFlowEvent_VZERO;
+    if(fCutsEvent)              delete fCutsEvent;
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskJetFlow::LocalInit()
@@ -95,14 +106,36 @@ void AliAnalysisTaskJetFlow::UserCreateOutputObjects()
     fHistAnalysisSummary->SetBinContent(3, fCentralityMax);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "pt bias");
     fOutputList->Add(fHistAnalysisSummary);
+    if(fDoTestFlowAnalysis) { // set up the binning for the test flow analysis
+        if((!fPtBins) && fVParticleAnalysis) {
+            Double_t pt[51];
+            for(Int_t i(0); i < 51; i++) pt[i] = .1*i;
+            fPtBins = new TArrayD(51, pt);  // assume they will be charged particles
+        } else if (!fPtBins) {
+            Double_t pt[] = {0., 10., 20., 30., 40., 50., 80., 110., 140., 170., 200.};
+            fPtBins = new TArrayD(sizeof(pt)/sizeof(pt[0]), pt);     // assuming jets
+        }
+        Double_t bounds[fPtBins->GetSize()];
+        for(Int_t i(0); i < fPtBins->GetSize(); i++) bounds[i] = fPtBins->At(i);
+        fv2VZEROA = new TProfile("v2_EP_VZEROA", "v2_EP_VZEROA", fPtBins->GetSize()-1, bounds);
+        fv2VZEROC = new TProfile("v2_EP_VZEROC", "v2_EP_VZEROC", fPtBins->GetSize()-1, bounds);
+        fv2VZEROA->GetXaxis()->SetTitle("Pt [GeV/c]");
+        fv2VZEROA->GetYaxis()->SetTitle("v_{2}^{obs}");
+        fv2VZEROC->GetXaxis()->SetTitle("Pt [GeV/c]");
+        fv2VZEROC->GetYaxis()->SetTitle("v_{2}^{obs}");
+        fOutputList->Add(fv2VZEROA);
+        fOutputList->Add(fv2VZEROC);
+    }
     // qa
     fCentralitySelection = new TH1F("fCentralitySelection", "fCentralitySelection", 100, 0, 100);
     fOutputList->Add(fCentralitySelection);
-
     PostData(1, fOutputList);
     // create the flow event and configure the static cc object
-    fFlowEvent = new AliFlowEvent(1000);
-    PostData(2, fFlowEvent);
+
+    (fVParticleAnalysis) ? fFlowEvent_VZERO = new AliFlowEvent(10000) : fFlowEvent_VZERO = new AliFlowEvent(100);
+    PostData(2, fFlowEvent_VZERO);
+    fFlowEvent_TPC = new AliFlowEvent(10000);
+    PostData(3, fFlowEvent_TPC);
     AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
     cc->SetPtMax(fCCMaxPt+fPtBump);
     cc->SetNbinsPt(fCCBinsInPt);
@@ -112,7 +145,7 @@ void AliAnalysisTaskJetFlow::UserExec(Option_t *)
 {
     // user exec
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(!InputEvent() || !fCutsNull || !fCutsRP) return; // coverity (and sanity)
+    if(!InputEvent() || !fCutsNull || !fCutsRP_TPC || !fCutsRP_VZERO) return; // coverity (and sanity)
     if(!fInitialized) { 
         if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
         else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
@@ -129,11 +162,16 @@ void AliAnalysisTaskJetFlow::UserExec(Option_t *)
             return;
         }
         // prepare the track selection for RP's and the flow event
-        fCutsRP->SetEvent(InputEvent(), MCEvent());
+        fCutsRP_VZERO->SetEvent(InputEvent(), MCEvent());
+        fCutsRP_TPC->SetEvent(InputEvent(), MCEvent());
         fCutsNull->SetEvent(InputEvent(), MCEvent());
-        fFlowEvent->ClearFast();
-        fFlowEvent->Fill(fCutsRP, fCutsNull);
-        fFlowEvent->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
+        fFlowEvent_VZERO->ClearFast();
+        fFlowEvent_TPC->ClearFast();
+        // the event is filled with rp's only, poi's will be added manually
+        fFlowEvent_VZERO->Fill(fCutsRP_VZERO, fCutsNull);
+        fFlowEvent_TPC->Fill(fCutsRP_TPC, fCutsNull);
+        fFlowEvent_VZERO->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
+        fFlowEvent_TPC->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
         // loop over jets and inject them as POI's
         if(fVParticleAnalysis) {
             for(Int_t i(0); i < iJets; i++) {
@@ -148,14 +186,15 @@ void AliAnalysisTaskJetFlow::UserExec(Option_t *)
                     flowTrack->SetPt(jet->Pt() + fPtBump);
                     flowTrack->SetForPOISelection(kTRUE);
                     flowTrack->SetForRPSelection(kFALSE);
-                    fFlowEvent->InsertTrack(flowTrack);
+                    fFlowEvent_TPC->InsertTrack(flowTrack);
+                    fFlowEvent_VZERO->InsertTrack(flowTrack);
                 }
             }
         } else {
             for(Int_t i(0); i < iJets; i++) {
                 AliEmcalJet* jet = static_cast<AliEmcalJet*>(jets->At(i));
                 if(jet) {
-                    if(jet->PtSub() + fPtBump <= fPOIPtMin || jet->Pt() > fPOIPtMax) {
+                    if(jet->PtSub() + fPtBump <= fPOIPtMin || jet->PtSub() > fPOIPtMax) {
                         fHistAnalysisSummary->SetBinContent(4, 1);
                         continue;
                     }
@@ -164,7 +203,8 @@ void AliAnalysisTaskJetFlow::UserExec(Option_t *)
                     flowTrack->SetPt(jet->PtSub() + fPtBump);
                     flowTrack->SetForPOISelection(kTRUE);
                     flowTrack->SetForRPSelection(kFALSE);
-                    fFlowEvent->InsertTrack(flowTrack);
+                    fFlowEvent_TPC->InsertTrack(flowTrack);
+                    fFlowEvent_VZERO->InsertTrack(flowTrack);
                 }
             }
         }
@@ -174,14 +214,13 @@ void AliAnalysisTaskJetFlow::UserExec(Option_t *)
         if(fDebug > 0) printf(" > No accepted jets in event ! < " );
         return;
     }
+    fFlowEvent_TPC->TagSubeventsInEta(-10, 0, 0, 10);
+    fFlowEvent_VZERO->TagSubeventsInEta(-10, 0, 0, 10);
+    if(fDoTestFlowAnalysis) DoTestFlowAnalysis();
     fCentralitySelection->Fill(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
     PostData(1, fOutputList);
-    PostData(2, fFlowEvent);
-    if(fDebug>0) {
-        printf(" Event summary \n");
-        printf(" > number of POI's (jets) %i ", fFlowEvent->NumberOfTracks() - fFlowEvent->GetNumberOfRPs());
-        printf(" > number of RP's %i \n", fFlowEvent->GetNumberOfRPs());
-    }
+    PostData(2, fFlowEvent_VZERO);
+    PostData(3, fFlowEvent_TPC);
 }
 //_____________________________________________________________________________
 Bool_t AliAnalysisTaskJetFlow::PassesCuts(AliVEvent* event)
@@ -206,6 +245,45 @@ Bool_t AliAnalysisTaskJetFlow::PassesCuts(AliVEvent* event)
     return (cent <= fCentralityMin || cent > fCentralityMax || TMath::Abs(cent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) ? kFALSE : kTRUE;
 }
 //_____________________________________________________________________________
+void AliAnalysisTaskJetFlow::DoTestFlowAnalysis()
+{
+    // get a crude estimate of v2 based on the event plane method FIXME not tested !!!
+    if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
+    Double_t _a(0), _b(0), _c(0), _d(0);        // dummmy's
+    Double_t Q2a(InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, _a, _b));
+    Double_t Q2c(InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, _c, _d));
+    TProfile* a = (TProfile*)fv2VZEROA->Clone("temp_a");
+    TProfile* c = (TProfile*)fv2VZEROC->Clone("temp_c");
+    if(!(a||c)) return; // coverity
+    TClonesArray* jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName.Data()));
+    if(jets) {
+        Int_t iJets = jets->GetEntriesFast();
+        if(fVParticleAnalysis) {
+            for(Int_t i(0); i < iJets; i++) {
+                AliVParticle* jet = static_cast<AliVParticle*>(jets->At(i));
+                if(jet && jet->Pt() + fPtBump >= fPOIPtMin && jet->Pt() < fPOIPtMax) {
+                    a->Fill(jet->Pt(), TMath::Cos(2.*(jet->Phi()-Q2a)));
+                    c->Fill(jet->Pt(), TMath::Cos(2.*(jet->Phi()-Q2c)));
+                }
+            }
+        } else {
+            for(Int_t i(0); i < iJets; i++) {
+                AliEmcalJet* jet = static_cast<AliEmcalJet*>(jets->At(i));
+                if(jet && jet->PtSub() + fPtBump >= fPOIPtMin && jet->PtSub() < fPOIPtMax) {
+                    a->Fill(jet->PtSub(), TMath::Cos(2.*(jet->Phi()-Q2a)));
+                    c->Fill(jet->PtSub(), TMath::Cos(2.*(jet->Phi()-Q2c)));
+                }
+            }
+        }
+        for(Int_t i(0); i < fv2VZEROA->GetXaxis()->GetNbins(); i++) {
+            fv2VZEROA->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., a->GetBinContent(i+1));
+            fv2VZEROC->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., c->GetBinContent(i+1));
+        }
+        delete a;
+        delete c;
+    }
+}
+//_____________________________________________________________________________
 void AliAnalysisTaskJetFlow::Terminate(Option_t *)
 {
     // terminate