PWGLF/FORWARD/analysis2/scripts/UnfoldMult.C

   1 /**
   2  * @file   Unfolding.C
   3  * @author Valentina Zaccolo
   4  * @date   Tue Mar 05 12:56:56 2013
   5  *
   6  * @brief  Bayesian Method Unfolding Script
   7  *
   8  */
   9 #include <iostream>
  10 #include <TList.h>
  11 #include <TFile.h>
  12 #include <TObject.h>
  13 #include <TDirectory.h>
  14 #include "TH1D.h"
  15 #include "TH2D.h"
  16 #include <TROOT.h>
  17 #include <TSystem.h>
  18 #include "RooUnfoldBayes.h"
  19 #include "RooUnfoldResponse.h"
  20 #include <TMath.h>
  21
  22 /**
  23  * Get an object from a directory
  24  *
  25  * @param d     Directory
  26  * @param name  Name of object
  27  *
  28  * @return Pointer to object or null if not found
  29  */
  30 TObject* GetObject(TDirectory* d, const char* name)
  31 {
  32   TObject* o = d->Get(name);
  33   if (!o) {
  34     Error("GetCollection", "%s not found in %s", name, d->GetName());
  35     return 0;
  36   }
  37   return o;
  38 }
  39 /**
  40  * Get an object from a collection
  41  *
  42  * @param d     Collection
  43  * @param name  Name of object
  44  *
  45  * @return Pointer to object or null if not found
  46  */
  47 TObject* GetObject(TCollection* d, const char* name)
  48 {
  49   TObject* o = d->FindObject(name);
  50   if (!o) {
  51     Error("GetCollection", "%s not found in %s", name, d->GetName());
  52     d->ls();
  53     return 0;
  54   }
  55   return o;
  56 }
  57 /**
  58  * Get a collection from a directory
  59  *
  60  * @param d     Parent directory
  61  * @param name  Name of collection
  62  *
  63  * @return Pointer to collection or null
  64  */
  65 TCollection* GetCollection(TDirectory* d, const char* name)
  66 {
  67   return static_cast<TCollection*>(GetObject(d, name));
  68 }
  69 /**
  70  * Get a collection from a collection
  71  *
  72  * @param d     Parent collection
  73  * @param name  Name of collection
  74  *
  75  * @return Pointer to collection or null
  76  */
  77 TCollection* GetCollection(TCollection* d, const char* name)
  78 {
  79   return static_cast<TCollection*>(GetObject(d, name));
  80 }
  81 /**
  82  * Get a 1D-histogram from a directory
  83  *
  84  * @param d     Parent directory
  85  * @param name  Name of histogram
  86  *
  87  * @return Pointer to histogram or null
  88  */
  89 TH1* GetH1(TDirectory* d, const char* name)
  90 {
  91   return static_cast<TH1*>(GetObject(d, name));
  92 }
  93 /**
  94  * Get a 1D-histogram from a collection
  95  *
  96  * @param d     Parent collectoin
  97  * @param name  Name of histogram
  98  *
  99  * @return Pointer to histogram or null
 100  */
 101 TH1* GetH1(TCollection* d, const char* name)
 102 {
 103   return static_cast<TH1*>(GetObject(d, name));
 104 }
 105 /**
 106  * Get a 2D-histogram from a directory
 107  *
 108  * @param d     Parent directory
 109  * @param name  Name of histogram
 110  *
 111  * @return Pointer to histogram or null
 112  */
 113 TH2* GetH2(TDirectory* d, const char* name)
 114 {
 115   return static_cast<TH2*>(GetObject(d, name));
 116 }
 117 /**
 118  * Get a 2D-histogram from a collection
 119  *
 120  * @param d     Parent collection
 121  * @param name  Name of histogram
 122  *
 123  * @return Pointer to histogram or null
 124  */
 125 TH2* GetH2(TCollection* d, const char* name)
 126 {
 127   return static_cast<TH2*>(GetObject(d, name));
 128 }
 129 /**
 130  * Open a file
 131  *
 132  * @param filename      File name
 133  * @param readNotWrite  If true, open read-only
 134  *
 135  * @return Pointer to file, or null
 136  */
 137 TFile* OpenFile(const char* filename, Bool_t readNotWrite)
 138 {
 139   TFile* ret = TFile::Open(filename, readNotWrite ? "READ" : "RECREATE");
 140   if (!ret) {
 141     Error("OpenFile", "Failed to open the file \"%s\"", filename);
 142     return 0;
 143   }
 144   return ret;
 145 }
 146 /** A type definition to make it easier to switch to TDirectory structure */
 147 typedef TCollection Dir;
 148
 149 /**
 150  * Get histograms for @f$\eta@f$-bin from @a lowLim to @a highLim and
 151  * process it .
 152  *
 153  * @param lowLim   Lower @f$\eta@f$ limit
 154  * @param highLim  Upper @f$\eta@f$ limit
 155  * @param resp     Directory containing response matrices
 156  * @param data     Directory containing data histogram
 157  * @param out      Output directory
 158  */
 159 void GetHistos(Double_t    lowLim,
 160                Double_t    highLim,
 161                Dir*        resp,
 162                Dir*        data,
 163                TDirectory* out);
 164 /**
 165  * Do the actual unfolding and corrections
 166  *
 167  * @param response   Response matrix
 168  * @param measured   Measured distribution
 169  * @param mcHist     MC truth distribution
 170  * @param esdHist    MC distribution after ESD event selection
 171  */
 172 void ProcessUnfold(TH2* response, TH1* measured, TH1* MCHist, TH1* ESDHist);
 173 /**
 174  * Caclulate the trigger bias correction
 175  *
 176  * @param mcHist     MC truth distribution
 177  * @param esdHist    MC distribution after ESD event selection
 178  *
 179  * @return  Histogram of trigger bias correction
 180  */
 181 TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist);
 182 /**
 183  * Apply the trigger bias correction to the unfolded result
 184  *
 185  * @param unfoldBefore Unfolded measured distribution
 186  * @param triggerBias  Trigger bias correction
 187  * @param mcHist       MC truth distribution
 188  */
 189 void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1* mcHist);
 190
 191
 192 /**
 193  * Main Loop to call eta range
 194  *
 195  * @param respFile File containing response matrices
 196  * @param dataFile File containing measured distributions
 197  * @param outFile  Output file
 198  */
 199 void UnfoldMult(const char* respFile="forward_response.root",
 200                 const char* dataFile="forward_multiplicy.root",
 201                 const char* outFile="unfolded.root")
 202 {
 203   TFile* output = OpenFile(outFile,false);
 204   TFile* resp   = OpenFile(respFile, true);
 205   TFile* data   = OpenFile(dataFile, true);
 206   if (!output || !resp || !data) return;
 207
 208   Dir* topResp  = GetCollection(resp, "ResponseMatricesSums");
 209   Dir* topData  = GetCollection(data, "MultSums");
 210   if (!topData || !topResp) return;
 211
 212   Double_t limits[] = {5.1, 3.4, 3.0, 2.4, 2.0, 1.5, 1.0, 0.5, 0. };
 213   Double_t* limit = limits;
 214
 215   while ((*limit) > 0.1) {
 216     if ((*limit) >5.) GetHistos(-3.4,      (*limit), topResp, topData, output);
 217     else              GetHistos(-(*limit), (*limit), topResp, topData, output);
 218     limit++;
 219     output->cd();
 220   }
 221   output->Close();
 222   resp->Close();
 223   data->Close();
 224 }
 225
 226 //____________________________________________________________________
 227 void GetHistos(Double_t    lowLim,
 228                Double_t    highLim,
 229                Dir*        topResp,
 230                Dir*        topData,
 231                TDirectory* out)
 232 {
 233   // Format the name of the bin
 234   TString sLow( TString::Format("%+03d",int(10*lowLim)));
 235   TString sHigh(TString::Format("%+03d",int(10*highLim)));
 236   sLow .ReplaceAll("+", "plus");
 237   sLow .ReplaceAll("-", "minus");
 238   sHigh.ReplaceAll("+", "plus");
 239   sHigh.ReplaceAll("-", "minus");
 240   TString     name = TString::Format("%s_%s", sLow.Data(), sHigh.Data());
 241   TDirectory* dir = out->mkdir(name);
 242   dir->cd();
 243
 244   // Get our collections
 245   Dir* listResp = GetCollection(topResp, name);
 246   Dir* listData = GetCollection(topData,name);
 247   if(!listResp || !listData) return;
 248
 249   // Get the MC histograms
 250   TH2* response = GetH2(listResp, "responseMatrix");
 251   TH1* mcHist   = GetH1(listResp, "fMCNSD");
 252   TH1* esdHist  = GetH1(listResp, "fESDNSD");
 253   if (!response || !mcHist || !esdHist) return;
 254
 255   // Get the data histogram
 256   TH1* measured = GetH1(listData, "mult");
 257   if(!measured) return;
 258
 259   // Now do the unfolding
 260   ProcessUnfold(response, measured, mcHist, esdHist);
 261 }
 262
 263
 264
 265 //____________________________________________________________________
 266 void ProcessUnfold(TH2* response, TH1* measured, TH1* mcHist, TH1* esdHist)
 267 {
 268   Int_t nX   = response->GetNbinsX();
 269   Int_t nY   = response->GetNbinsY();
 270   TH1* projY = static_cast<TH1*>(response->ProjectionY("projY",1,nX,""));
 271   TH1* projX = static_cast<TH1*>(response->ProjectionX("projX",1,nY,""));
 272   projX->SetDirectory(0);
 273   projY->SetDirectory(0);
 274
 275   TH2* responseTranspose = static_cast<TH2*>(response->Clone("response"));
 276   for(Int_t i = 1; i <= nX; i++) {
 277     for(Int_t j = 1; j <= nY; j++) {
 278       responseTranspose->SetBinContent(j,i, response->GetBinContent(i,j));
 279       responseTranspose->SetBinError(j,i, response->GetBinError(i,j));
 280     }
 281   }
 282
 283   RooUnfoldResponse* responseObject = new RooUnfoldResponse(projY, projX,
 284                                                             responseTranspose,
 285                                                             "", "");
 286   RooUnfold*         unfold         = new RooUnfoldBayes(responseObject,
 287                                                          measured, 5);
 288   TH1*               unfolded       = static_cast<TH1*>(unfold->Hreco());
 289   TH1*               triggerBias    = GetTriggerBias(mcHist, esdHist);
 290   DoCorrection(unfolded, TriggerBias, mcHist);
 291
 292   Double_t scale_unf = 1/unfolded->Integral();
 293   unfolded->Scale(scale_unf);
 294   unfolded->Write("Unfolded");
 295
 296   Double_t scale_raw = 1/measured->Integral();
 297   measured->Scale(scale_raw);
 298   measured->Write("Raw");
 299   triggerBias->Write("TriggerBiasCorr");
 300 }
 301
 302
 303
 304 //____________________________________________________________________
 305 TH1* GetTriggerBias(TH1* mcHist, TH1* esdHist)
 306 {
 307   mcHist->Sumw2();
 308   esdHist->Sumw2();
 309
 310   TH1* corr = new TH1D("corr", "corr", 40, -0.5, 39.5);
 311   for (Int_t j=1; j<corr->GetNbinsX(); j++) {
 312     Double_t errSqr = 0.;
 313     Double_t cESD   = esdHist->GetBinContent(j);
 314     Double_t cMS    = mcHist->GetBinContent(j);
 315     Double_t c      = (cMC  == 0 ? 1 :
 316                        cESD == 0 ? 1/cMC : cESD/cMC);
 317     corr->SetBinContent(j, c);
 318     corr->SetBinError(j, c*TMath::Sqrt(errSqr));
 319   }
 320   return corr;
 321 }
 322
 323 //____________________________________________________________________
 324 void DoCorrection(TH1* unfoldBefore, TH1* triggerBias, TH1*)
 325 {
 326   for (Int_t k=1; k<=35; k++) {
 327     Double_t tb = triggerBias->GetBinContent(k);
 328     Double_t ub = unfoldBefore->GetBinContent(k);
 329     if (tb > 1e-5 && ub > 0) {
 330       unfoldBefore->SetBinContent(k, ub / tb);
 331       Double_t eub = UnfoldBefore->GetBinError(k);
 332       Double_t etb = TriggerBias->GetBinError(k);
 333       unfoldBefore->SetBinError(k, TMath::Sqrt(TMath::Power(eub/ub)+
 334                                                TMath::Power(etb/tb))*ub/tb);
 335     }
 336   }
 337 }
 338 //
 339 // EOF
 340 //
 341
 342
 343