PWGUD/dNdPt/macros/ApplyCorrections.C

   1 Double_t ApplyCorrections(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
   2 //void ApplyCorrections()
   3 {
   4
   5 // tmp setting
   6 //     const char* datafile = "/lustre/alice/train/V005.PbPb/2010-11-28_0258.4230/mergedPeriods/data/PbPb/LHC10h.pass1/mknichel_dNdPtPbPb_gt_v0_c0.root";
   7 //     const char* datatask = "mknichel_dNdPtPbPb_gt_v0_c0";
   8 //     const char* corrfile = "/u/mknichel/alice/dNdPt/2010-11-28/corrMatr_gt_v0_c0.root";
   9 //     const char* idstring = "gt_v0_c0";
  10 //     const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-28/finalSpectra_v0_c0.root";
  11 //     const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-28";
  12
  13     // settings vor zVertex cut (event and track level)
  14     Double_t zVert = 10.0;
  15
  16     // setting on eta cut (track level)
  17     Double_t eta = 0.8;
  18
  19     //load required libraries
  20     //load required libraries
  21     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");
  22
  23   gSystem->Load("libCore");
  24   gSystem->Load("libPhysics");
  25   gSystem->Load("libMinuit");
  26   gSystem->Load("libGui");
  27   gSystem->Load("libXMLParser");
  28
  29   gSystem->Load("libGeom");
  30   gSystem->Load("libVMC");
  31
  32   gSystem->Load("libNet");
  33
  34   gSystem->Load("libSTEERBase");
  35   gSystem->Load("libESD");
  36   gSystem->Load("libCDB");
  37   gSystem->Load("libRAWDatabase");
  38   gSystem->Load("libRAWDatarec");
  39   gSystem->Load("libANALYSIS");
  40
  41
  42
  43     gSystem->Load("libANALYSIS.so");
  44     gSystem->Load("libANALYSISalice.so");
  45     gSystem->Load("libTENDER.so");
  46     gSystem->Load("libCORRFW.so");
  47     gSystem->Load("libPWG0base.so");
  48     gSystem->Load("libPWG0dep");
  49     gSystem->Load("libPWG0selectors.so");
  50
  51
  52     // make plots nicer
  53     gROOT->SetStyle("Plain");
  54     gStyle->SetPalette(1);
  55
  56     // array for all histograms to be saves
  57     TObjArray* Hists = new TObjArray();
  58
  59     // open file with correction matrices
  60     TFile *fcorr = TFile::Open(corrfile,"READ");
  61     if (!fcorr) return -2;
  62
  63     // load data
  64     TFile* fdata = TFile::Open(datafile,"READ");
  65     if (!fdata) return -1;
  66     TList* ldata = dynamic_cast<TList*>(fdata->Get(datatask));
  67     if (!ldata) return -1;
  68     AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(ldata->FindObject("dNdPtAnalysisPbPb"));
  69     if (!obj) return -1;
  70
  71     //Event statistics
  72     THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events
  73     TH2D* h2RecEvent2All = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
  74     fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  75     TH2D* h2RecEvent2 = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
  76
  77     Double_t ReconstructedEvents = h2RecEvent2->Integral();
  78     Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
  79
  80     Hists->Add(h2RecEvent2);
  81     Hists->Add(h2RecEvent2All);
  82
  83     printf("=== Number of events from DATA                      %lf ===\n",ReconstructedEvents);
  84     printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);
  85
  86     TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
  87     TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
  88
  89     h1ReconstructedEvents->Fill(0,ReconstructedEvents);
  90     h1ReconstructedEvents->SetEntries(ReconstructedEvents);
  91     h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
  92     h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
  93
  94     Hists->Add(h1ReconstructedEvents);
  95     Hists->Add(h1ReconstructedEventsAll);
  96
  97      // retrieve tracks
  98     THnSparse* fRecTrackHist2 = obj->GetRecTrackHist2(2); //after all cuts (2)
  99     fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
 100     fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta
 101
 102     TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
 103     TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
 104     TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
 105
 106     Hists->Add(h1RecTrackHist2_zv);
 107     Hists->Add(h1RecTrackHist2_pt);
 108     Hists->Add(h1RecTrackHist2_eta);
 109
 110     // retrieve correction matrices for tracking efficiency (note different binning!)
 111     TH1D* h1TrackCorr_pt  = (TH1D*)fcorr->Get("h1TrackCorr_pt");
 112     TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
 113
 114     // retrieve correction matrices for secondaries (note different binning!)
 115     TH1D* h1SecCorr_pt  = (TH1D*)fcorr->Get("h1SecCorr_pt");
 116     TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");
 117
 118     // create corrected spectra (as clone of raw data)
 119     TH1D* h1Corrected_pt  = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
 120     TH1D* h1Corrected_eta  = h1RecTrackHist2_eta->Clone("h1Corrected_eta");
 121
 122     // secondaries correction for pt spectrum
 123     for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
 124         Double_t pt = h1Corrected_pt->GetBinCenter(i);
 125         Double_t val = h1Corrected_pt->GetBinContent(i);
 126         Double_t err = h1Corrected_pt->GetBinError(i);
 127         if (pt >= 19) { pt = 19; }  // above 20 GeV corr matr have low statistics
 128         Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
 129         Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
 130         Double_t cval = val*secCorr;
 131         Double_t cerr = TMath::Sqrt(val*val*secCorrErr*secCorrErr + err*err*secCorr*secCorr);
 132         h1Corrected_pt->SetBinContent(i,cval);
 133         h1Corrected_pt->SetBinError(i,cerr);
 134     }
 135
 136     // tracking efficiency correction pt spectrum
 137     for (int i=1; i <= h1Corrected_pt->GetNbinsX() ;i++) {
 138         Double_t pt = h1Corrected_pt->GetBinCenter(i);
 139         Double_t val = h1Corrected_pt->GetBinContent(i);
 140         Double_t err = h1Corrected_pt->GetBinError(i);
 141         if (pt >= 19) { pt = 19; } // above 20 GeV corr matr have low statistics
 142         Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
 143         Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
 144
 145         Double_t cval = val*effCorr;
 146         Double_t cerr = TMath::Sqrt(val*val*effCorrErr*effCorrErr + err*err*effCorr*effCorr);
 147         h1Corrected_pt->SetBinContent(i,cval);
 148         h1Corrected_pt->SetBinError(i,cerr);
 149     }
 150
 151     // for eta the correction is simpler because of same binning
 152     h1Corrected_eta->Multiply(h1SecCorr_eta);
 153     h1Corrected_eta->Multiply(h1TrackCorr_eta);
 154
 155     Hists->Add(h1Corrected_pt);
 156     Hists->Add(h1Corrected_eta);
 157
 158     // create final spectra (as clone of corrected data)
 159     TH1D* dNdPt   = h1Corrected_pt->Clone("dNdPt");
 160     TH1D* dNdEta  = h1Corrected_eta->Clone("dNdEta");
 161
 162     // also uncorrected spectra (as clone of raw data)
 163     TH1D* dNdPt_raw   = h1RecTrackHist2_pt->Clone("dNdPt_raw");
 164     TH1D* dNdEta_raw  = h1RecTrackHist2_eta->Clone("dNdEta_raw");
 165
 166
 167 //TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
 168 //TF1 *fcent  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
 169 TString id = TString(idstring);
 170 if ( id.Contains("c0") || id.Contains("c5") ) {
 171   TF1 *fun  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
 172   cout << "++++++++++++++++++++++++++++++++++++++" << endl;
 173   cout << "CENTRAL pt-resolution correction used!" << endl;
 174   cout << "++++++++++++++++++++++++++++++++++++++" << endl;
 175 } else {
 176 TF1 *fun  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
 177   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 178   cout << "PERIPHERAL pt-resolution correction used!" << endl;
 179   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 180 }
 181
 182     // normalization and finalization
 183     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 184    for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
 185         Double_t pt = dNdPt->GetBinCenter(i);
 186         Double_t width = dNdPt->GetBinWidth(i);
 187         Double_t val = dNdPt->GetBinContent(i);
 188         Double_t err = dNdPt->GetBinError(i);
 189         Double_t corrPtResol = 1.0;
 190         corrPtResol = fun->Eval(pt);
 191         if (pt < 5.) { corrPtResol = 1.0; }
 192         Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 193         Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 194         dNdPt->SetBinContent(i,cval);
 195         dNdPt->SetBinError(i,cerr);
 196     }
 197     // for dndeta again simpler
 198     dNdEta->Scale(1,"width");
 199     dNdEta->Scale(1./ReconstructedEvents);
 200
 201     // normalization and finalization
 202     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 203    for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
 204         Double_t pt = dNdPt_raw->GetBinCenter(i);
 205         Double_t width = dNdPt_raw->GetBinWidth(i);
 206         Double_t val = dNdPt_raw->GetBinContent(i);
 207         Double_t err = dNdPt_raw->GetBinError(i);
 208         Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 209         Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
 210         dNdPt_raw->SetBinContent(i,cval);
 211         dNdPt_raw->SetBinError(i,cerr);
 212     }
 213     // for dndeta again simpler
 214     dNdEta_raw->Scale(1,"width");
 215     dNdEta_raw->Scale(1./ReconstructedEvents);
 216
 217     dNdEta->SetMarkerStyle(21);
 218     dNdPt->SetMarkerStyle(21);
 219     dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 220     dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 221
 222     dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 223     dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 224
 225     Hists->Add(dNdEta_raw);
 226     Hists->Add(dNdPt_raw);
 227     Hists->Add(dNdEta);
 228     Hists->Add(dNdPt);
 229
 230
 231    // plot pictures and save to gifdir
 232     for (i=0; i < Hists->LastIndex(); i++) {
 233         TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
 234         if (gifdir && ctmp) {
 235             TString gif(gifdir);
 236             gif += '/';
 237             gif += ctmp->GetName();
 238             gif += ".gif";
 239             ctmp->SaveAs(gif.Data(),"gif");
 240             delete ctmp;
 241         }
 242     }
 243
 244     // save all correction matrices and control histograms to file
 245     if (!outfile) { return; }
 246     TFile *out = TFile::Open(outfile,"RECREATE");
 247     Hists->Write();
 248     out->Close();
 249
 250     return ReconstructedEvents;
 251
 252 }
 253
 254 //___________________________________________________________________________
 255 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 256 {
 257     TH1* h = dynamic_cast<TH1*>(hobj);
 258     if (!h) return 0;
 259     if (h->GetDimension() > 2) return 0;
 260     h->SetStats(0);
 261     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 262     TString t(label);
 263     if (label) t += "_";
 264     t += h->GetName();
 265     h->SetTitle(t.Data());
 266     TCanvas* c = new TCanvas(t.Data(),t.Data());
 267     if (h->GetDimension() >= 1) {
 268         TString xlabel(h->GetXaxis()->GetTitle());
 269         if (xlabel.Contains("Pt")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }
 270         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }
 271     }
 272     if (h->GetDimension() == 2) {
 273         TString ylabel(h->GetYaxis()->GetTitle());
 274         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
 275         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
 276         h->Draw("COLZ");
 277     }
 278     if (h->GetDimension() == 1) {
 279         h->Draw();
 280     }
 281     return c;
 282
 283 }