PWGLF/SPECTRA/ChargedHadrons/dNdPt/AlidNdPtAnalysisPbPbAOD.h

   1 #ifndef AlidNdPtAnalysisPbPbAOD_H
   2 #define AlidNdPtAnalysisPbPbAOD_H
   3
   4
   5 //------------------------------------------------------------------------------
   6 // AlidNdPtAnalysisPbPbAOD class used for dNdPt analysis in PbPb collision
   7 // via AODs
   8 //
   9 // Author: P. Luettig, 15.05.2013
  10 // last modified: 10.06.2014
  11 //------------------------------------------------------------------------------
  12
  13
  14
  15 class iostream;
  16
  17 #include "AliAnalysisTaskSE.h"
  18 #include "TObject.h"
  19 #include "TList.h"
  20 #include "TFile.h"
  21 #include "TH1.h"
  22 #include "TH2.h"
  23 #include "TH3.h"
  24 #include "THnSparse.h"
  25 #include "THn.h"
  26 #include "TClonesArray.h"
  27 #include "TString.h"
  28 #include "TProfile.h"
  29 #include "TVector2.h"
  30
  31 #include "TParticlePDG.h"
  32 #include "TDatabasePDG.h"
  33
  34 #include "AliLog.h"
  35 #include "AliCentrality.h"
  36 #include "AliAODEvent.h"
  37 #include "AliVEvent.h"
  38
  39 #include "AliInputEventHandler.h"
  40 #include "AliAODInputHandler.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliGenHijingEventHeader.h"
  46 #include "AliGenPythiaEventHeader.h"
  47 #include "AliExternalTrackParam.h"
  48 #include "AliESDtrack.h"
  49 #include "AliEventplane.h"
  50
  51 #include "TSystem.h"
  52 #include "TROOT.h"
  53
  54 class AlidNdPtAnalysisPbPbAOD : public AliAnalysisTaskSE {
  55   public :
  56     enum CheckQuantity { cqCrossedRows = 0, cqNcluster = 1, cqChi = 2, cqLength = 3, cqRowsOverFindable = 4 };
  57     enum KinematicQuantity { kqPt = 0, kqEta = 1, kqPhi = 2 };
  58     enum MaxCheckQuantity { cqMax = 5 };
  59     enum MaxKinematicQuantity { kqMax = 3 };
  60
  61     AlidNdPtAnalysisPbPbAOD(const char *name = "dNdPtPbPbAOD");
  62     ~AlidNdPtAnalysisPbPbAOD();
  63
  64     virtual void UserCreateOutputObjects();
  65     virtual void UserExec(Option_t *option);
  66     virtual void Terminate(Option_t *);
  67
  68     // Set binning for Histograms (if not set default binning is used)
  69     void SetBinsMult(Int_t nbins, Double_t* edges)                      { Printf("[I] Setting Mult Bins"); fMultNbins = nbins; fBinsMult = GetArrayClone(nbins,edges); }
  70     void SetBinsPt(Int_t nbins, Double_t* edges)                        { Printf("[I] Setting pT Bins"); fPtNbins = nbins; fBinsPt = GetArrayClone(nbins,edges); }
  71     void SetBinsPtCorr(Int_t nbins, Double_t* edges)            { Printf("[I] Setting pTcorr Bins"); fPtCorrNbins = nbins; fBinsPtCorr = GetArrayClone(nbins,edges); }
  72     void SetBinsPtCheck(Int_t nbins, Double_t* edges)           { Printf("[I] Setting pTcheck Bins"); fPtCheckNbins = nbins; fBinsPtCheck = GetArrayClone(nbins,edges); }
  73     void SetBinsEta(Int_t nbins, Double_t* edges)                       { Printf("[I] Setting Eta Bins"); fEtaNbins = nbins; fBinsEta = GetArrayClone(nbins,edges); }
  74     void SetBinsEtaCheck(Int_t nbins, Double_t* edges)          { Printf("[I] Setting EtaCheck Bins"); fEtaCheckNbins = nbins; fBinsEtaCheck = GetArrayClone(nbins,edges); }
  75     void SetBinsZv(Int_t nbins, Double_t* edges)                        { Printf("[I] Setting Zv Bins"); fZvNbins = nbins; fBinsZv= GetArrayClone(nbins,edges); }
  76     void SetBinsCentrality(Int_t nbins, Double_t* edges)        { Printf("[I] Setting Cent Bins"); fCentralityNbins = nbins; fBinsCentrality = GetArrayClone(nbins,edges); }
  77     void SetBinsPhi(Int_t nbins, Double_t* edges)                       { Printf("[I] Setting Phi Bins"); fPhiNbins = nbins; fBinsPhi = GetArrayClone(nbins,edges); }
  78     void SetBinsRunNumber(Int_t nbins, Double_t* edges)         { Printf("[I] Setting RunNumber Bins"); fRunNumberNbins = nbins; fBinsRunNumber = GetArrayClone(nbins,edges); }
  79
  80     // set event cut variables
  81     void SetCutMaxZVertex( Double_t d)                                      { fCutMaxZVertex = d; }
  82     Double_t GetCutMaxZVertex()                                                 { return fCutMaxZVertex; }
  83
  84     // set track kinematic cut parameters
  85     void SetCutPtRange(Double_t ptmin, Double_t ptmax)          { fCutPtMin = ptmin; fCutPtMax = ptmax; }
  86     Double_t GetCutPtMin()                                                          { return fCutPtMin; }
  87     Double_t GetCutPtMax()                                                          { return fCutPtMax; }
  88
  89     void SetCutEtaRange(Double_t etamin, Double_t etamax)       { fCutEtaMin = etamin; fCutEtaMax = etamax; }
  90     Double_t GetCutEtaMin()                                                         { return fCutEtaMin; }
  91     Double_t GetCutEtaMax()                                                         { return fCutEtaMax; }
  92
  93     void EnableRelativeCuts()                                                           { Printf("[I] Relative Cuts enabled"); fUseRelativeCuts = kTRUE; }
  94     Bool_t AreRelativeCutsEnabled()                                                     { return fUseRelativeCuts; }
  95
  96     // setter and getter track quality cut parameters
  97     void SetFilterBit(Int_t b)                                                          { fFilterBit = b; };
  98     Int_t GetFilterBit()                                                                        { return fFilterBit; }
  99
 100     void SetCutRequireTPCRefit(Bool_t *b)                                       { fCutRequireTPCRefit = b; }
 101     Bool_t IsTPCRefitRequired()                                                         { return fCutRequireTPCRefit; }
 102
 103     void SetCutRequireITSRefit(Bool_t *b)                                       { fCutRequireITSRefit = b; }
 104     Bool_t IsITSRefitRequired()                                                         { return fCutRequireITSRefit; }
 105
 106     void SetCutMinNClustersTPC(Double_t d)                                      { fCutMinNumberOfClusters = d; }
 107     Double_t GetCutMinNClustersTPC()                                            { return fCutMinNumberOfClusters; }
 108
 109     void SetCutPercMinNClustersTPC(Double_t d)                          { Printf("[I] Take only %.2f%% tracks with most clusters", d*100.); fCutPercMinNumberOfClusters = d; }
 110     Double_t GetCutPercMinNClustersTPC()                                        { return fCutPercMinNumberOfClusters; }
 111
 112     void SetCutMinNCrossedRowsTPC(Double_t d)                           { fCutMinNumberOfCrossedRows = d; }
 113     Double_t GetCutMinNCrossedRowsTPC()                                         { return fCutMinNumberOfCrossedRows; }
 114
 115     void SetCutPercMinNCrossedRowsTPC(Double_t d)                       { Printf("[I] Take only %.2f%% tracks with most crossedRows", d*100.); fCutPercMinNumberOfCrossedRows = d; }
 116     Double_t GetCutPercMinNCrossedRowsTPC()                                     { return fCutPercMinNumberOfCrossedRows; }
 117
 118     void SetCutMinRatioCrossedRowsOverFindableClustersTPC(Double_t d)   { fCutMinRatioCrossedRowsOverFindableClustersTPC = d; }
 119     Double_t GetCutMinRatioCrossedRowsOverFindableClustersTPC()                 { return fCutMinRatioCrossedRowsOverFindableClustersTPC; }
 120
 121     void SetCutLengthInTPCPtDependent(Bool_t b)                         { fCutLengthInTPCPtDependent = b; }
 122     Bool_t DoCutLengthInTPCPtDependent()                                        { return fCutLengthInTPCPtDependent; }
 123
 124     void SetPrefactorLengthInTPCPtDependent(Double_t d)         { fPrefactorLengthInTPCPtDependent = d; }
 125     Double_t GetPrefactorLengthInTPCPtDependent()                       { return fPrefactorLengthInTPCPtDependent; }
 126
 127     void SetCutMaxChi2PerClusterTPC(Double_t d)                         { fCutMaxChi2PerClusterTPC = d; }
 128     void SetCutMaxFractionSharedTPCClusters(Double_t d)         { fCutMaxFractionSharedTPCClusters = d; }
 129     void SetCutMaxDCAToVertexZ(Double_t d)                                      { fCutMaxDCAToVertexZ = d; }
 130     void SetCutMaxDCAToVertexXY(Double_t d)                             { fCutMaxDCAToVertexXY = d; }
 131     void SetCutMaxChi2PerClusterITS(Double_t d)                         { fCutMaxChi2PerClusterITS = d; }
 132     void SetCutDCAToVertex2D(Bool_t *b)                                         { fCutDCAToVertex2D = b; }
 133     void SetCutRequireSigmaToVertex(Bool_t *b)                          { fCutRequireSigmaToVertex = b; }
 134     void SetCutMaxDCAToVertexXYPtDep(Double_t d0, Double_t d1, Double_t d2)
 135     {
 136       fCutMaxDCAToVertexXYPtDepPar0 = d0;
 137       fCutMaxDCAToVertexXYPtDepPar1 = d1;
 138       fCutMaxDCAToVertexXYPtDepPar2 = d2;
 139     }
 140     void SetCutAcceptKinkDaughters(Bool_t *b)                           { fCutAcceptKinkDaughters = b; }
 141     void SetCutMaxChi2TPCConstrainedGlobal(Double_t d)          { fCutMaxChi2TPCConstrainedGlobal = d; }
 142
 143     // fill function for cross check histos
 144     Bool_t FillDebugHisto(Double_t *dCrossCheckVar, Double_t *dKineVar, Double_t dCentrality, Bool_t bIsAccepted);
 145
 146     // fill function for cut settings
 147     void StoreCutSettingsToHistogram();
 148
 149     // getter for DCA
 150     Bool_t GetDCA(const AliAODTrack *track, AliAODEvent *evt, Double_t d0z0[2]);
 151
 152     THnSparseF * GetHistZvPtEtaCent() const { return fZvPtEtaCent; }
 153     TH1F * GetHistEventStatistics() const { return fEventStatistics; }
 154
 155     const char * GetParticleName(Int_t pdg);
 156
 157     AliGenHijingEventHeader* GetHijingEventHeader(AliAODMCHeader *header);
 158     AliGenPythiaEventHeader* GetPythiaEventHeader(AliAODMCHeader *header);
 159
 160         Double_t RotatePhi(Double_t phiTrack, Double_t phiEP);
 161 //      Double_t MoveEventplane(Double_t dMCEP);
 162
 163     Bool_t SetRelativeCuts(AliAODEvent *event);
 164
 165     Bool_t IsTrackAccepted(AliAODTrack *tr, Double_t dCentrality, Double_t bMagZ);
 166     Bool_t IsMCTrackAccepted(AliAODMCParticle *part);
 167
 168     Bool_t IsHijingParticle(const AliAODMCParticle *part, AliGenHijingEventHeader* hijingGenHeader);
 169     Bool_t IsPythiaParticle(const AliAODMCParticle *part, AliGenPythiaEventHeader* pythiaGenHeader);
 170
 171     static Double_t* GetArrayClone(Int_t n, Double_t* source);
 172
 173         void SetEventplaneSelector(char *c) { fEPselector = c; }
 174         TString GetEventplaneSelector() { return fEPselector; }
 175
 176   private :
 177
 178     // Output List
 179     TList               *fOutputList;
 180
 181     // Histograms
 182     TH1F            *fPt; // simple pT histogramm
 183     TH1F            *fMCPt; // simple pT truth histogramm
 184     THnSparseF  *fZvPtEtaCent; //-> Zv:Pt:Eta:Cent
 185     THnSparseF  *fDeltaphiPtEtaCent; //-> DeltaPhi:Pt:Eta:Cent
 186     THnSparseF  *fPtResptCent; //-> 1/pt:ResolutionPt:Cent
 187     THnSparseF  *fMCRecPrimZvPtEtaCent; //-> MC Zv:Pt:Eta:Cent
 188     THnSparseF  *fMCGenZvPtEtaCent; //-> MC Zv:Pt:Eta:Cent
 189     THnSparseF  *fMCRecSecZvPtEtaCent; //-> MC Zv:Pt:Eta:Cent, only secondaries
 190     THnSparseF  *fMCRecPrimDeltaphiPtEtaCent; //-> MC Phi:Pt:Eta:Cent
 191     THnSparseF  *fMCGenDeltaphiPtEtaCent; //-> MC Phi:Pt:Eta:Cent
 192     THnSparseF  *fMCRecSecDeltaphiPtEtaCent; //-> MC Phi:Pt:Eta:Cent, only secondaries
 193     TH1F            *fEventStatistics; // contains statistics of number of events after each cut
 194     TH1F        *fEventStatisticsCentrality; // contains number of events vs centrality, events need to have a track in kinematic range
 195     TH1F            *fMCEventStatisticsCentrality; // contains MC number of events vs centrality, events need to have a track in kinematic range
 196     TH1F            *fAllEventStatisticsCentrality; // contains number of events vs centrality, events need to be triggered
 197     TH2F            *fEventStatisticsCentralityTrigger; // contains number of events vs centrality in 1% bins vs trigger
 198     THnSparseF  *fZvMultCent; // Zv:Mult:Cent
 199     TH1F            *fTriggerStatistics; // contains number of events per trigger
 200     TH1F            *fCharge; // charge distribution in data
 201     TH1F            *fMCCharge; // charge distribution in MC
 202     THnSparseF  *fDCAPtAll; //control histo: DCAz:DCAxy:pT:eta:phi for all reconstructed tracks
 203     THnSparseF  *fDCAPtAccepted; //control histo: DCAz:DCAxy:pT:eta:phi for all accepted reco tracks
 204     THnSparseF  *fMCDCAPtSecondary; //control histo: DCAz:DCAxy:pT:eta:phi for all accepted reco track, which are secondaries (using MC info)
 205     THnSparseF  *fMCDCAPtPrimary; //control histo: DCAz:DCAxy:pT:eta:phi for all accepted reco track, which are primaries (using MC info)
 206     THnF            *fCrossCheckAll[5]; //control histo: {CrossedRows,Ncluster,Chi,Length,CrossedRows/Findable} vs pT,eta,phi,Centrality for all tracks
 207     THnF            *fCrossCheckAcc[5]; //control histo: {CrossedRows,Ncluster,Chi,Length,CrossedRows/Findable} vs pT,eta,phi,Centrality after cuts
 208     TH1F            *fCutPercClusters; // control histo: number of clusters, where the relative cut has been set e-by-e
 209     TH1F            *fCutPercCrossed; // control histo: number of crossed rows, where the relative cut has been set e-by-e
 210     TH2F            *fCrossCheckRowsLength; // control histo: number of crossed rows vs length in TPC
 211     TH2F            *fCrossCheckClusterLength; // control histo: number of clusters vs length in TPC
 212     TH2F            *fCrossCheckRowsLengthAcc; // control histo: number of crossed rows vs length in TPC for all accepted tracks
 213     TH2F            *fCrossCheckClusterLengthAcc; // control histo: number of clusters vs length in TPC for all accepted tracks
 214     TH2F                *fCrossCheckPtresLength; // control histo: relative pt resolution in 1/pt vs lenght in TPC
 215     TH2F                *fCrossCheckPtresRows; // control histo: relative pt resolution in 1/pt vs number of crossed rows in TPC
 216     TH1F        *fCutSettings; // control histo: cut settings
 217
 218     TH1F                *fEventplaneDist; // event plane distribution in phi
 219     TH2F                *fEventplaneRunDist; // event plane distribution in phi
 220     TH1F                *fMCEventplaneDist; // MC event plane distribution in phi
 221     TH2F                *fCorrelEventplaneMCDATA; // correlation between data and MC eventplane
 222     THnSparseF  *fCorrelEventplaneDefaultCorrected; // correlation between default and corrected (== subtraction of current track) eventplane
 223     TH2F                *fEventplaneSubtractedPercentage; // percentage of subtracted tracks
 224
 225         // cross check for event plane resolution
 226         TH2F            *fEPDistCent; // event plane distribution vs centrality
 227         TH2F            *fPhiCent;      // particle phi distribution vs centrality
 228         TProfile        *fPcosEPCent; // < cos 2 psi_ep > vs centrality
 229         TProfile        *fPsinEPCent; // < sin 2 psi_ep > vs centrality
 230         TProfile        *fPcosPhiCent; // < cos 2 phi > vs centrality
 231         TProfile        *fPsinPhiCent; // < sin 2 phi > vs centrality
 232
 233         // cross check for event plane determination
 234         TH2F            *fDeltaPhiCent; // DeltaPhi:Cent - DeltaPhi in the range from -pi to pi
 235
 236         THnSparseF      *fCrossCheckFilterBitPhiCent; // FilterBit:Phi:Centrality
 237
 238         // global variables
 239     Bool_t fIsMonteCarlo;
 240
 241         TString fEPselector;
 242
 243     // event cut variables
 244     Double_t fCutMaxZVertex;
 245
 246     // track kinematic cut variables
 247     Double_t fCutPtMin;
 248     Double_t fCutPtMax;
 249     Double_t fCutEtaMin;
 250     Double_t fCutEtaMax;
 251
 252     // track quality cut variables
 253     Int_t           fFilterBit;
 254     Bool_t          fUseRelativeCuts;
 255     Bool_t      fCutRequireTPCRefit;
 256     Bool_t          fCutRequireITSRefit;
 257     Double_t    fCutMinNumberOfClusters;
 258     Double_t    fCutPercMinNumberOfClusters;
 259     Double_t    fCutMinNumberOfCrossedRows;
 260     Double_t    fCutPercMinNumberOfCrossedRows;
 261     Double_t    fCutMinRatioCrossedRowsOverFindableClustersTPC;
 262     Double_t    fCutMaxChi2PerClusterTPC;
 263     Double_t    fCutMaxFractionSharedTPCClusters;
 264     Double_t    fCutMaxDCAToVertexZ;
 265     Double_t    fCutMaxDCAToVertexXY;
 266     Double_t    fCutMaxChi2PerClusterITS;
 267     Bool_t      fCutDCAToVertex2D;
 268     Bool_t          fCutRequireSigmaToVertex;
 269     Double_t    fCutMaxDCAToVertexXYPtDepPar0;
 270     Double_t    fCutMaxDCAToVertexXYPtDepPar1;
 271     Double_t    fCutMaxDCAToVertexXYPtDepPar2;
 272     Bool_t          fCutAcceptKinkDaughters;
 273     Double_t    fCutMaxChi2TPCConstrainedGlobal;
 274     Bool_t          fCutLengthInTPCPtDependent;
 275     Double_t    fPrefactorLengthInTPCPtDependent;
 276
 277     //binning for THNsparse
 278     Int_t       fMultNbins;
 279     Int_t       fPtNbins;
 280     Int_t       fPtCorrNbins;
 281     Int_t       fPtCheckNbins;
 282     Int_t       fEtaNbins;
 283     Int_t       fEtaCheckNbins;
 284     Int_t       fZvNbins;
 285     Int_t       fCentralityNbins;
 286     Int_t       fPhiNbins;
 287         Int_t           fRunNumberNbins;
 288     Double_t*   fBinsMult; //[fMultNbins]
 289     Double_t*   fBinsPt; //[fPtNbins]
 290     Double_t*   fBinsPtCorr; //[fPtCorrNbins]
 291     Double_t*   fBinsPtCheck; //[fPtCheckNbins]
 292     Double_t*   fBinsEta; //[fEtaNbins]
 293     Double_t*   fBinsEtaCheck; //[fEtaCheckNbins]
 294     Double_t*   fBinsZv; //[fZvNbins]
 295     Double_t*   fBinsCentrality; //[fCentralityNbins]
 296     Double_t*   fBinsPhi; //[fPhiNbins]
 297     Double_t*   fBinsRunNumber; //[fRunNumberNbins]
 298
 299     AlidNdPtAnalysisPbPbAOD(const AlidNdPtAnalysisPbPbAOD&); // not implemented
 300     AlidNdPtAnalysisPbPbAOD& operator=(const AlidNdPtAnalysisPbPbAOD&); // not implemented
 301
 302     ClassDef(AlidNdPtAnalysisPbPbAOD,13); // has to be at least 1, otherwise not streamable...
 303 };
 304
 305 #endif