PWG/FLOW/Tasks/AliJetFlowTools.h

   1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   2  * See cxx source for full Copyright notice */
   3  /* $Id$ */
   4
   5 /* Author: Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
   6  * see implementation for additional information */
   7
   8 #ifndef AliJetFlowTools_H
   9 #define AliJetFlowTools_H
  10
  11 // root forward declarations
  12 class TF1;
  13 class TH1D;
  14 class TH2D;
  15 class TCanvas;
  16 class TLegend;
  17 class TString;
  18 class TArrayD;
  19 class TGraph;
  20 class TGraphErrors;
  21 class TObjArray;
  22 // aliroot forward declarations
  23 class AliAnaChargedJetResponseMaker;
  24 class AliUnfolding;
  25 // root includes
  26 #include "TMatrixD.h"
  27 #include "TList.h"
  28 #include "TDirectoryFile.h"
  29 #include "TFile.h"
  30 #include "TProfile.h"
  31 #include "TVirtualPad.h"
  32 //_____________________________________________________________________________
  33 class AliJetFlowTools {
  34     public:
  35         AliJetFlowTools();
  36     protected:
  37         ~AliJetFlowTools();     // not implemented (deliberately). object ownership is a bit messy in this class
  38                                 // since most (or all) of the objects are owned by the input and output files
  39     public:
  40         // enumerators
  41         enum unfoldingAlgorithm {       // type of unfolding alrogithm
  42             kChi2,                      // chi^2 unfolding, implemented in AliUnfolding
  43             kBayesian,                  // Bayesian unfolding, implemented in RooUnfold
  44             kBayesianAli,               // Bayesian unfolding, implemented in AliUnfolding
  45             kSVD,                       // SVD unfolding, implemented in RooUnfold
  46             kNone };                    // no unfolding
  47         enum prior {                    // prior that is used for unfolding
  48             kPriorChi2,                 // prior from chi^2 method
  49             kPriorMeasured };           // use measured spectrum as prior
  50         enum histoType {                // histogram identifier, only used internally
  51             kInPlaneSpectrum,
  52             kOutPlaneSpectrum,
  53             kUnfoldedSpectrum,
  54             kFoldedSpectrum,
  55             kMeasuredSpectrum,
  56             kEmpty };
  57         // setters, interface to the class
  58         void            SetSaveFull(Bool_t b)           {fSaveFull              = b;}
  59         void            SetInputList(TList* list)       {
  60             fInputList          = list;
  61             fRefreshInput       = kTRUE;
  62         }
  63         void            SetOutputFileName(TString name) {fOutputFileName        = name;}
  64         void            CreateOutputList(TString name) {
  65             // create a new output list and add it to the full output
  66             if(!fOutputFile) fOutputFile = new TFile(fOutputFileName.Data(), "RECREATE");
  67             fOutputFile->cd();  // avoid nested dirs
  68             if(name.EqualTo(fActiveString)) {
  69                 printf(Form(" > Warning, duplicate output list, renaming %s to %s_2 ! < \n", name.Data(), name.Data()));
  70                 name+="_2";
  71             }
  72             fActiveString = name;
  73             fActiveDir = new TDirectoryFile(fActiveString.Data(), fActiveString.Data());
  74             fActiveDir->cd();
  75         }
  76         void            SetCentralityBin(Int_t bin)             {fCentralityBin         = bin;}
  77         void            SetDetectorResponse(TH2D* dr)           {fDetectorResponse      = dr;}
  78         void            SetJetFindingEfficiency(TH1D* e)        {fJetFindingEff         = e;}
  79         void            SetBinsTrue(TArrayD* bins)              {fBinsTrue              = bins;}
  80         void            SetBinsRec(TArrayD* bins)               {fBinsRec               = bins;}
  81         void            SetBinsTruePrior(TArrayD* bins)         {fBinsTruePrior         = bins;}
  82         void            SetBinsRecPrior(TArrayD* bins)          {fBinsRecPrior          = bins;}
  83         void            SetSVDReg(Int_t r)                      {fSVDRegIn              = r; fSVDRegOut         = r;}
  84         void            SetSVDReg(Int_t in, Int_t out)          {fSVDRegIn              = in; fSVDRegOut        = out;}
  85         void            SetSVDToy(Bool_t b, Float_t r)          {fSVDToy                = b; fJetRadius         = r;}
  86         void            SetBeta(Double_t b)                     {fBetaIn                = b; fBetaOut           = b;}
  87         void            SetBeta(Double_t i, Double_t o)         {fBetaIn                = i; fBetaOut           = o;}
  88         void            SetBayesianIter(Int_t i)                {fBayesianIterIn        = i; fBayesianIterOut    = i;}
  89         void            SetBayesianIter(Int_t i, Int_t o)       {fBayesianIterIn        = i; fBayesianIterOut    = o;}
  90         void            SetBayesianSmooth(Float_t s)            {fBayesianSmoothIn      = s; fBayesianSmoothOut = s;}
  91         void            SetBayesianSmooth(Float_t i, Float_t o) {fBayesianSmoothIn      = i; fBayesianSmoothOut = o;}
  92         void            SetAvoidRoundingError(Bool_t r)         {fAvoidRoundingError    = r;}
  93         void            SetUnfoldingAlgorithm(unfoldingAlgorithm ua)    {fUnfoldingAlgorithm                    = ua;}
  94         void            SetPrior(prior p)                       {fPrior                 = p;}
  95         void            SetNormalizeSpectra(Bool_t b)           {fNormalizeSpectra      = b;}
  96         void            SetNormalizeSpectra(Int_t e)            { // use to normalize to this no of events
  97             fEventCount         = e;
  98             fNormalizeSpectra   = kFALSE;
  99         }
 100         void            SetSmoothenPrior(Bool_t b, Float_t min = 50., Float_t max = 100., Float_t start= 75., Bool_t counts = kTRUE) {
 101             fSmoothenPrior      = b;
 102             fFitMin             = min;
 103             fFitMax             = max;
 104             fFitStart           = start;
 105             fSmoothenCounts     = counts;
 106         }
 107         void            SetTestMode(Bool_t t)                   {fTestMode              = t;}
 108         void            SetEventPlaneResolution(Double_t r)     {fEventPlaneRes         = r;}
 109         void            SetUseDetectorResponse(Bool_t r)        {fUseDetectorResponse   = r;}
 110         void            SetUseDptResponse(Bool_t r)             {fUseDptResponse        = r;}
 111         void            SetTrainPowerFit(Bool_t t)              {fTrainPower            = t;}
 112         void            SetDphiUnfolding(Bool_t i)              {fDphiUnfolding         = i;}
 113         void            SetDphiDptUnfolding(Bool_t i)           {fDphiDptUnfolding      = i;}
 114         void            SetExLJDpt(Bool_t i)                    {fExLJDpt               = i;}
 115         void            SetWeightFunction(TF1* w)               {fResponseMaker->SetRMMergeWeightFunction(w);}
 116         void            Make();
 117         void            MakeAU();       // test function, use with caution (09012014)
 118         void            Finish() {
 119             fOutputFile->cd();
 120             if(fRMSSpectrumIn)  fRMSSpectrumIn->Write();
 121             if(fRMSSpectrumOut) fRMSSpectrumOut->Write();
 122             if(fRMSRatio)       fRMSRatio->Write();
 123             fOutputFile->Close();}
 124         void            PostProcess(
 125                 TString def,
 126                 Int_t columns = 4,
 127                 TString in = "UnfoldedSpectra.root",
 128                 TString out = "ProcessedSpectra.root") const;
 129         Bool_t          SetRawInput (
 130                 TH2D* detectorResponse, // detector response matrix
 131                 TH1D* jetPtIn,          // in plane jet spectrum
 132                 TH1D* jetPtOut,         // out of plane jet spectrum
 133                 TH1D* dptIn,            // in plane delta pt distribution
 134                 TH1D* dptOut,           // out of plane delta pt distribution
 135                 Int_t eventCount = 0);  // event count (optional)
 136         // static const helper functions, mainly histogram manipulation
 137         static TH1D*    ResizeXaxisTH1D(TH1D* histo, Int_t low, Int_t up, TString suffix = "");
 138         static TH2D*    ResizeYaxisTH2D(TH2D* histo, TArrayD* x, TArrayD* y, TString suffix = "");
 139         static TH2D*    NormalizeTH2D(TH2D* histo, Bool_t noError = kTRUE);
 140         static TH1D*    RebinTH1D(TH1D* histo, TArrayD* bins, TString suffix = "", Bool_t kill = kTRUE);
 141         TH2D*           RebinTH2D(TH2D* histo, TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
 142         static TH2D*    MatrixMultiplication(TH2D* a, TH2D* b, TString name = "CombinedResponse");
 143         static TH1D*    NormalizeTH1D(TH1D* histo, Double_t scale = 1.);
 144         static TGraphErrors*    GetRatio(TH1 *h1 = 0x0, TH1* h2 = 0x0, TString name = "", Bool_t appendFit = kFALSE, Int_t xmax = -1);
 145         static TGraphErrors*    GetV2(TH1* h1 = 0x0, TH1* h2 = 0x0, Double_t r = .7, TString name = "");
 146         static void     WriteObject(TObject* object, TString suffix = "", Bool_t kill = kTRUE);
 147         static TH2D*    ConstructDPtResponseFromTH1D(TH1D* dpt, Bool_t AvoidRoundingError);
 148         static TH2D*    GetUnityResponse(TArrayD* binsTrue, TArrayD* binsRec, TString suffix = "");
 149         void            SaveConfiguration(Bool_t convergedIn, Bool_t convergedOut) const;
 150         static TMatrixD*        CalculatePearsonCoefficients(TMatrixD* covmat);
 151         static TH1D*    SmoothenPrior(TH1D* spectrum, TF1* function, Double_t min, Double_t max, Double_t start, Bool_t kill = kTRUE, Bool_t counts = kTRUE);
 152         // set styles
 153         static void     Style(TCanvas* c, TString style = "PEARSON");
 154         static void     Style(TVirtualPad* c, TString style = "SPECTRUM");
 155         static void     Style(TLegend* l);
 156         static void     Style(TH1* h, EColor col = kBlue, histoType = kEmpty);
 157         static void     Style(TGraph* h, EColor col = kBlue, histoType = kEmpty);
 158         static TLegend* AddLegend(TVirtualPad* p)       {return p->BuildLegend(.27, .61, .88, .88);}
 159         static void     SavePadToPDF(TVirtualPad* pad)  {pad->SaveAs(Form("%s.pdf", pad->GetName()));}
 160         // interface to AliUnfolding, not necessary but nice to have all parameters in one place
 161         static void     SetMinuitStepSize(Float_t s)    {AliUnfolding::SetMinuitStepSize(s);}
 162         static void     SetMinuitPrecision(Float_t s)   {AliUnfolding::SetMinuitPrecision(s);}
 163         static void     SetMinuitPrecision(Int_t i)     {AliUnfolding::SetMinuitMaxIterations(i);}
 164         static void     SetMinuitStrategy(Double_t s)   {AliUnfolding::SetMinuitStrategy(s);}
 165         static void     SetDebug(Int_t d)               {AliUnfolding::SetDebug(d);}
 166     private:
 167         Bool_t          PrepareForUnfolding();
 168         Bool_t          PrepareForUnfolding(Int_t low, Int_t up);
 169         TH1D*           GetPrior(                       const TH1D* measuredJetSpectrum,
 170                                                         const TH2D* resizedResponse,
 171                                                         const TH1D* kinematicEfficiency,
 172                                                         const TH1D* measuredJetSpectrumTrueBins,
 173                                                         const TString suffix,
 174                                                         const TH1D* jetFindingEfficiency);
 175         TH1D*           UnfoldWrapper(                  const TH1D* measuredJetSpectrum,
 176                                                         const TH2D* resizedResponse,
 177                                                         const TH1D* kinematicEfficiency,
 178                                                         const TH1D* measuredJetSpectrumTrueBins,
 179                                                         const TString suffix,
 180                                                         const TH1D* jetFindingEfficiency = 0x0);
 181         TH1D*           UnfoldSpectrumChi2(             const TH1D* measuredJetSpectrum,
 182                                                         const TH2D* resizedResponse,
 183                                                         const TH1D* kinematicEfficiency,
 184                                                         const TH1D* measuredJetSpectrumTrueBins,
 185                                                         const TString suffix,
 186                                                         const TH1D* jetFindingEfficiency = 0x0);
 187         TH1D*           UnfoldSpectrumSVD(              const TH1D* measuredJetSpectrum,
 188                                                         const TH2D* resizedResponse,
 189                                                         const TH1D* kinematicEfficiency,
 190                                                         const TH1D* measuredJetSpectrumTrueBins,
 191                                                         const TString suffix,
 192                                                         const TH1D* jetFindingEfficiency = 0x0);
 193         TH1D*           UnfoldSpectrumBayesianAli(      const TH1D* measuredJetSpectrum,
 194                                                         const TH2D* resizedResponse,
 195                                                         const TH1D* kinematicEfficiency,
 196                                                         const TH1D* measuredJetSpectrumTrueBins,
 197                                                         const TString suffix,
 198                                                         const TH1D* jetFindingEfficiency = 0x0);
 199         TH1D*           UnfoldSpectrumBayesian(         const TH1D* measuredJetSpectrum,
 200                                                         const TH2D* resizedResponse,
 201                                                         const TH1D* kinematicEfficiency,
 202                                                         const TH1D* measuredJetSpectrumTrueBins,
 203                                                         const TString suffix,
 204                                                         const TH1D* jetFindingEfficiency = 0x0);
 205         static void     ResetAliUnfolding();
 206         // give object a unique name via the 'protect heap' functions.
 207         // may seem redundant, but some internal functions of root (e.g.
 208         // ProjectionY()) check for existing objects by name and re-use them
 209         TH1D*           ProtectHeap(TH1D* protect, Bool_t kill = kTRUE, TString suffix = "");
 210         TH2D*           ProtectHeap(TH2D* protect, Bool_t kill = kTRUE, TString suffix = "");
 211         TGraphErrors*   ProtectHeap(TGraphErrors* protect, Bool_t kill = kTRUE, TString suffix = "");
 212         // members, accessible via setters
 213         AliAnaChargedJetResponseMaker*  fResponseMaker; // utility object
 214         TF1*                    fPower;                 // smoothening fit
 215         Bool_t                  fSaveFull;              // save all generated histograms to file
 216         TString                 fActiveString;          // identifier of active output
 217         TDirectoryFile*         fActiveDir;             // active directory
 218         TList*                  fInputList;             // input list
 219         Bool_t                  fRefreshInput;          // re-read the input (called automatically if input list changes)
 220         TString                 fOutputFileName;        // output file name
 221         TFile*                  fOutputFile;            // output file
 222         Int_t                   fCentralityBin;         // centrality bin
 223         TH2D*                   fDetectorResponse;      // detector response
 224         TH1D*                   fJetFindingEff;         // jet finding efficiency
 225         Double_t                fBetaIn;                // regularization strength, in plane unfolding
 226         Double_t                fBetaOut;               // regularization strength, out of plane unfoldign
 227         Int_t                   fBayesianIterIn;        // bayesian regularization parameter, in plane unfolding
 228         Int_t                   fBayesianIterOut;       // bayesian regularization parameter, out plane unfolding
 229         Float_t                 fBayesianSmoothIn;      // bayesian smoothening parameter (AliUnfolding)
 230         Float_t                 fBayesianSmoothOut;     // bayesian smoothening parameter (AliUnfolding)
 231         Bool_t                  fAvoidRoundingError;    // set dpt to zero for small values far from the diagonal
 232         unfoldingAlgorithm      fUnfoldingAlgorithm;    // algorithm used for unfolding
 233         prior                   fPrior;                 // prior for unfolding
 234         TArrayD*                fBinsTrue;              // pt true bins
 235         TArrayD*                fBinsRec;               // pt rec bins
 236         TArrayD*                fBinsTruePrior;         // holds true bins for the chi2 prior for SVD. setting this is optional
 237         TArrayD*                fBinsRecPrior;          // holds rec bins for the chi2 prior for SVD. setting this is optional
 238         Int_t                   fSVDRegIn;              // svd regularization (a good starting point is half of the number of bins)
 239         Int_t                   fSVDRegOut;             // svd regularization out of plane
 240         Bool_t                  fSVDToy;                // use toy to estimate coveriance matrix for SVD method
 241         Float_t                 fJetRadius;             // jet radius (for SVD toy)
 242         Int_t                   fEventCount;            // number of events
 243         Bool_t                  fNormalizeSpectra;      // normalize spectra to event count
 244         Bool_t                  fSmoothenPrior;         // smoothen the tail of the measured spectrum using a powerlaw fit
 245         Float_t                 fFitMin;                // lower bound of smoothening fit
 246         Float_t                 fFitMax;                // upper bound of smoothening fit
 247         Float_t                 fFitStart;              // from this value, use smoothening
 248         Bool_t                  fSmoothenCounts;        // fill smoothened spectrum with counts
 249         Bool_t                  fTestMode;              // unfold with unity response for testing
 250         Bool_t                  fRawInputProvided;      // input histograms provided, not read from file
 251         Double_t                fEventPlaneRes;         // event plane resolution for current centrality
 252         Bool_t                  fUseDetectorResponse;   // add detector response to unfolding
 253         Bool_t                  fUseDptResponse;        // add dpt response to unfolding
 254         Bool_t                  fTrainPower;            // don't clear the params of fPower for call to Make
 255                                                         // might give more stable results, but possibly introduces
 256                                                         // a bias / dependency on previous iterations
 257         Bool_t                  fDphiUnfolding;         // do the unfolding in in and out of plane orientation
 258         Bool_t                  fDphiDptUnfolding;      // do the unfolding in dphi and dpt bins (to fit v2)
 259         Bool_t                  fExLJDpt;               // exclude randon cones with leading jet
 260         // members, set internally
 261         TProfile*               fRMSSpectrumIn;         // rms of in plane spectra of converged unfoldings
 262         TProfile*               fRMSSpectrumOut;        // rms of out of plane spectra of converged unfoldings
 263         TProfile*               fRMSRatio;              // rms of ratio of converged unfolded spectra
 264         TProfile*               fRMSV2;                 // rms of v2 of converged unfolded spectra
 265         TH2D*                   fDeltaPtDeltaPhi;       // delta pt delta phi distribution
 266         TH2D*                   fJetPtDeltaPhi;         // jet pt delta phi distribution
 267         TH1D*                   fSpectrumIn;            // in plane jet pt spectrum
 268         TH1D*                   fSpectrumOut;           // out of plane jet pt spectrum
 269         TH1D*                   fDptInDist;             // in plane dpt distribution
 270         TH1D*                   fDptOutDist;            // out of plane dpt distribution
 271         TH2D*                   fDptIn;                 // in plane dpt matrix
 272         TH2D*                   fDptOut;                // out plane dpt matrix
 273         TH2D*                   fFullResponseIn;        // full response matrix, in plane
 274         TH2D*                   fFullResponseOut;       // full response matrix, out of plane
 275         // copy and assignment
 276         AliJetFlowTools(const AliJetFlowTools&);             // not implemented
 277         AliJetFlowTools& operator=(const AliJetFlowTools&);  // not implemented
 278 };
 279 #endif
 280 //_____________________________________________________________________________