PWGHF/vertexingHF/charmFlow/AliEventPlaneResolutionHandler.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 #include <TMath.h>
  17 #include <TH1F.h>
  18 #include <TFile.h>
  19 #include "AliVertexingHFUtils.h"
  20 #include "AliEventPlaneResolutionHandler.h"
  21 #include "AliLog.h"
  22
  23 /* $Id$ */
  24
  25 ///////////////////////////////////////////////////////////////////
  26 //                                                               //
  27 // Class with functions useful for different D2H analyses        //
  28 // - event plane resolution                                      //
  29 // - <pt> calculation with side band subtraction                 //
  30 // - tracklet multiplicity calculation                            //
  31 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  32 //                                                               //
  33 ///////////////////////////////////////////////////////////////////
  34
  35 ClassImp(AliEventPlaneResolutionHandler)
  36
  37 //______________________________________________________________________
  38 AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler():TObject(),
  39   fEventPlane(kTPCPosEta),
  40   fResolOption(kTwoEtaSub),
  41   fMinCent(30.),
  42   fMaxCent(50.),
  43   fCorrHistoName1(""),
  44   fCorrHistoName2(""),
  45   fCorrHistoName3(""),
  46   fNsubevents(2),
  47   fExtrapToFull(kFALSE),
  48   fWeight(0),
  49   fHistoAB(0x0),
  50   fHistoBC(0x0),
  51   fHistoAC(0x0),
  52   fRootFileName("$ALICE_ROOT/PWGHF/vertexingHF/charmFlow/EventPlaneResolutionHistos.root")
  53 {
  54   // Default contructor
  55   InitializeWeights();
  56 }
  57 //______________________________________________________________________
  58 AliEventPlaneResolutionHandler::AliEventPlaneResolutionHandler(TString filename):TObject(),
  59   fEventPlane(kTPCPosEta),
  60   fResolOption(kTwoEtaSub),
  61   fMinCent(30.),
  62   fMaxCent(50.),
  63   fCorrHistoName1(""),
  64   fCorrHistoName2(""),
  65   fCorrHistoName3(""),
  66   fNsubevents(2),
  67   fExtrapToFull(kFALSE),
  68   fWeight(0),
  69   fHistoAB(0x0),
  70   fHistoBC(0x0),
  71   fHistoAC(0x0),
  72   fRootFileName(filename.Data())
  73 {
  74   // Standard contructor
  75   InitializeWeights();
  76 }
  77 //______________________________________________________________________
  78 Double_t AliEventPlaneResolutionHandler::GetEventPlaneResolution(Double_t minCent, Double_t maxCent){
  79   // method to compute the event plane resolution starting from correlations
  80
  81   TDirectory *current = gDirectory;
  82
  83   TFile* inputFile=new TFile(fRootFileName.Data());
  84   if(!inputFile || (inputFile && !inputFile->IsOpen())){
  85     printf("Root file not opened --- Return dummy resolution value\n");
  86     return -1.;
  87   }
  88   Bool_t isOk=SetHistoNames();
  89   if(!isOk){
  90     printf("Return dummy resolution value\n");
  91     return -1.;
  92   }
  93
  94   Int_t minCentrTimesTen=(Int_t)(minCent*10);
  95   Int_t maxCentrTimesTen=(Int_t)(maxCent*10);
  96   Int_t ncBin=minCentrTimesTen/25;
  97   Int_t whichWei=fWeight;
  98   if(minCent>50.5 || maxCent>50.5) fWeight=0;
  99   Double_t weight=-1.;
 100   if(fWeight==1) weight=fNcoll[ncBin];
 101   // protection: weights from yield computed for 30-50, 10-30 and 0-10 centralities
 102   // cannot be used on wider ranges
 103   else if(fWeight==2 && minCent>29.5 && maxCent<50.5)  weight=fYield24[ncBin];
 104   else if(fWeight==3 && minCent>29.5 && maxCent<50.5)  weight=fYield46[ncBin];
 105   else if(fWeight==4 && minCent>29.5 && maxCent<50.5)  weight=fYield612[ncBin];
 106   else if(fWeight==2 && minCent>9.5 && maxCent<30.5)  weight=fYield24[ncBin];
 107   else if(fWeight==3 && minCent>9.5 && maxCent<30.5)  weight=fYield46[ncBin];
 108   else if(fWeight==4 && minCent>9.5 && maxCent<30.5)  weight=fYield612[ncBin];
 109   else if(fWeight==2 && maxCent<10.5)  weight=fYield24[ncBin];
 110   else if(fWeight==3 && maxCent<10.5)  weight=fYield46[ncBin];
 111   else if(fWeight==4 && maxCent<10.5)  weight=fYield612[ncBin];
 112   TH1F* hevpls2=0x0;
 113   TH1F* hevpls3=0x0;
 114   TH1F* hevpls1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,minCentrTimesTen+25));
 115   hevpls1->SetName(Form("%sCentr%d_%d",fCorrHistoName1.Data(),minCentrTimesTen,maxCentrTimesTen));
 116   if(weight>0.) hevpls1->Scale(weight);
 117   if(fNsubevents==3){
 118     hevpls2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,minCentrTimesTen+25));
 119     hevpls2->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
 120     hevpls3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),minCentrTimesTen,minCentrTimesTen+25));
 121     hevpls3->SetName(Form("%sCentr%d_%d",fCorrHistoName2.Data(),minCentrTimesTen,maxCentrTimesTen));
 122     if(weight>0.){
 123       hevpls2->Scale(weight);
 124       hevpls3->Scale(weight);
 125     }
 126   }
 127   for(Int_t iCentr=minCentrTimesTen+25; iCentr<maxCentrTimesTen; iCentr+=25){
 128     ncBin=iCentr/25;
 129     weight=-1.;
 130     if(fWeight==1) weight=fNcoll[ncBin];
 131     else if(fWeight==2 && minCent>29.5 && maxCent<50.5)  weight=fYield24[ncBin];
 132     else if(fWeight==3 && minCent>29.5 && maxCent<50.5)  weight=fYield46[ncBin];
 133     else if(fWeight==4 && minCent>29.5 && maxCent<50.5)  weight=fYield612[ncBin];
 134     else if(fWeight==2 && minCent>9.5 && maxCent<30.5)  weight=fYield24[ncBin];
 135     else if(fWeight==3 && minCent>9.5 && maxCent<30.5)  weight=fYield46[ncBin];
 136     else if(fWeight==4 && minCent>9.5 && maxCent<30.5)  weight=fYield612[ncBin];
 137     else if(fWeight==2 && maxCent<10.5)  weight=fYield24[ncBin];
 138     else if(fWeight==3 && maxCent<10.5)  weight=fYield46[ncBin];
 139     else if(fWeight==4 && maxCent<10.5)  weight=fYield612[ncBin];
 140     TH1F* htmp1=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName1.Data(),iCentr,iCentr+25));
 141     if(weight>0.) htmp1->Scale(weight);
 142     hevpls1->Add(htmp1);
 143     if(fNsubevents==3){
 144       TH1F* htmp2=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName2.Data(),iCentr,iCentr+25));
 145       TH1F* htmp3=(TH1F*)inputFile->Get(Form("%sCentr%d_%d",fCorrHistoName3.Data(),iCentr,iCentr+25));
 146       if(weight>0.){
 147         htmp2->Scale(weight);
 148         htmp3->Scale(weight);
 149       }
 150       hevpls2->Add(htmp2);
 151       hevpls3->Add(htmp3);
 152     }
 153   }
 154   fWeight=whichWei;
 155   Double_t resol=1.;
 156   if(fNsubevents==2){
 157     if(fExtrapToFull) resol=AliVertexingHFUtils::GetFullEvResol(hevpls1);
 158     else resol=AliVertexingHFUtils::GetSubEvResol(hevpls1);
 159     current->cd();
 160     if(fHistoAB) delete fHistoAB;
 161     fHistoAB=(TH1F*)hevpls1->Clone("hCorrelSubAB");
 162   }else if(fNsubevents==3){
 163     Double_t correl1=hevpls1->GetMean();
 164     Double_t correl2=hevpls2->GetMean();
 165     Double_t correl3=hevpls3->GetMean();
 166     resol=TMath::Sqrt(correl1*correl2/correl3);
 167     current->cd();
 168     if(fHistoAB) delete fHistoAB;
 169     fHistoAB=(TH1F*)hevpls1->Clone("hCorrelSubAB");
 170     if(fHistoBC) delete fHistoBC;
 171     fHistoBC=(TH1F*)hevpls2->Clone("hCorrelSubBC");
 172     if(fHistoAC) delete fHistoAC;
 173     fHistoAC=(TH1F*)hevpls3->Clone("hCorrelSubAC");
 174   }
 175   inputFile->Close();
 176   return resol;
 177 }
 178 //______________________________________________________________________
 179 Bool_t AliEventPlaneResolutionHandler::SetHistoNames(){
 180   // set names for histograms to be read
 181   if(fEventPlane==kTPCPosEta && fResolOption==kTwoRandSub){
 182     fCorrHistoName1="hCorrelRandSubTPCPosEta";
 183     fNsubevents=2;
 184     fExtrapToFull=kTRUE;
 185   }else if(fEventPlane==kTPCPosEta && fResolOption==kTwoChargeSub){
 186     fCorrHistoName1="hCorrelChargeSubTPCPosEta";
 187     fNsubevents=2;
 188     fExtrapToFull=kTRUE;
 189   }else if(fEventPlane==kTPCPosEta && fResolOption==kTwoEtaSub){
 190     fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
 191     fNsubevents=2;
 192     fExtrapToFull=kFALSE;
 193   }else if(fEventPlane==kTPCPosEta && fResolOption==kThreeSub){
 194     fCorrHistoName1="hCorrelTPCPosEtaV0A";
 195     fCorrHistoName2="hCorrelTPCPosEtaV0C";
 196     fCorrHistoName3="hCorrelV0AV0C";
 197     fNsubevents=3;
 198     fExtrapToFull=kFALSE;
 199   }else if(fEventPlane==kTPCFullEta && fResolOption==kTwoRandSub){
 200     fCorrHistoName1="hCorrelRandSubTPCFullEta";
 201     fNsubevents=2;
 202     fExtrapToFull=kTRUE;
 203   }else if(fEventPlane==kTPCFullEta && fResolOption==kTwoEtaSub){
 204     fCorrHistoName1="hCorrelTPCEtaPosTPCEtaNeg";
 205     fNsubevents=2;
 206     fExtrapToFull=kTRUE;
 207   }else if(fEventPlane==kTPCFullEta && fResolOption==kThreeSub){
 208     fCorrHistoName1="hCorrelTPCFullEtaV0A";
 209     fCorrHistoName2="hCorrelTPCFullEtaV0C";
 210     fCorrHistoName3="hCorrelV0AV0C";
 211     fNsubevents=3;
 212     fExtrapToFull=kFALSE;
 213   }else if(fEventPlane==kVZERO && fResolOption==kThreeSub){
 214     fCorrHistoName1="hCorrelTPCEtaPosV0";
 215     fCorrHistoName2="hCorrelTPCEtaNegV0";
 216     fCorrHistoName3="hCorrelTPCEtaPosTPCEtaNeg";
 217     fNsubevents=3;
 218     fExtrapToFull=kFALSE;
 219   }else if(fEventPlane==kVZERO && fResolOption==kThreeSubTPCGap){
 220     fCorrHistoName1="hCorrelTPCEtaGt02V0";
 221     fCorrHistoName2="hCorrelTPCEtaLt02V0";
 222     fCorrHistoName3="hCorrelTPCEtaGt02TPCEtaLt02";
 223     fNsubevents=3;
 224     fExtrapToFull=kFALSE;
 225   }else if(fEventPlane==kVZEROA && fResolOption==kThreeSub){
 226     fCorrHistoName1="hCorrelTPCFullEtaV0A";
 227     fCorrHistoName2="hCorrelV0AV0C";
 228     fCorrHistoName3="hCorrelTPCFullEtaV0C";
 229     fNsubevents=3;
 230     fExtrapToFull=kFALSE;
 231   }else if(fEventPlane==kVZEROC && fResolOption==kThreeSub){
 232     fCorrHistoName1="hCorrelTPCFullEtaV0C";
 233     fCorrHistoName2="hCorrelV0AV0C";
 234     fCorrHistoName3="hCorrelTPCFullEtaV0A";
 235     fNsubevents=3;
 236     fExtrapToFull=kFALSE;
 237   }else{
 238     printf("Not possible to compute the resolution with the required option\n");
 239     return kFALSE;
 240   }
 241   return kTRUE;
 242 }
 243 //______________________________________________________________________
 244 void AliEventPlaneResolutionHandler::InitializeWeights(){
 245   // Ncoll values in 2.5% wide classes
 246   Double_t ncoll[20]={1790.77,1578.44,1394.82,1236.17,1095.08,
 247                       969.86,859.571,759.959,669.648,589.588,
 248                       516.039,451.409,392.853,340.493,294.426,
 249                       252.385,215.484,183.284,155.101,130.963};
 250
 251   Double_t y24[20]={544.155,491.803,533.516,449.618,
 252                     213.097,179.483,194.977,156.73,118.84,152.306,92.7,122.538,
 253                     36.64,33.01,27.79,22.57,17.90,15.25,15.03,8.57};
 254
 255   Double_t y46[20]={158.386,129.294,135.683,112.203,
 256                     59.425,61.506,65.779,43.521,35.530,41.014,26.744,33.593,
 257                     12.24,12.32,10.44,8.55,7.20,5.85,5.80,4.37};
 258
 259   Double_t y612[20]={76.010,85.944,99.105,84.779,
 260                      45.404,25.059,36.734,31.177,31.587,46.088,19.367,39.945,
 261                      7.31,7.55,6.47,5.38,4.68,4.86,3.10,4.78};
 262
 263   for(Int_t i=0; i<20; i++){
 264     fNcoll[i]=ncoll[i];
 265     fYield24[i]=y24[i];
 266     fYield46[i]=y46[i];
 267     fYield612[i]=y612[i];
 268   }
 269 }