PWGHF/hfe/macros/CorrectForEfficiencyPbPb.C

   1 void CorrectForEfficiencyPbPb(const char *filedata,const char *fileMC,const char *NameList,Bool_t smoothingon=kFALSE,Bool_t hadroncontaminationsubtracted=kFALSE,Bool_t unsetCorrelatedErrors=kTRUE);
   2 TList *GetResults(const Char_t *testfile,const Char_t *plus="");
   3 TList *GetQA(const Char_t *testfile,const Char_t *plus="");
   4 TObject* GetSpectrum(AliCFContainer *c, Int_t step);
   5 THnSparseF* GetHadronEffbyIPcut(TList *hfeqa);
   6 void CorrectFromTheWidth(TH1D *h1);
   7
   8 void CorrectForEfficiencyPbPb(const char *filedata,const char *fileMC,const char *NameList,Bool_t smoothingon,Bool_t hadroncontaminationsubtracted,Bool_t unsetCorrelatedErrors) {
   9
  10   ///////////////////////////////////////////////////////////////////////////////////////////////////////
  11   // Inputs of the macro:
  12   // filedata and fileMC: path of the output over data and MC
  13   // NameList: name of the list in this output ("HFE_Results" for example)
  14   // smoothingon: kTRUE means smoothing is used, should be kFALSE in general
  15   // hadroncontaminationsubtracted: kTRUE means hadron contamination is subtracted
  16   // unsetCorrelatedErrors: kTRUE means that the errors may be not properly calculated, nevertheless in 2D (centrality,pt) the proper calculation of the errors gave crazy results in the past, therefore kTRUE per default
  17   //
  18   // In the macro the centrality bins for the final corrected spectra are choosen
  19   // Look at:
  20   // spectrum->SetNCentralityBinAtTheEnd(2);
  21   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(0,2,0); // 0-10%
  22   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(7,9,1); // 60-80%
  23   // Per definition:
  24   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(0,1,0); // 0-5%
  25   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(1,2,0); // 5-10%
  26   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(2,3,0); // 10-20%
  27   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(3,4,0); // 20-30%
  28   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(4,5,0); // 30-40%
  29   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(5,6,0); // 40-50%
  30   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(6,7,0); // 50-60%
  31   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(7,8,0); // 60-70%
  32   // spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(8,9,0); // 70-80%
  33   //
  34   // Outputs:
  35   //  TFile*         finalSpectrum.root
  36   //
  37   //  Minimum-bias results not relevant for the analysis
  38   //
  39   //  TGraphErrors     UnfoldingCorrectedSpectrum;1
  40   //  TGraphErrors     AlltogetherSpectrum;1
  41   //  TH1D     RatioUnfoldingAlltogetherSpectrum;1
  42   //  THnSparseT<TArrayF>      UnfoldingCorrectedNotNormalizedSpectrum;1       step1 projection centralitypt
  43   //  AliCFDataGrid    AlltogetherCorrectedNotNormalizedSpectrum;1     dataGrid
  44   //
  45   // Results:
  46   //
  47   // Unfolded results: can be used for minimum-bias and enhanced sample
  48   //
  49   // TH1D     Unfolded_Notnormalized_centrality_bin_[0_2[;1    // Not normalized to the number of events
  50   // TGraphErrors     Unfolded_normalized_centrality_bin_[0_2[;1  // Normalized to the number of events ---> usual output
  51   //
  52   // Direct correction (ptesd/ptMC): can be used ONLY for minimum-bias sample
  53   //
  54   // TH1D     Dirrectcorrected_Notnormalized_centrality_bin_[0_2[;1   // Not normalized to the number of events
  55   // TGraphErrors     Dirrectedcorrected_normalized_centrality_bin_[0_2[;1 // Normalized to the number of events --> usual output
  56   //
  57   // Important:
  58   //
  59   // TPC pid efficiency is not included --> results have to be multiplied by 2
  60   //
  61   // For minimum-bias Dirrectcorrected_normalized_centrality_bin_[0_2[ and Unfolded_normalized_centrality_bin_[0_2[ should be similar
  62   // Not for enhance sample: Dirrectcorrected is wrong.
  63   //
  64   // Will produce finalSpectrum.root with results inside
  65   // TGraphErrors UnfoldingCorrectedSpectrum -> unfolding procedure (doesn't work in 2D with trunk)
  66   // TGraphErrors AlltogetherSpectrum -> corrected spectrum with other method
  67   // TH1D RatioUnfoldingAlltogetherSpectrum
  68   // THnSparseF UnfoldingCorrectedNotNormalizedSpectrum -> unfolding procedure not yet normalized
  69   // AliCFDataGrid AlltogetherCorrectedNotNormalizedSpectrum -> other method not yet normalized
  70   ///////////////////////////////////////////////////////////////////////////////////////////////////////
  71
  72   gStyle->SetPalette(1);
  73   gStyle->SetOptStat(1111);
  74   gStyle->SetPadBorderMode(0);
  75   gStyle->SetCanvasColor(10);
  76   gStyle->SetPadLeftMargin(0.13);
  77   gStyle->SetPadRightMargin(0.13);
  78
  79   /////////////////////
  80   // Take the stuff
  81   /////////////////////
  82
  83   TList *resultsdata = GetResults(filedata,NameList);
  84   if(!resultsdata){
  85     printf("No output objects for data: Calculation will terminate here\n");
  86     return;
  87   }
  88
  89   ///////////////////////////////////////////////////////////
  90   // Event container for the normalization to CINB1
  91   // Normalize to the number events after event selection
  92   ////////////////////////////////////////////////////////////
  93
  94   AliCFContainer *containerdata = (AliCFContainer *) resultsdata->FindObject("eventContainer");
  95   if(!containerdata) {
  96     printf("No container \n");
  97     return;
  98   }
  99   AliCFDataGrid *dataGrid = (AliCFDataGrid *) GetSpectrum(containerdata,AliHFEcuts::kEventStepReconstructed);
 100   THnSparseF *eventcontainer = (THnSparseF*) dataGrid->GetGrid();
 101   TAxis *cenaxisb = eventcontainer->GetAxis(2);
 102   Int_t nbbin = cenaxisb->GetNbins();
 103   printf("There are %d centrality bins!!!!!!!!!!!\n",nbbin);
 104   Double_t *nbeventstotal = new Double_t[nbbin];
 105
 106   for(Int_t binc = 0; binc < nbbin; binc++){
 107     cenaxisb->SetRange(binc+1,binc+1);
 108     TH1D *spectrum1Dc = (TH1D *) dataGrid->Project(0);
 109     nbeventstotal[binc] = spectrum1Dc->Integral();
 110   }
 111
 112   Float_t *numberofentries = new Float_t[nbbin];
 113   Float_t numberofentriessum = 0.0;
 114   for(Int_t binc = 0; binc < nbbin; binc++) {
 115     numberofentries[binc] = nbeventstotal[binc];
 116     numberofentriessum += numberofentries[binc];
 117     printf("Number %f for centrality bin %d\n",numberofentries[binc],binc);
 118   }
 119
 120   //////////////////////////////
 121   // Take more stuff
 122   ///////////////////////////////
 123
 124
 125   TList *resultsmc = GetResults(fileMC,NameList);
 126   if(!resultsmc){
 127     printf("No output objects for mc: Calculation will terminate here\n");
 128     return;
 129   }
 130
 131   AliHFEcontainer *datahfecontainer = (AliHFEcontainer *) resultsdata->FindObject("trackContainer");
 132   if(!datahfecontainer) {
 133     printf("No container for data \n");
 134     return;
 135   }
 136   AliCFContainer *sumcontainer = datahfecontainer->GetCFContainer("recTrackContReco");
 137
 138
 139   /////////////////////////////
 140   // Check number of events
 141   /////////////////////////////
 142
 143
 144   Int_t numberOfEventsafterallcuts = (Int_t) datahfecontainer->GetNumberOfEvents();
 145
 146   AliHFEcontainer *mchfecontainer = (AliHFEcontainer *) resultsmc->FindObject("trackContainer");
 147   if(!mchfecontainer) {
 148     printf("No mc container \n");
 149     return;
 150   }
 151
 152
 153   //////////////
 154   // Correct
 155   /////////////
 156
 157   AliHFEspectrum *spectrum = new AliHFEspectrum("HFEspectrum");
 158   spectrum->SetNbDimensions(1);
 159   // If you want to correct positive (0) or negative (1) separately
 160   //spectrum->SetChargeChoosen(0);
 161   spectrum->SetDebugLevel(1);
 162   // Give the number of events per centrality bins
 163   for(Int_t binc = 0; binc < nbbin; binc++) {
 164     spectrum->SetNumberOfEvents((Int_t)numberofentries[binc],binc);
 165   }
 166   // Calculation of the errors
 167   spectrum->SetUnSetCorrelatedErrors(unsetCorrelatedErrors);
 168   /////////////////////////////////////////////////////////////////////////
 169   // Number of centrality bin we want for the final corrected spectra
 170   ///////////////////////////////////////////////////////////////////////
 171   spectrum->SetNCentralityBinAtTheEnd(2);
 172   spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(0,2,0); // 0-10%
 173   spectrum->SetLowHighBoundaryCentralityBinAtTheEnd(7,9,1); // 60-80%
 174
 175   spectrum->SetDumpToFile(kTRUE);
 176   spectrum->SetSmoothing(smoothingon);
 177   // True step in MC (events in range +- 10 cm and no pile up)
 178   spectrum->SetMCTruthStep(AliHFEcuts::kStepMCGeneratedEventCut);
 179   // Step where we correct from MC (tracking + TOF)
 180   spectrum->SetMCEffStep(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD +1);
 181   // Step from data we correct for
 182   spectrum->SetStepToCorrect(AliHFEcuts::kStepHFEcutsTRD + 2);
 183   // PbPb flag
 184   spectrum->SetPbPbAnalysis(kTRUE);
 185
 186   // Give everything (data container, mc container and V0 data container)
 187   spectrum->Init(datahfecontainer,mchfecontainer);
 188
 189   // kTRUE means subtract hadron background, kFALSE means do not subtract hadron background
 190   spectrum->Correct(hadroncontaminationsubtracted);
 191
 192
 193
 194
 195
 196 }
 197
 198 //_____________________________________________________________________
 199 TList *GetResults(const Char_t *testfile,const Char_t *plus){
 200   //
 201   // read output
 202   //
 203   TFile *f = TFile::Open(testfile);
 204   if(!f || f->IsZombie()){
 205     printf("File not readable\n");
 206     return 0x0;
 207   }
 208   TString name(plus);
 209   printf("Name of TList %s\n",(const char*)name);
 210   TList *l = dynamic_cast<TList *>(f->Get((const char*)name));
 211   if(!l){
 212     printf("Results list was not found\n");
 213     f->Close(); delete f;
 214     return 0x0;
 215   }
 216   TList *returnlist = dynamic_cast<TList *>(l->Clone());
 217   f->Close(); delete f;
 218   return returnlist;
 219 }
 220
 221
 222 //_____________________________________________________________________
 223 TList *GetQA(const Char_t *testfile,const Char_t *plus){
 224   //
 225   // read output
 226   //
 227   TFile *f = TFile::Open(testfile);
 228   if(!f || f->IsZombie()){
 229     printf("File not readable\n");
 230     return 0x0;
 231   }
 232   TString name("HFE_QA");
 233   name += plus;
 234   printf("Name of TList %s\n",(const char*)name);
 235   TList *l = dynamic_cast<TList *>(f->Get((const char*)name));
 236   if(!l){
 237     printf("QA list was not found\n");
 238     f->Close(); delete f;
 239     return 0x0;
 240   }
 241   TList *returnlist = dynamic_cast<TList *>(l->Clone());
 242   f->Close(); delete f;
 243   return returnlist;
 244 }
 245
 246 //_____________________________________________________________________
 247 THnSparseF* GetHadronEffbyIPcut(TList *hfeqa){
 248
 249
 250   // get hadron reduction factors due to the IP cuts
 251   TList* tl=(TList*)hfeqa->FindObject("list_TaskQA");
 252   TH1F* hbefore=(TH1F*)tl->FindObject("hadronsBeforeIPcut");
 253   TH1F* hafter=(TH1F*)tl->FindObject("hadronsAfterIPcut");
 254   TH1F* hreduction= (TH1F*)hafter->Clone("hreduction");
 255   hbefore->Sumw2();
 256   hafter->Sumw2();
 257   hreduction->Sumw2();
 258   hreduction->Divide(hbefore);
 259
 260   Double_t* binEdges[0];
 261   Int_t hnBin = hreduction->GetNbinsX();
 262   Int_t nBin[1] = {hnBin};
 263
 264   for(int i=0; i<nBin[0]; i++){
 265     binEdges[0][i] = hreduction->GetBinLowEdge(i+1);
 266   }
 267   binEdges[0][nBin[0]] = hreduction->GetBinLowEdge(nBin[0]) + hreduction->GetBinWidth(nBin[0]);
 268
 269   THnSparseF* hsreduction = new THnSparseF("hadroncontamin", "hadroncontamin; pt[GeV/c]", 1, nBin);
 270   hsreduction->SetBinEdges(0, binEdges[0]);
 271
 272   Double_t dataE[1];
 273   Double_t yval;
 274   for(int i=0; i<nBin[0]; i++){
 275     dataE[0]=hreduction->GetBinCenter(i+1);
 276     yval=hreduction->GetBinContent(i+1);
 277     hsreduction->Fill(dataE,yval);
 278   }
 279
 280   Int_t* ibins;
 281   ibins = new Int_t[nBin[0] + 1];
 282   hsreduction->SetBinError(ibins,0);
 283
 284
 285   return hsreduction;
 286 }
 287
 288 //_________________________________________________________________________
 289 TObject* GetSpectrum(AliCFContainer *c, Int_t step) {
 290   AliCFDataGrid* data = new AliCFDataGrid("data","",*c,step);
 291   //data->SetMeasured(step);
 292   return data;
 293 }
 294 //___________________________________________________________________________
 295 void CorrectFromTheWidth(TH1D *h1) {
 296   //
 297   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
 298   //
 299
 300   TAxis *axis = h1->GetXaxis();
 301   Int_t nbinX = h1->GetNbinsX();
 302
 303   for(Int_t i = 1; i <= nbinX; i++) {
 304
 305     Double_t width = axis->GetBinWidth(i);
 306     Double_t content = h1->GetBinContent(i);
 307     Double_t error = h1->GetBinError(i);
 308     h1->SetBinContent(i,content/width);
 309     h1->SetBinError(i,error/width);
 310   }
 311
 312 }