PWGJE/EMCALJetTasks/AliAnalysisTaskLocalRhoDev.h

   1 #ifndef AliAnalysisTaskLocalRhoDev_H
   2 #define AliAnalysisTaskLocalRhoDev_H
   3
   4 // $Id$
   5
   6 #include <AliAnalysisTaskEmcalJet.h>
   7 #include <AliEmcalJet.h>
   8 #include <AliVEvent.h>
   9 #include <AliVTrack.h>
  10 #include <AliVCluster.h>
  11 #include <TClonesArray.h>
  12 #include <TMath.h>
  13 #include <TRandom3.h>
  14 #include <AliLog.h>
  15 #include <AliJetContainer.h>
  16 #include <AliParticleContainer.h>
  17 #include <TF1.h>
  18 #include <TH1.h>
  19
  20 class THF1;
  21 class THF2;
  22 class TProfile;
  23 class AliLocalRhoParameter;
  24 class TArrayI;
  25
  26 class AliAnalysisTaskLocalRhoDev : public AliAnalysisTaskEmcalJet {
  27  public:
  28   // enumerators
  29   enum fitModulationType  { kNoFit, kV2, kV3, kCombined, kFourierSeries, kIntegratedFlow, kQC2, kQC4 }; // fit type
  30   enum fitGoodnessTest    { kChi2ROOT, kChi2Poisson, kLinearFit };
  31   enum detectorType       { kTPC, kVZEROA, kVZEROC, kVZEROComb};  // detector that was used
  32   enum qcRecovery         { kFixedRho, kNegativeVn, kTryFit };    // how to deal with negative cn value for qcn value
  33   enum runModeType        { kLocal, kGrid };                      // run mode type
  34   // constructors, destructor
  35   AliAnalysisTaskLocalRhoDev();
  36   AliAnalysisTaskLocalRhoDev(const char *name, runModeType type);
  37   virtual                 ~AliAnalysisTaskLocalRhoDev();
  38   // setting up the task and technical aspects
  39   void                    ExecOnce();
  40   Bool_t                  InitializeAnalysis();
  41   virtual void            UserCreateOutputObjects();
  42   TH1F*                   BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c = -1, Bool_t append = kTRUE);
  43   TH2F*                   BookTH2F(const char* name, const char* x, const char* y, Int_t binsx, Double_t minx, Double_t maxx,
  44                                    Int_t binsy, Double_t miny, Double_t maxy, Int_t c = -1, Bool_t append = kTRUE);
  45   virtual Bool_t          Run();
  46   /* inline */   Double_t PhaseShift(Double_t x) const {
  47     while (x>=TMath::TwoPi())x-=TMath::TwoPi();
  48     while (x<0.)x+=TMath::TwoPi();
  49     return x; }
  50   /* inline */   Double_t PhaseShift(Double_t x, Double_t n) const {
  51     x = PhaseShift(x);
  52     if(TMath::Nint(n)==2) while (x>TMath::Pi()) x-=TMath::Pi();
  53     if(TMath::Nint(n)==3) {
  54       if(x>2.*TMath::TwoPi()/n) x = TMath::TwoPi() - x;
  55       if(x>TMath::TwoPi()/n) x = TMath::TwoPi()-(x+TMath::TwoPi()/n);
  56     }
  57     return x; }
  58
  59   /* inline */    Double_t ChiSquarePDF(Int_t ndf, Double_t x) const {
  60     Double_t n(ndf/2.), denom(TMath::Power(2, n)*TMath::Gamma(n));
  61     if (denom!=0)  return ((1./denom)*TMath::Power(x, n-1)*TMath::Exp(-x/2.));
  62     return -999; }
  63   //note that the cdf of the chisquare distribution is the normalized lower incomplete gamma function
  64   /*inline */    Double_t ChiSquareCDF(Int_t ndf, Double_t x) const { return TMath::Gamma(ndf/2., x/2.); }
  65   /*inline */    Double_t ChiSquare(TH1& histo, TF1* func) const {
  66     // evaluate the chi2 using a poissonian error estimate on bins
  67     Double_t chi2(0.);
  68     for(Int_t i(0); i < histo.GetXaxis()->GetNbins(); i++) {
  69       if(histo.GetBinContent(i+1) <= 0.) continue;
  70       chi2 += TMath::Power((histo.GetBinContent(i+1)-func->Eval(histo.GetXaxis()->GetBinCenter(1+i))), 2)/histo.GetBinContent(i+1);
  71     }
  72     return chi2;
  73   }
  74   /*inline*/ Double_t KolmogorovTest(TH1F& histo, TF1* func, Bool_t toy = kTRUE) const {
  75     // return the probability from a Kolmogorov test
  76     TH1F test(histo);       // stack copy of test statistic
  77     for(Int_t i(0); i < test.GetXaxis()->GetNbins(); i++) test.SetBinContent(i+1, func->Eval(test.GetXaxis()->GetBinCenter(1+i)));
  78     if (toy) return histo.TH1::KolmogorovTest((&test), "X");
  79     else return histo.TH1::KolmogorovTest((&test));
  80   }
  81
  82   // setters - setup how to run
  83   void                    SetDebugMode(Int_t d)                           {fDebug = d;}
  84   void                    SetCentralityClasses(TArrayD* c)                {fCentralityClasses = c;}
  85   void                    SetAttachToEvent(Bool_t a)                      {fAttachToEvent = a;}
  86   void                    SetUseScaledRho(Bool_t s)                       {fUseScaledRho = s;}
  87   void                    SetFillHistograms(Bool_t b)                     {fFillHistograms = b;}
  88   // setters - analysis details
  89   void                    SetNoEventWeightsForQC(Bool_t e)                {fNoEventWeightsForQC = e;}
  90   void                    SetIntegratedFlow(TH1F* i, TH1F* j)             {fUserSuppliedV2 = i; fUserSuppliedV3 = j; }
  91   void                    SetOnTheFlyResCorrection(TH1F* r2, TH1F* r3)    {fUserSuppliedR2 = r2; fUserSuppliedR3 = r3; }
  92   void                    SetModulationFit(TF1* fit);
  93   void                    SetModulationFitMinMaxP(Float_t m, Float_t n)   {fMinPvalue = m; fMaxPvalue = n; }
  94   void                    SetExpectedRuns(TArrayI* r)                     {fExpectedRuns = r; }
  95   void                    SetExpectedSemiGoodRuns(TArrayI* r)             {fExpectedSemiGoodRuns = r;}
  96   void                    SetNameRhoSmall(TString s)                      {fNameSmallRho = s;}
  97   void                    SetModulationFitType(fitModulationType type)    {fFitModulationType = type; }
  98   void                    SetControlFunction(TF1* func)                   {fFitControl = func; }
  99   void                    SetFitGoodnessType(fitGoodnessTest test)        {fFitGoodnessTest = test; }
 100   void                    SetQCnRecoveryType(qcRecovery type)             {fQCRecovery = type; }
 101   void                    SetModulationFitOptions(TString opt)            {fFitModulationOptions = opt; }
 102   void                    SetReferenceDetector(detectorType type)         {fDetectorType = type; }
 103   void                    SetUsePtWeight(Bool_t w)                        {fUsePtWeight = w; }
 104   void                    SetUsePtWeightErrorPropagation(Bool_t w)        {fUsePtWeightErrorPropagation = w;}
 105   void                    SetRunModeType(runModeType type)                {fRunModeType = type; }
 106   void                    SetForceAbsVnHarmonics(Bool_t f)                {fAbsVnHarmonics = f; }
 107   void                    SetExcludeLeadingJetsFromFit(Float_t n)         {fExcludeLeadingJetsFromFit = n; }
 108   void                    SetRebinSwapHistoOnTheFly(Bool_t r)             {fRebinSwapHistoOnTheFly = r; }
 109   void                    SetSaveThisPercentageOfFits(Float_t p)          {fPercentageOfFits = p; }
 110   void                    SetUseV0EventPlaneFromHeader(Bool_t h)          {fUseV0EventPlaneFromHeader = h;}
 111   void                    SetSoftTrackMinMaxPt(Float_t min, Float_t max)  {fSoftTrackMinPt = min; fSoftTrackMaxPt = max;}
 112   // getters
 113   TString                 GetLocalRhoName() const                         {return fLocalRhoName; }
 114   // numerical evaluations
 115   void                    CalculateEventPlaneVZERO(Double_t vzero[2][2]) const;
 116   void                    CalculateEventPlaneTPC(Double_t* tpc);
 117   void                    CalculateEventPlaneCombinedVZERO(Double_t* comb) const;
 118   Double_t                CalculateQC2(Int_t harm);
 119   Double_t                CalculateQC4(Int_t harm);
 120   // helper calculations for the q-cumulant analysis, also used by AliAnalyisTaskJetFlow
 121   void                    QCnQnk(Int_t n, Int_t k, Double_t &reQ, Double_t &imQ);
 122   Double_t                QCnS(Int_t i, Int_t j);
 123   Double_t                QCnM();
 124   Double_t                QCnM11();
 125   Double_t                QCnM1111();
 126   Bool_t                  QCnRecovery(Double_t psi2, Double_t psi3);
 127   // analysis details
 128   Bool_t                  CorrectRho(Double_t psi2, Double_t psi3);
 129   void                    FillEventPlaneHistograms(Double_t psi2, Double_t psi3) const;
 130   void                    FillAnalysisSummaryHistogram() const;
 131   // track selection
 132   /* inline */     Bool_t PassesCuts(AliVTrack* track) const { return AcceptTrack(track, 0);}
 133   /* inline */     Bool_t PassesCuts(AliEmcalJet* jet) { return AcceptJet(jet, 0);}
 134   /* inline*/      Bool_t PassesSimpleCuts(AliEmcalJet* jet) {
 135   Float_t minPhi(GetJetContainer()->GetJetPhiMin()), maxPhi(GetJetContainer()->GetJetPhiMax());
 136   return (jet && jet->Pt() > 1 && jet->Eta() < .9-GetJetRadius() && jet->Eta() > -.9+GetJetRadius() && jet->Phi() > minPhi && jet->Phi() < maxPhi && jet->Area() > .557*GetJetRadius()*GetJetRadius()*TMath::Pi());
 137   }
 138
 139   /* inline */     AliEmcalJet* GetLeadingJet() {
 140       Int_t iJets(fJets->GetEntriesFast());
 141       Double_t pt(0);
 142       AliEmcalJet* leadingJet(0x0);
 143       for(Int_t i(0); i < iJets; i++) {
 144           AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 145           if(!PassesSimpleCuts(jet)) continue;
 146           if(jet->Pt() > pt) {
 147               leadingJet = jet;
 148               pt = leadingJet->Pt();
 149           }
 150       }
 151       return leadingJet;
 152   }
 153   Bool_t                  PassesCuts(AliVEvent* event);
 154   virtual void            Terminate(Option_t* option);
 155
 156  private:
 157   Int_t                   fDebug;                 // debug level (0 none, 1 fcn calls, 2 verbose)
 158   Bool_t                  fInitialized;           //! is the analysis initialized?
 159   Bool_t                  fAttachToEvent;         // attach local rho to the event
 160   Bool_t                  fFillHistograms;        // fill qa histograms
 161   Bool_t                  fNoEventWeightsForQC;   // don't store event weights for qc analysis
 162   Bool_t                  fUseScaledRho;          // use scaled rho
 163   TArrayD*                fCentralityClasses;     // centrality classes (maximum 10) used for QA
 164   TH1F*                   fUserSuppliedV2;        // histo with integrated v2
 165   TH1F*                   fUserSuppliedV3;        // histo with integrated v3
 166   TH1F*                   fUserSuppliedR2;        // correct the extracted v2 with this r
 167   TH1F*                   fUserSuppliedR3;        // correct the extracted v3 with this r
 168   Int_t                   fNAcceptedTracks;       //! number of accepted tracks
 169   Int_t                   fNAcceptedTracksQCn;    //! accepted tracks for QCn
 170   Int_t                   fInCentralitySelection; //! centrality bin, only for QA plots
 171   fitModulationType       fFitModulationType;     // fit modulation type
 172   fitGoodnessTest         fFitGoodnessTest;       // goodness of fit test
 173   qcRecovery              fQCRecovery;            // recovery type for e-by-e qc method
 174   Bool_t                  fUsePtWeight;           // use dptdphi instead of dndphi
 175   Bool_t                  fUsePtWeightErrorPropagation; // recalculate the bin error on the dpt dphi histogram
 176   detectorType            fDetectorType;          // type of detector used for modulation fit
 177   TString                 fFitModulationOptions;  // fit options for modulation fit
 178   runModeType             fRunModeType;           // run mode type
 179   TF1*                    fFitModulation;         // modulation fit for rho
 180   TF1*                    fFitControl;            // control function
 181   Float_t                 fMinPvalue;             // minimum value of p
 182   Float_t                 fMaxPvalue;             // maximum value of p
 183   TArrayI*                fExpectedRuns;          // list of known run numbers with default cuts
 184   TArrayI*                fExpectedSemiGoodRuns;  // list of runs that are marked as semi-good in the rct
 185   Int_t                   fRunNumber;             //! current run number
 186   AliRhoParameter*        fCachedRho;             //! cached rho object
 187   TString                 fNameSmallRho;          // name of rho object for semi-good tpc runs
 188   // additional jet cuts (most are inherited)
 189   Float_t                 fLocalJetMinEta;        // local eta cut for jets
 190   Float_t                 fLocalJetMaxEta;        // local eta cut for jets
 191   Float_t                 fLocalJetMinPhi;        // local phi cut for jets
 192   Float_t                 fLocalJetMaxPhi;        // local phi cut for jets
 193   Float_t                 fSoftTrackMinPt;        // min pt for soft tracks
 194   Float_t                 fSoftTrackMaxPt;        // max pt for soft tracks
 195   Double_t                fSemiGoodJetMinPhi;     // min phi for semi good tpc runs
 196   Double_t                fSemiGoodJetMaxPhi;     // max phi for semi good tpc runs
 197   Double_t                fSemiGoodTrackMinPhi;   // min phi for semi good tpc runs
 198   Double_t                fSemiGoodTrackMaxPhi;   // max phi for semi good tpc runs
 199
 200   // general qa histograms
 201   TH1F*                   fHistPvalueCDFROOT;     //! pdf value of chisquare p
 202   TH2F*                   fHistPvalueCDFROOTCent; //! p value versus centrlaity from root
 203   TH2F*                   fHistChi2ROOTCent;      //! reduced chi2 from ROOT, centrality correlation
 204   TH2F*                   fHistPChi2Root;         //! correlation p value and reduced chi2
 205   TH1F*                   fHistPvalueCDF;         //! cdf value of chisquare p
 206   TH2F*                   fHistPvalueCDFCent;     //! p value vs centrality
 207   TH2F*                   fHistChi2Cent;          //! reduced chi2, centrlaity correlation
 208   TH2F*                   fHistPChi2;             //! correlation p value and reduced chi2
 209   TH2F*                   fHistRhoStatusCent;     //! status of rho as function of centrality
 210
 211   // general settings
 212   Bool_t                  fAbsVnHarmonics;        // force postive local rho
 213   Float_t                 fExcludeLeadingJetsFromFit;    // exclude n leading jets from fit
 214   Bool_t                  fRebinSwapHistoOnTheFly;       // rebin swap histo on the fly
 215   Float_t                 fPercentageOfFits;      // save this percentage of fits
 216   Bool_t                  fUseV0EventPlaneFromHeader;    // use the vzero event plane from the header
 217   // transient object pointers
 218   TList*                  fOutputList;            //! output list
 219   TList*                  fOutputListGood;        //! output list for local analysis
 220   TList*                  fOutputListBad;         //! output list for local analysis
 221   TH1F*                   fHistSwap;              //! swap histogram
 222   TH1F*                   fHistAnalysisSummary;   //! flags
 223   TProfile*               fProfV2;                //! extracted v2
 224   TProfile*               fProfV2Cumulant;        //! v2 cumulant
 225   TProfile*               fProfV3;                //! extracted v3
 226   TProfile*               fProfV3Cumulant;        //! v3 cumulant
 227   TH1F*                   fHistPsi2[10];          //! psi 2
 228   TH1F*                   fHistPsi3[10];          //! psi 3
 229
 230   AliAnalysisTaskLocalRhoDev(const AliAnalysisTaskLocalRhoDev&);                  // not implemented
 231   AliAnalysisTaskLocalRhoDev& operator=(const AliAnalysisTaskLocalRhoDev&);       // not implemented
 232
 233   ClassDef(AliAnalysisTaskLocalRhoDev, 6);
 234 };
 235 #endif