PWGLF/SPECTRA/ChargedHadrons/dNdPt/macros/ApplyCorrections_PbPb.C

   1 Double_t ApplyCorrections_PbPb(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
   2 //Double_t ApplyCorrections_PbPb()
   3 {
   4
   5 // tmp setting
   6 //      const char* datafile = "/lustre/alice/train/V006.PbPb/2011-03-15_0009.5917/mergedPeriods/PbPb/2.76ATeV/LHC10h.pass1/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";
   7 //      const char* datatask = "mknichel_dNdPtPbPb_TPCITS_VZERO";
   8 //      const char* corrfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/corrMatr_LHC10h8_TPCITS.root";
   9 //      const char* idstring = "c0";
  10 //      const char* outfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/finalSpectra_LHC10h.pass1_TPCITS.root";
  11 //      const char* gifdir = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/LHC10h.pass1";
  12
  13
  14 Int_t c_first = 1;
  15 Int_t c_last = 11;
  16
  17 TString id = TString(idstring);
  18 if ( id.Contains("c0-5") ) { c_first = c_last = 1; }
  19 if ( id.Contains("c5-10") ) { c_first = c_last = 2; }
  20 if ( id.Contains("c10-20") ) { c_first = c_last = 3; }
  21 if ( id.Contains("c20-30") ) { c_first = c_last = 4; }
  22 if ( id.Contains("c30-40") ) { c_first = c_last = 5; }
  23 if ( id.Contains("c40-50") ) { c_first = c_last = 6; }
  24 if ( id.Contains("c50-60") ) { c_first = c_last = 7; }
  25 if ( id.Contains("c60-70") ) { c_first = c_last = 8; }
  26 if ( id.Contains("c70-80") ) { c_first = c_last = 9; }
  27 if ( id.Contains("c80-90") ) { c_first = c_last = 10; }
  28 if ( id.Contains("c90-100") ) { c_first = c_last = 11; }
  29
  30     // settings vor zVertex cut (event and track level)
  31     Double_t zVert = 10.0;
  32
  33     // setting on eta cut (track level)
  34     Double_t eta = 0.8;
  35
  36     //load required libraries
  37     //load required libraries
  38     gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/ -I$ALICE_ROOT/include -I$ALICE_ROOT/STEER  -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWG0 -I$ALICE_ROOT/PWGPP -I$ALICE_ROOT/PWG2 -I$ALICE_ROOT/PWG3 -I$ALICE_ROOT/PWG3/vertexingHF -I$ALICE_ROOT/PWG4 -I$ALICE_ROOT/CORRFW -I$ALICE_ROOT/TPC -I$ALICE_ROOT/TRD -I$ALICE_ROOT/PWG3/muon -I$ALICE_ROOT/JETAN -I$ALICE_ROOT/ANALYSIS/Tender");
  39
  40   gSystem->Load("libCore");
  41   gSystem->Load("libPhysics");
  42   gSystem->Load("libMinuit");
  43   gSystem->Load("libGui");
  44   gSystem->Load("libXMLParser");
  45
  46   gSystem->Load("libGeom");
  47   gSystem->Load("libVMC");
  48
  49   gSystem->Load("libNet");
  50
  51   gSystem->Load("libSTEERBase");
  52   gSystem->Load("libESD");
  53   gSystem->Load("libCDB");
  54   gSystem->Load("libRAWDatabase");
  55   gSystem->Load("libRAWDatarec");
  56   gSystem->Load("libANALYSIS");
  57
  58
  59
  60     gSystem->Load("libANALYSIS");
  61     gSystem->Load("libANALYSISalice");
  62     gSystem->Load("libTender");
  63     gSystem->Load("libCORRFW");
  64     gSystem->Load("libPWG0base");
  65     gSystem->Load("libPWG0dep");
  66     gSystem->Load("libPWG0selectors");
  67
  68
  69     // make plots nicer
  70     gROOT->SetStyle("Plain");
  71     gStyle->SetPalette(1);
  72
  73     // array for all histograms to be saved
  74     TObjArray* Hists = new TObjArray();
  75
  76     // open file with correction matrices
  77     TFile *fcorr = TFile::Open(corrfile,"READ");
  78     if (!fcorr) return -2;
  79
  80     // load data
  81     TFile* fdata = TFile::Open(datafile,"READ");
  82     if (!fdata) return -1;
  83     TList* ldata = dynamic_cast<TList*>(fdata->Get(datatask));
  84     if (!ldata) return -1;
  85     AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(ldata->FindObject("dNdPtAnalysisPbPb"));
  86     if (!obj) return -1;
  87
  88     //Event statistics
  89     THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events
  90     fRecEventHist2->GetAxis(2)->SetRange(c_first,c_last); // select centrality
  91     TH2D* h2RecEvent2All = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
  92     fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  93     TH2D* h2RecEvent2 = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
  94
  95     Double_t ReconstructedEvents = h2RecEvent2->Integral();
  96     Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
  97
  98     Hists->Add(h2RecEvent2);
  99     Hists->Add(h2RecEvent2All);
 100
 101     printf("=== Number of events from DATA                      %lf ===\n",ReconstructedEvents);
 102     printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);
 103
 104     TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
 105     TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
 106
 107     h1ReconstructedEvents->Fill(0,ReconstructedEvents);
 108     h1ReconstructedEvents->SetEntries(ReconstructedEvents);
 109     h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
 110     h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
 111
 112     Hists->Add(h1ReconstructedEvents);
 113     Hists->Add(h1ReconstructedEventsAll);
 114
 115      // retrieve tracks
 116     THnSparse* fRecTrackHist2 = obj->GetRecTrackHist2(2); //after all cuts (2)
 117     fRecTrackHist2->GetAxis(3)->SetRange(c_first,c_last); // select centrality
 118     fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
 119     fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta
 120
 121     TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
 122     TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
 123     TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
 124
 125     Hists->Add(h1RecTrackHist2_zv);
 126     Hists->Add(h1RecTrackHist2_pt);
 127     Hists->Add(h1RecTrackHist2_eta);
 128
 129     // retrieve correction matrices for tracking efficiency (note different binning!)
 130     TH1D* h1TrackCorr_pt  = (TH1D*)fcorr->Get("h1TrackCorr_pt");
 131     TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
 132
 133     // retrieve correction matrices for secondaries (note different binning!)
 134     TH1D* h1SecCorr_pt  = (TH1D*)fcorr->Get("h1SecCorr_pt");
 135     TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");
 136
 137     // create corrected spectra (as clone of raw data)
 138     TH1D* h1Corrected_pt  = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
 139     TH1D* h1Corrected_eta  = h1RecTrackHist2_eta->Clone("h1Corrected_eta");
 140
 141     // secondaries correction for pt spectrum
 142     for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
 143         Double_t pt = h1Corrected_pt->GetBinCenter(i);
 144         Double_t val = h1Corrected_pt->GetBinContent(i);
 145         Double_t err = h1Corrected_pt->GetBinError(i);
 146         if (pt >= 50) { pt = 49.5; }  // above 50 GeV corr matr have low statistics
 147         Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
 148         Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
 149         Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
 150         Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
 151         Double_t corr    = effCorr*secCorr;
 152         Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated
 153         Double_t cval = val*corr;
 154         Double_t cerr = TMath::Sqrt(val*val*corrErr*corrErr + err*err*corr*corr);
 155         h1Corrected_pt->SetBinContent(i,cval);
 156         h1Corrected_pt->SetBinError(i,cerr);
 157     }
 158
 159
 160     // for eta the correction is simpler because of same binning
 161     h1Corrected_eta->Multiply(h1SecCorr_eta);
 162     h1Corrected_eta->Multiply(h1TrackCorr_eta);
 163
 164     Hists->Add(h1Corrected_pt);
 165     Hists->Add(h1Corrected_eta);
 166
 167     // create final spectra (as clone of corrected data)
 168     TH1D* dNdPt   = h1Corrected_pt->Clone("dNdPt");
 169     TH1D* dNdEta  = h1Corrected_eta->Clone("dNdEta");
 170
 171     // also uncorrected spectra (as clone of raw data)
 172     TH1D* dNdPt_raw   = h1RecTrackHist2_pt->Clone("dNdPt_raw");
 173     TH1D* dNdEta_raw  = h1RecTrackHist2_eta->Clone("dNdEta_raw");
 174
 175     TH1D* dNdPt_nores   = h1Corrected_pt->Clone("dNdPt_nores");
 176
 177
 178 //TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
 179 //TF1 *fcent  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
 180 TFile* fptcorr = TFile::Open("ptcorrPbPb_150511.root");
 181
 182 TF1 * fun = 0;
 183
 184   TF1 *fcent  = new TF1("fcent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,50.);
 185   TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,50.);
 186 Double_t downscale = 1.;
 187             if (c_first != c_last) {
 188              cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
 189              cout << "WARNING: pt resolution correction error!" << endl;
 190              cout << " (works only for single centraliy bins) " << endl;
 191              cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
 192         }
 193         if (c_first == 1) {
 194                 fun = (TF1*) fptcorr->Get("ptcorr_c0");
 195                 downscale = 0.8;
 196                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 197                 cout << "0-5% central pt-resolution correction used!" << endl;
 198                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 199         } else if (c_first == 2) {
 200                 downscale = 0.8;
 201                 fun =(TF1*) fptcorr->Get("ptcorr_c5");
 202                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 203                 cout << "5-10% central pt-resolution correction used!" << endl;
 204                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 205         } else if (c_first == 3) {
 206                 downscale = 0.8;
 207                 fun =(TF1*) fptcorr->Get("ptcorr_c10");
 208                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 209                 cout << "10-20% central pt-resolution correction used!" << endl;
 210                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 211         } else if (c_first == 4) {
 212                 downscale = 0.9;
 213                 fun =(TF1*) fptcorr->Get("ptcorr_c20");
 214                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 215                 cout << "20-30% central pt-resolution correction used!" << endl;
 216                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 217         } else if (c_first == 5) {
 218                 fun =(TF1*) fptcorr->Get("ptcorr_c30");
 219                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 220                 cout << "30-40% central pt-resolution correction used!" << endl;
 221                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 222         } else if (c_first == 6) {
 223                 fun =(TF1*) fptcorr->Get("ptcorr_c40");
 224                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 225                 cout << "40-50% central pt-resolution correction used!" << endl;
 226                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 227         } else if (c_first == 7) {
 228                 fun =(TF1*) fptcorr->Get("ptcorr_c50");
 229                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 230                 cout << "50-60% central pt-resolution correction used!" << endl;
 231                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 232         } else if (c_first == 8) {
 233                 fun =(TF1*) fptcorr->Get("ptcorr_c60");
 234                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 235                 cout << "60-70% central pt-resolution correction used!" << endl;
 236                 cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 237         } else {
 238                 fun =(TF1*) fptcorr->Get("ptcorr_c70");
 239                 cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
 240                 cout << "70-80% central pt-resolution correction used!" << endl;
 241                 cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;
 242         }
 243
 244
 245     // normalization and finalization
 246     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 247
 248    for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
 249         Double_t pt = dNdPt->GetBinCenter(i);
 250         Double_t width = dNdPt->GetBinWidth(i);
 251         Double_t val = dNdPt->GetBinContent(i);
 252         Double_t err = dNdPt->GetBinError(i);
 253         Double_t corrPtResol = 1.0;
 254         corrPtResol = 1.-((1.-fun->Eval(pt))*downscale );
 255         if (pt < 10.) { corrPtResol = 1.0; }
 256         Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 257         Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 258         dNdPt->SetBinContent(i,cval);
 259         dNdPt->SetBinError(i,cerr);
 260     }
 261     // for dndeta again simpler
 262     dNdEta->Scale(1,"width");
 263     dNdEta->Scale(1./ReconstructedEvents);
 264
 265     // normalization and finalization without resolution correction
 266     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 267    for (int i=1; i <= dNdPt_nores->GetNbinsX() ;i++) {
 268         Double_t pt = dNdPt_nores->GetBinCenter(i);
 269         Double_t width = dNdPt_nores->GetBinWidth(i);
 270         Double_t val = dNdPt_nores->GetBinContent(i);
 271         Double_t err = dNdPt_nores->GetBinError(i);
 272         Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 273         Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 274         dNdPt_nores->SetBinContent(i,cval);
 275         dNdPt_nores->SetBinError(i,cerr);
 276     }
 277
 278     // normalization and finalization
 279     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 280    for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
 281         Double_t pt = dNdPt_raw->GetBinCenter(i);
 282         Double_t width = dNdPt_raw->GetBinWidth(i);
 283         Double_t val = dNdPt_raw->GetBinContent(i);
 284         Double_t err = dNdPt_raw->GetBinError(i);
 285         Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 286         Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 287         dNdPt_raw->SetBinContent(i,cval);
 288         dNdPt_raw->SetBinError(i,cerr);
 289     }
 290     // for dndeta again simpler
 291     dNdEta_raw->Scale(1,"width");
 292     dNdEta_raw->Scale(1./ReconstructedEvents);
 293
 294     dNdEta->SetMarkerStyle(21);
 295     dNdPt->SetMarkerStyle(21);
 296     dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T}) (GeV/c)^{-2}");
 297     dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 298
 299     dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T}) (GeV/c)^{-2}");
 300     dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 301
 302     Hists->Add(dNdEta_raw);
 303     Hists->Add(dNdPt_raw);
 304     Hists->Add(dNdEta);
 305     Hists->Add(dNdPt);
 306     Hists->Add(dNdPt_nores);
 307
 308
 309    // plot pictures and save to gifdir
 310     for (i=0; i < Hists->LastIndex(); i++) {
 311         TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
 312         if (gifdir && ctmp) {
 313             TString gif(gifdir);
 314             gif += '/';
 315             gif += ctmp->GetName();
 316             gif += ".gif";
 317             ctmp->SaveAs(gif.Data(),"gif");
 318             delete ctmp;
 319         }
 320     }
 321
 322     // save all correction matrices and control histograms to file
 323     if (!outfile) { return; }
 324     TFile *out = TFile::Open(outfile,"RECREATE");
 325     Hists->Write();
 326     out->Close();
 327
 328     return ReconstructedEvents;
 329
 330 }
 331
 332 //___________________________________________________________________________
 333 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 334 {
 335     TH1* h = dynamic_cast<TH1*>(hobj);
 336     if (!h) return 0;
 337     if (h->GetDimension() > 2) return 0;
 338     h->SetStats(0);
 339     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 340     TString t(label);
 341     if (label) t += "_";
 342     t += h->GetName();
 343     h->SetTitle(t.Data());
 344     TCanvas* c = new TCanvas(t.Data(),t.Data());
 345     if (h->GetDimension() >= 1) {
 346         TString xlabel(h->GetXaxis()->GetTitle());
 347         if (xlabel.Contains("Pt")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 348         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 349     }
 350     if (h->GetDimension() == 2) {
 351         TString ylabel(h->GetYaxis()->GetTitle());
 352         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 353         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 354         h->Draw("COLZ");
 355     }
 356     if (h->GetDimension() == 1) {
 357         h->Draw();
 358     }
 359     return c;
 360
 361 }