PWG/FLOW/Tasks/AliAnalysisTaskJetFlow.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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  * AliAnaysisTaskJetFlow
  18  * author: Redmer Alexander Bertens
  19  * rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
  20  *
  21  * Interface task between EMCal jet framework and flow package
  22  *
  23  * This task expects POI's in a TClonesArray  (e.g. from running it after
  24  * $ALICE_ROOT/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation() )
  25  * and connects them into an AliFlowEvent which is filled with either VZERO tracks or
  26  * TPC trakcs.
  27  *
  28  * POI's can be supplied as AliEmcalJets or as AliVParticles, note the different behavior
  29  * with respect to the Pt value: for AliEmcalJets subtracted Pt is used by default, for VParticles
  30  * Pt is taken directly.
  31  *
  32  * */
  33
  34 // root includes
  35 #include <TChain.h>
  36 #include <TH1F.h>
  37 #include <TArrayD.h>
  38 #include <TProfile.h>
  39 #include <TString.h>
  40 #include <TList.h>
  41 // aliroot includes
  42 #include <AliAODEvent.h>
  43 #include <AliESDEvent.h>
  44 #include <AliEmcalJet.h>
  45 #include <AliAnalysisTask.h>
  46 #include <AliAnalysisManager.h>
  47 #include <AliFlowEventSimple.h>
  48 #include <AliFlowEvent.h>
  49 #include <AliFlowTrack.h>
  50 #include <AliFlowTrackCuts.h>
  51 #include <AliFlowEventCuts.h>
  52 #include <AliFlowCommonConstants.h>
  53 // local includes
  54 #include "AliAnalysisTaskJetFlow.h"
  55
  56 class AliAnalysisTaskJetFlow;
  57
  58 using namespace std;
  59
  60 ClassImp(AliAnalysisTaskJetFlow)
  61
  62 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow() : AliAnalysisTaskSE(),
  63     fDebug(-1), fExplicitOutlierCut(-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), fTempA(0), fTempC(0)
  64 { /* default constructor for ROOT IO */ }
  65 //_____________________________________________________________________________
  66 AliAnalysisTaskJetFlow::AliAnalysisTaskJetFlow(const char* name) : AliAnalysisTaskSE(name),
  67     fDebug(-1), fExplicitOutlierCut(-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), fTempA(0), fTempC(0)
  68 {
  69     // constructor
  70     DefineInput(0, TChain::Class());
  71     DefineOutput(1, TList::Class());
  72     DefineOutput(2, AliFlowEventSimple::Class());
  73     DefineOutput(3, AliFlowEventSimple::Class());
  74 }
  75 //_____________________________________________________________________________
  76 AliAnalysisTaskJetFlow::~AliAnalysisTaskJetFlow()
  77 {
  78     // destructor
  79     if(fOutputList)             delete fOutputList;
  80     if(fFlowEvent_TPC)          delete fFlowEvent_TPC;
  81     if(fFlowEvent_VZERO)        delete fFlowEvent_VZERO;
  82     if(fCutsEvent)              delete fCutsEvent;
  83     if(fCutsRP_VZERO)           delete fCutsRP_VZERO;
  84     if(fCutsRP_TPC)             delete fCutsRP_TPC;
  85     if(fCutsNull)               delete fCutsNull;
  86     if(fCutsPOI)                delete fCutsPOI;
  87     if(fPtBins)                 delete fPtBins;
  88     if(fTempA)                  delete fTempA;
  89     if(fTempC)                  delete fTempC;
  90 }
  91 //_____________________________________________________________________________
  92 void AliAnalysisTaskJetFlow::LocalInit()
  93 {
  94     // executed once
  95     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
  96     fCutsEvent = new AliFlowEventCuts();
  97     fCutsEvent->SetRefMultMethod(AliESDtrackCuts::kTrackletsITSTPC);
  98 }
  99 //_____________________________________________________________________________
 100 void AliAnalysisTaskJetFlow::UserCreateOutputObjects()
 101 {
 102     // create output objects
 103     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 104     fOutputList = new TList();
 105     fOutputList->SetOwner(kTRUE);
 106     // histograms
 107     fHistAnalysisSummary = new TH1F("fHistAnalysisSummary", "fHistAnalysisSummary", 4, -.5, 3.5);
 108     fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fDataType");
 109     fHistAnalysisSummary->SetBinContent(1, (int)fDataType);
 110     fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fCentralityMin");
 111     fHistAnalysisSummary->SetBinContent(2, fCentralityMin);
 112     fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fCentralityMax");
 113     fHistAnalysisSummary->SetBinContent(3, fCentralityMax);
 114     fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "pt bias");
 115     fOutputList->Add(fHistAnalysisSummary);
 116     if(fDoTestFlowAnalysis) { // set up the binning for the test flow analysis
 117         if((!fPtBins) && fVParticleAnalysis) {
 118             Double_t pt[51];
 119             for(Int_t i(0); i < 51; i++) pt[i] = .1*i;
 120             fPtBins = new TArrayD(51, pt);  // assume they will be charged particles
 121         } else if (!fPtBins) {
 122             Double_t pt[] = {0., 10., 20., 30., 40., 50., 80., 110., 140., 170., 200.};
 123             fPtBins = new TArrayD(sizeof(pt)/sizeof(pt[0]), pt);     // assuming jets
 124         }
 125         Double_t bounds[fPtBins->GetSize()];
 126         for(Int_t i(0); i < fPtBins->GetSize(); i++) bounds[i] = fPtBins->At(i);
 127         fv2VZEROA = new TProfile("v2_EP_VZEROA", "v2_EP_VZEROA", fPtBins->GetSize()-1, bounds);
 128         fv2VZEROC = new TProfile("v2_EP_VZEROC", "v2_EP_VZEROC", fPtBins->GetSize()-1, bounds);
 129         fv2VZEROA->GetXaxis()->SetTitle("Pt [GeV/c]");
 130         fv2VZEROA->GetYaxis()->SetTitle("v_{2}^{obs}");
 131         fv2VZEROC->GetXaxis()->SetTitle("Pt [GeV/c]");
 132         fv2VZEROC->GetYaxis()->SetTitle("v_{2}^{obs}");
 133         fOutputList->Add(fv2VZEROA);
 134         fOutputList->Add(fv2VZEROC);
 135         // bookkeeping histo's, will not be stored
 136         fTempA = (TProfile*)fv2VZEROA->Clone("temp_a");
 137         fTempC = (TProfile*)fv2VZEROC->Clone("temp_c");
 138     }
 139     // qa
 140     fCentralitySelection = new TH1F("fCentralitySelection", "fCentralitySelection", 100, 0, 100);
 141     fOutputList->Add(fCentralitySelection);
 142     PostData(1, fOutputList);
 143     // create the flow event and configure the static cc object
 144
 145     (fVParticleAnalysis) ? fFlowEvent_VZERO = new AliFlowEvent(10000) : fFlowEvent_VZERO = new AliFlowEvent(100);
 146     PostData(2, fFlowEvent_VZERO);
 147     fFlowEvent_TPC = new AliFlowEvent(10000);
 148     PostData(3, fFlowEvent_TPC);
 149     AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
 150     cc->SetPtMax(fCCMaxPt+fPtBump);
 151     cc->SetNbinsPt(fCCBinsInPt);
 152 }
 153 //_____________________________________________________________________________
 154 void AliAnalysisTaskJetFlow::UserExec(Option_t *)
 155 {
 156     // user exec
 157     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 158     if(!InputEvent() || !fCutsNull || !fCutsRP_TPC || !fCutsRP_VZERO) return; // coverity (and sanity)
 159     if(!fInitialized) {
 160         if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
 161         else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
 162         fInitialized = kTRUE;
 163     }
 164     if(!PassesCuts(InputEvent())) return;               // check the event cuts
 165     // get the jet array, which is added as an extension to the AliVEvent by the jetfinder
 166     TClonesArray* jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName.Data()));
 167     Int_t nAcceptedJets(0);
 168     if(jets) {
 169         Int_t iJets = jets->GetEntriesFast();
 170         if(iJets <= 0) {
 171             if(fDebug>0) printf(" > Retrieved empty AliVEvent extension, aborting ! < \n ");
 172             return;
 173         }
 174         // prepare the track selection for RP's and the flow event
 175         fCutsRP_VZERO->SetEvent(InputEvent(), MCEvent());
 176         fCutsRP_TPC->SetEvent(InputEvent(), MCEvent());
 177         fCutsNull->SetEvent(InputEvent(), MCEvent());
 178         fFlowEvent_VZERO->ClearFast();
 179         fFlowEvent_TPC->ClearFast();
 180         // the event is filled with rp's only, poi's will be added manually
 181         fFlowEvent_VZERO->Fill(fCutsRP_VZERO, fCutsNull);
 182         fFlowEvent_TPC->Fill(fCutsRP_TPC, fCutsNull);
 183         fFlowEvent_VZERO->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
 184         fFlowEvent_TPC->SetReferenceMultiplicity(fCutsEvent->RefMult(InputEvent(), MCEvent()));
 185         // loop over jets and inject them as POI's
 186         if(fVParticleAnalysis) {
 187             for(Int_t i(0); i < iJets; i++) {
 188                 AliVParticle* jet = static_cast<AliVParticle*>(jets->At(i));
 189                 if(jet) {
 190                     if(jet->Pt() + fPtBump <= fPOIPtMin || jet->Pt() > fPOIPtMax) {
 191                         fHistAnalysisSummary->SetBinContent(4, 1);
 192                         continue;
 193                     }
 194                     nAcceptedJets++;
 195                     // AliFlowTracks are created on the stack
 196                     AliFlowTrack flowTrack = AliFlowTrack(jet);
 197                     flowTrack.SetPt(jet->Pt() + fPtBump);
 198                     flowTrack.SetForPOISelection(kTRUE);
 199                     flowTrack.SetForRPSelection(kFALSE);
 200                     fFlowEvent_TPC->InsertTrack(&flowTrack);
 201                     fFlowEvent_VZERO->InsertTrack(&flowTrack);
 202                 }
 203             }
 204         } else {
 205             for(Int_t i(0); i < iJets; i++) {
 206                 AliEmcalJet* jet = static_cast<AliEmcalJet*>(jets->At(i));
 207                 if(jet) {
 208                     if(jet->PtSub() + fPtBump <= fPOIPtMin || jet->PtSub() > fPOIPtMax) {
 209                         fHistAnalysisSummary->SetBinContent(4, 1);
 210                         continue;
 211                     }
 212                     nAcceptedJets++;
 213                     AliFlowTrack flowTrack = AliFlowTrack(jet);
 214                     flowTrack.SetPt(jet->PtSub() + fPtBump);
 215                     flowTrack.SetForPOISelection(kTRUE);
 216                     flowTrack.SetForRPSelection(kFALSE);
 217                     fFlowEvent_TPC->InsertTrack(&flowTrack);
 218                     fFlowEvent_VZERO->InsertTrack(&flowTrack);
 219                 }
 220             }
 221         }
 222     }
 223     else if(fDebug > 0 ) printf(" Failed to find TClones Jet array while using name %s \n ", fJetsName.Data());
 224     if(nAcceptedJets < 1) {
 225         if(fDebug > 0) printf(" > No accepted jets in event ! < " );
 226         return;
 227     }
 228     fFlowEvent_TPC->TagSubeventsInEta(-10, 0, 0, 10);
 229     fFlowEvent_VZERO->TagSubeventsInEta(-10, 0, 0, 10);
 230     if(fDoTestFlowAnalysis) DoTestFlowAnalysis();
 231     fCentralitySelection->Fill(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 232     PostData(1, fOutputList);
 233     PostData(2, fFlowEvent_VZERO);
 234     PostData(3, fFlowEvent_TPC);
 235 }
 236 //_____________________________________________________________________________
 237 Bool_t AliAnalysisTaskJetFlow::PassesCuts(AliVEvent* event)
 238 {
 239     // event cuts
 240     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 241     if(!event) return kFALSE;
 242     if(TMath::Abs(InputEvent()->GetPrimaryVertex()->GetZ()) > 10.) return kFALSE;
 243     // aod and esd specific checks
 244     switch (fDataType) {
 245        case kESD: {
 246             AliESDEvent* esdEvent = static_cast<AliESDEvent*>(InputEvent());
 247             if( (!esdEvent) || (TMath::Abs(esdEvent->GetPrimaryVertexSPD()->GetZ() - esdEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
 248        } break;
 249        case kAOD: {
 250             AliAODEvent* aodEvent = static_cast<AliAODEvent*>(InputEvent());
 251             if( (!aodEvent) || (TMath::Abs(aodEvent->GetPrimaryVertexSPD()->GetZ() - aodEvent->GetPrimaryVertex()->GetZ()) > .5) ) return kFALSE;
 252        } break;
 253        default: break;
 254     }
 255     Float_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 256     if(fExplicitOutlierCut == 2010 || fExplicitOutlierCut == 2011) {
 257        if(!PassesCuts(fExplicitOutlierCut)) return kFALSE;
 258     }
 259     return (cent <= fCentralityMin || cent > fCentralityMax || TMath::Abs(cent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) ? kFALSE : kTRUE;
 260 }
 261 //_____________________________________________________________________________
 262 Bool_t AliAnalysisTaskJetFlow::PassesCuts(Int_t year)
 263 {
 264     // additional centrality cut based on relation between tpc and global multiplicity
 265     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 266     AliAODEvent* event(dynamic_cast<AliAODEvent*>(InputEvent()));
 267     if(!event) return kFALSE;
 268     Int_t multTPC(0), multGlob(0), nTracks(InputEvent()->GetNumberOfTracks());
 269     for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
 270         AliAODTrack* track = event->GetTrack(iTracks);
 271         if(!track) continue;
 272         if (!track || track->Pt() < .2 || track->Pt() > 5.0 || TMath::Abs(track->Eta()) > .8 || track->GetTPCNcls() < 70 || !track->GetDetPid() || track->GetDetPid()->GetTPCsignal() < 10.0)  continue;  // general quality cut
 273         if (track->TestFilterBit(1) && track->Chi2perNDF() > 0.2) multTPC++;
 274         if (!track->TestFilterBit(16) || track->Chi2perNDF() < 0.1) continue;
 275         Double_t b[2] = {-99., -99.};
 276         Double_t bCov[3] = {-99., -99., -99.};
 277         if (track->PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov) && TMath::Abs(b[0]) < 0.3 && TMath::Abs(b[1]) < 0.3) multGlob++;
 278     }
 279     if(year == 2010 && multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob)) return kTRUE;
 280     if(year == 2011  && multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob)) return kTRUE;
 281     return kFALSE;
 282 }
 283 //_____________________________________________________________________________
 284 void AliAnalysisTaskJetFlow::DoTestFlowAnalysis()
 285 {
 286     // get a crude estimate of v2 based on the event plane method FIXME not tested !!!
 287     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 288     Double_t _a(0), _b(0), _c(0), _d(0);        // dummmy's
 289     Double_t Q2a(InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, _a, _b));
 290     Double_t Q2c(InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, _c, _d));
 291     fTempA->Reset();    // clear the containers for a new iteration
 292     fTempC->Reset();
 293     TClonesArray* jets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName.Data()));
 294     if(jets) {
 295         Int_t iJets = jets->GetEntriesFast();
 296         if(fVParticleAnalysis) {
 297             for(Int_t i(0); i < iJets; i++) {
 298                 AliVParticle* jet = static_cast<AliVParticle*>(jets->At(i));
 299                 if(jet && jet->Pt() + fPtBump >= fPOIPtMin && jet->Pt() < fPOIPtMax) {
 300                     fTempA->Fill(jet->Pt(), TMath::Cos(PhaseShift(2.*(jet->Phi()-Q2a), 2)));
 301                     fTempC->Fill(jet->Pt(), TMath::Cos(PhaseShift(2.*(jet->Phi()-Q2c), 2)));
 302                 }
 303             }
 304         } else {
 305             for(Int_t i(0); i < iJets; i++) {
 306                 AliEmcalJet* jet = static_cast<AliEmcalJet*>(jets->At(i));
 307                 if(jet && jet->PtSub() + fPtBump >= fPOIPtMin && jet->PtSub() < fPOIPtMax) {
 308                     fTempA->Fill(jet->PtSub(), TMath::Cos(PhaseShift(2.*(jet->Phi()-Q2a), 2)));
 309                     fTempC->Fill(jet->PtSub(), TMath::Cos(PhaseShift(2.*(jet->Phi()-Q2c), 2)));
 310                 }
 311             }
 312         }
 313         for(Int_t i(0); i < fv2VZEROA->GetXaxis()->GetNbins(); i++) {
 314             fv2VZEROA->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., fTempA->GetBinContent(i+1));
 315             fv2VZEROC->Fill(fPtBins->At(i)+(fPtBins->At(i)+fPtBins->At(1+i))/2., fTempC->GetBinContent(i+1));
 316         }
 317     }
 318 }
 319 //_____________________________________________________________________________
 320 void AliAnalysisTaskJetFlow::Terminate(Option_t *)
 321 {
 322     // terminate
 323     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 324 }
 325 //_____________________________________________________________________________