PWG0/dNdPt/macros/ApplyCorrections_pp.C

   1 Double_t ApplyCorrections_pp(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
   2 //Double_t ApplyCorrections_pp()
   3 {
   4
   5 // tmp setting
   6 // const char* mcfile = "/d/alice09/mknichel/train/V007.MC_pp/2011-05-05_2347.7147/mergedPeriods/MC_pp/7TeV/LHC10f6a/mknichel_dNdPtpp_TPCITS.root";
   7 // const char* mctask = "mknichel_dNdPtpp_TPCITS";
   8 // const char* idstring = "TPCITS";
   9 // const char* gifdir = "/u/mknichel/alice/dNdPt_pp/temp";
  10 // const char* datafile = "/d/alice09/mknichel/train/V007.MC_pp/2011-05-05_2347.7147/mergedPeriods/MC_pp/7TeV/LHC10f6a/mknichel_dNdPtpp_TPCITS.root";
  11 // const char* outfile = "/u/mknichel/alice/dNdPt_pp/temp/FinalSpectra.root";
  12 // const char* corrfile = "/u/mknichel/alice/dNdPt_pp/temp/corrMatr_LHC10f6a_TPCITS.root";
  13 // const char* gifdir = "/u/mknichel/alice/dNdPt_pp/temp";
  14
  15 TString fname (datafile);
  16
  17 TString id (idstring);
  18 TString taskname = "mknichel_dNdPtpp_" + id;
  19 TString objname = "dNdPtAnalysis_" + id;
  20
  21 // tmp setting
  22  //     const char* datatask = "jotwinow_dNdPtAnalysis_TPCITS";
  23 //      const char* corrfile = "/u/mknichel/alice/dNdPt_pp/2011-04-04_1251/corrMatr_LHC10f6a_TPCITS.root";
  24 //      const char* idstring = "";
  25 //      const char* outfile = "/u/mknichel/alice/dNdPt_pp/2011-04-04_1251/finalSpectra_LHC11a_TPCITS.root";
  26 //      const char* gifdir = "/u/mknichel/alice/dNdPt_pp/2011-04-04_1251/plots";
  27
  28     // settings vor zVertex cut (event and track level)
  29     Double_t zVert = 10.0;
  30
  31     // setting on eta cut (track level)
  32     Double_t eta = 0.8;
  33     /*
  34     //load required libraries
  35     //load required libraries
  36     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");
  37
  38   gSystem->Load("libCore");
  39   gSystem->Load("libPhysics");
  40   gSystem->Load("libMinuit");
  41   gSystem->Load("libGui");
  42   gSystem->Load("libXMLParser");
  43
  44   gSystem->Load("libGeom");
  45   gSystem->Load("libVMC");
  46
  47   gSystem->Load("libNet");
  48
  49   gSystem->Load("libSTEERBase");
  50   gSystem->Load("libESD");
  51   gSystem->Load("libCDB");
  52   gSystem->Load("libRAWDatabase");
  53   gSystem->Load("libRAWDatarec");
  54   gSystem->Load("libANALYSIS");
  55
  56
  57
  58     gSystem->Load("libANALYSIS.so");
  59     gSystem->Load("libOADB.so");
  60     gSystem->Load("libANALYSISalice.so");
  61     gSystem->Load("libTENDER.so");
  62     gSystem->Load("libCORRFW.so");
  63     gSystem->Load("libPWG0base.so");
  64     gSystem->Load("libPWG0dep");
  65     gSystem->Load("libPWG0selectors.so");
  66     */
  67
  68     // make plots nicer
  69     gROOT->SetStyle("Plain");
  70     gStyle->SetPalette(1);
  71
  72     // array for all histograms to be saves
  73     TObjArray* Hists = new TObjArray();
  74
  75     // open file with correction matrices
  76     TFile *fcorr = TFile::Open(corrfile,"READ");
  77     if (!fcorr) return -2;
  78
  79     TH1D* dNdPt_MC      =  (TH1D*) fcorr->Get("dNdPt_MC");
  80     TH1D* dNdEta_MC      = (TH1D*) fcorr->Get("dNdEta_MC");
  81     Hists->Add(dNdPt_MC);
  82     Hists->Add(dNdEta_MC);
  83
  84     // load data
  85     TFile* fdata = TFile::Open(datafile,"READ");
  86     if (!fdata) return -1;
  87     TList* ldata = dynamic_cast<TList*>(fdata->Get(taskname.Data()));
  88     if (!ldata) return -1;
  89     AlidNdPtAnalysis *obj = dynamic_cast<AlidNdPtAnalysis*>(ldata->FindObject(objname.Data()));
  90     if (!obj) return -1;
  91
  92     //Event statistics
  93     obj->GetRecEventHist2()->GetAxis(0)->SetRange();
  94     THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events
  95     TH1D* h1RecEventHist2_zv = (TH1D*) fRecEventHist2->Projection(0)->Clone("h1RecEventHist2_zv");  // zvertex distribution
  96     TH2D* h2RecEvent2All = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
  97     fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
  98     TH2D* h2RecEvent2 = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
  99     TH2D* h2RecEvent2Corrected = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2Corrected"); //will be corrected
 100
 101     THnSparse* fEventCount = obj->GetEventCount();
 102     Hists->Add(fEventCount);
 103     Double_t TriggeredEventsNoVertex = fEventCount->Projection(0)->GetBinContent(2) -  fEventCount->Projection(1)->GetBinContent(2); // all triggered events without rec. vertex
 104     Double_t AllTriggeredEvents = fEventCount->Projection(0)->GetBinContent(2); // all
 105
 106     Double_t ReconstructedEvents = h2RecEvent2->Integral();
 107     Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
 108     Double_t TriggerEffInel = 0.864; //numeber from michele floris for 7TeV (I use it also for 2.76)
 109     TriggerEffInel = 0.9379; //mc value from mc
 110     if ( fname.Contains("900GeV") ) { TriggerEffInel = 0.9156; }
 111     if ( fname.Contains("2.76TeV") ) { TriggerEffInel = 0.883; }
 112     if ( fname.Contains("7TeV") ) { TriggerEffInel = 0.8524; }
 113
 114     cout << "Using Trigger_to_Inel efficiecy: " << TriggerEffInel <<endl;
 115
 116
 117     Double_t ReconstructedEventsFraction = (ReconstructedEventsAll-ReconstructedEvents) / ReconstructedEventsAll;
 118     Double_t SelectedEventsFraction = ReconstructedEvents / ReconstructedEventsAll;
 119     Double_t InelasticEventsSimple = AllTriggeredEvents/TriggerEffInel*SelectedEventsFraction;
 120
 121     Hists->Add(h2RecEvent2);
 122     Hists->Add(h2RecEvent2All);
 123     Hists->Add(h2RecEvent2Corrected);
 124
 125     printf("=== Number of events from DATA                      %lf ===\n",ReconstructedEvents);
 126     printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);
 127     printf("=== Number of events from DATA (all triggered)      %lf ===\n",AllTriggeredEvents);
 128
 129     TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
 130     TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
 131
 132     h1ReconstructedEvents->Fill(0,ReconstructedEvents);
 133     h1ReconstructedEvents->SetEntries(ReconstructedEvents);
 134     h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
 135     h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
 136
 137     Hists->Add(h1ReconstructedEvents);
 138     Hists->Add(h1ReconstructedEventsAll);
 139
 140      // retrieve corrections (event level)
 141      TH2D* h2EventTriggerEffAll   = (TH2D*)fcorr->Get("h2EventTriggerEffAll");
 142      TH2D* h2EventTriggerCorrAll  = (TH2D*)fcorr->Get("h2EventTriggerCorrAll");
 143      TH2D* h2EventTriggerEff      = (TH2D*)fcorr->Get("h2EventTriggerEff");
 144      TH2D* h2EventTriggerCorr     = (TH2D*)fcorr->Get("h2EventTriggerCorr");
 145      TH2D* h2EventRecEffAll       = (TH2D*)fcorr->Get("h2EventRecEffAll");
 146      TH2D* h2EventRecCorrAll      = (TH2D*)fcorr->Get("h2EventRecCorrAll");
 147      TH2D* h2EventRecEff          = (TH2D*)fcorr->Get("h2EventRecEff");
 148      TH2D* h2EventRecCorr         = (TH2D*)fcorr->Get("h2EventRecCorr");
 149      TH2D* h2EventEffAll          = (TH2D*)fcorr->Get("h2EventEffAll");
 150      TH2D* h2EventCorrAll         = (TH2D*)fcorr->Get("h2EventCorrAll");
 151      TH2D* h2EventEff             = (TH2D*)fcorr->Get("h2EventEff");
 152      TH2D* h2EventCorr            = (TH2D*)fcorr->Get("h2EventCorr");
 153      TH1D* h1TriggerEff_bin0_zv   = (TH1D*)fcorr->Get("h1TriggerEff_bin0_zv");
 154      TH1D* h1TriggerCorr_bin0_zv  = (TH1D*)fcorr->Get("h1TriggerCorr_bin0_zv");
 155      TH1D* h1Ratio_zv             = (TH1D*)fcorr->Get("h1Ratio_zv");
 156      TH1D* corr_shape_trig0_zv    = (TH1D*)fcorr->Get("corr_shape_trig0_zv");
 157      TH1D* corr_shape_notrig0_zv  = (TH1D*)fcorr->Get("corr_shape_notrig0_zv");
 158
 159
 160     Double_t corrTrackMatch = 1.0;
 161    //tracking efficiency is 2% to low in MC LHC10e13 900GeV as compared to data
 162    if ( fname.Contains("900GeV") ) {
 163        corrTrackMatch = 1.02;
 164        cout << "900 GeV: correct tracking matching efficiency " <<endl;
 165    }
 166
 167
 168     // correct bin0
 169     //h1RecEventHist2_zv->Scale(1./h1RecEventHist2_zv->Integral());
 170     /*
 171     TH1D* h1Bin0_zv =  (TH1D*) h1RecEventHist2_zv->Clone("h1Bin0_zv");
 172     h1Bin0_zv->Multiply(h1Ratio_zv);
 173     h1Bin0_zv->Scale(1./h1Bin0_zv->Integral());
 174     h1Bin0_zv->Scale(TriggeredEventsNoVertex);
 175     //h1Bin0_zv->Multiply(h1TriggerCorr_bin0_zv);
 176     h1Bin0_zv->GetXaxis()->SetRangeUser(-zVert, zVert-0.01); //zVertex
 177     Double_t bin0EventsCorrected =  h1Bin0_zv->Integral();
 178     */
 179     // this is what is used for normalization
 180     Double_t InelasticEventsAll =  AllTriggeredEvents/TriggerEffInel;
 181     Double_t Bin0EventsAll = TriggeredEventsNoVertex;
 182     Double_t UntriggeredEventsAll = InelasticEventsAll - AllTriggeredEvents;
 183     cout << "cross check: INEL Events 1 " << InelasticEventsAll <<endl;
 184     cout << "cross check: INEL Events 2 " << (ReconstructedEventsAll+Bin0EventsAll+UntriggeredEventsAll) <<endl;
 185
 186     // correct bin0
 187     TH1D* h1Bin0_zv =  (TH1D*) h1RecEventHist2_zv->Clone("h1Bin0_zv");
 188     h1Bin0_zv->Multiply(corr_shape_trig0_zv);     //correct for shape
 189     h1Bin0_zv->Scale(1./h1Bin0_zv->Integral());    //normalize
 190     h1Bin0_zv->Scale(Bin0EventsAll);
 191     Double_t Bin0Events = h1Bin0_zv->Integral(5,8);
 192
 193     //correct untriggered
 194     TH1D* h1NoTrig_zv =  (TH1D*) h1RecEventHist2_zv->Clone("h1NoTrig_zv");
 195     h1NoTrig_zv->Multiply(corr_shape_notrig0_zv);     //correct for shape
 196     h1NoTrig_zv->Scale(1./h1NoTrig_zv->Integral());    //normalize
 197     h1NoTrig_zv->Scale(UntriggeredEventsAll);
 198     Double_t UntriggeredEvents = h1NoTrig_zv->Integral(5,8);
 199
 200
 201     Double_t InelasticEvents = ReconstructedEvents + Bin0Events + UntriggeredEvents;
 202
 203     /*
 204     // my way (old)
 205     TH1D* h1Bin0_zv =  (TH1D*) h1RecEventHist2_zv->Clone("h1Bin0_zv");
 206     h1Bin0_zv->Multiply(h1Ratio_zv);
 207     h1Bin0_zv->Multiply(h1TriggerEff_bin0_zv);
 208     h1Bin0_zv->Scale(1./h1Bin0_zv->Integral());
 209     h1Bin0_zv->Scale(TriggeredEventsNoVertex);
 210     h1Bin0_zv->Multiply(h1TriggerCorr_bin0_zv);
 211     Double_t bin0EventsCorrected = h1Bin0_zv->Integral(5,8);
 212      */
 213
 214     /*
 215     //jaceks way (?)
 216     h1RecEventHist2_zv->Scale(1./h1RecEventHist2_zv->Integral());
 217     h1RecEventHist2_zv->Scale(TriggeredEventsNoVertex);
 218     Double_t bin0EventsCorrected = h1RecEventHist2_zv->Integral(5,8);
 219 */
 220      // retrieve tracks
 221     THnSparse* hSRecTrack; // = obj->GetRecTrackHist()->Clone("hsRecTrack"); //thnsparse to be corrected
 222     THnSparse* fRecTrackHist2 = obj->GetRecTrackHist(); //after all cuts (2)
 223     fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
 224     fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta
 225     Int_t dims[4] = {0,1,2,3};
 226     hSRecTrack = (THnSparse*) (fRecTrackHist2->Projection(4,dims,"e")->Clone("hSRecTrack"));
 227     hSRecTrackAllMult = (THnSparse*) (fRecTrackHist2->Projection(3,dims,"e")->Clone("hSRecTrackAllMult"));
 228
 229     TH3D* h3RecTrackHist2 = fRecTrackHist2->Projection(0,1,2)->Clone("h3RecTrackHist2");
 230
 231     TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
 232     TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
 233     TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
 234
 235     Hists->Add(h1RecTrackHist2_zv);
 236     Hists->Add(h1RecTrackHist2_pt);
 237     Hists->Add(h1RecTrackHist2_eta);
 238     Hists->Add(h3RecTrackHist2);
 239
 240     // generate corrections for events
 241     TH1D* h1MCGeneratedEvents     = (TH1D*)fcorr->Get("h1MCGeneratedEvents");
 242     TH1D* h1MCReconstructedEvents = (TH1D*)fcorr->Get("h1MCReconstructedEvents");
 243     TH1D* h1MCTriggeredEvents0mult = (TH1D*)fcorr->Get("h1MCTriggeredEvents0mult");
 244     TH1D* h1MCTriggeredEventsAll0mult = (TH1D*)fcorr->Get("h1MCTriggeredEventsAll0mult");
 245     TH1D* h1MCTriggeredEventsAll = (TH1D*)fcorr->Get("h1MCTriggeredEventsAll");
 246     TH1D* h1MCTriggeredEvents = (TH1D*)fcorr->Get("h1MCTriggeredEvents");
 247
 248
 249     Double_t MCGeneratedEvents = h1MCGeneratedEvents->GetEntries();
 250     Double_t MCReconstructedEvents = h1MCReconstructedEvents->GetEntries();
 251     Double_t CorrEvent = MCGeneratedEvents / MCReconstructedEvents;
 252
 253     Double_t MCTriggeredEvents = h1MCTriggeredEvents->GetEntries();
 254     Double_t MCTriggeredEventsAll = h1MCTriggeredEventsAll->GetEntries();
 255     Double_t MCTriggeredEvents0mult = h1MCTriggeredEvents0mult->GetEntries();
 256     Double_t MCTriggeredEventsAll0mult = h1MCTriggeredEventsAll0mult->GetEntries();
 257     Double_t CorrVtxEvent0mult    = MCTriggeredEvents / (MCTriggeredEvents-MCTriggeredEvents0mult); //this is used
 258     Double_t CorrVtxEventAll0mult = MCTriggeredEventsAll / (MCTriggeredEventsAll-MCTriggeredEventsAll0mult);
 259
 260     // correct for trigger/vertex inefficiencies
 261     for (int xbin=1; xbin <= h2RecEvent2Corrected->GetNbinsX(); xbin++) {
 262         for (int ybin=1; ybin <= h2RecEvent2Corrected->GetNbinsY(); ybin++) {
 263             Double_t x = h2RecEvent2Corrected->GetXaxis()->GetBinCenter(xbin);
 264             Double_t y = h2RecEvent2Corrected->GetYaxis()->GetBinCenter(ybin);
 265             Int_t bin  = h2EventCorr->FindBin(x,y);
 266             Double_t corr    = h2EventCorr->GetBinContent(bin);
 267             Double_t correrr = h2EventCorr->GetBinError(bin);
 268             if (corr < 0.01) { corr=1.; correrr=0.;} //bin empty in correction matrix
 269             Double_t val     = h2RecEvent2Corrected->GetBinContent(xbin,ybin);
 270             Double_t err     = h2RecEvent2Corrected->GetBinError(xbin,ybin);
 271             h2RecEvent2Corrected->SetBinContent(xbin,ybin,corr*val);
 272             h2RecEvent2Corrected->SetBinError(xbin,ybin,TMath::Sqrt(val*val*correrr*correrr + err*err*corr*corr));
 273
 274         }
 275     }
 276     //Double_t ReconstructedEventsCorrected = h2RecEvent2Corrected->Integral() * CorrVtxEvent0mult;
 277     //Double_t ReconstructedEventsCorrected = h2RecEvent2Corrected->Integral() + bin0EventsCorrected;
 278     printf("=== Number of events from DATA (final correction)   %lf ===\n",InelasticEvents);
 279     printf("=== Number of events from DATA (corr, no bin0)      %lf ===\n",h2RecEvent2Corrected->Integral());
 280     printf("=== Number of events from DATA (simple corr. 7TeV)  %lf ===\n",InelasticEventsSimple);
 281     printf("=== Number of events from DATA (simple correction)  %lf ===\n",ReconstructedEvents * CorrEvent);
 282     printf("=== Number of bin0events from DATA                  %lf ===\n",TriggeredEventsNoVertex);
 283     //printf("=== Number of bin0events from DATA (corrected)      %lf ===\n",bin0EventsCorrected);
 284
 285     // retrieve 3D correction matrices for tracking efficiency and secondaries (note different binning!)
 286     TH3D* h3TrackEff   = (TH3D*)fcorr->Get("h3TrackEff");
 287     TH3D* h3TrackCorr  = (TH3D*)fcorr->Get("h3TrackCorr");
 288     TH3D* h3SecCont    = (TH3D*)fcorr->Get("h3SecCont");
 289     TH3D* h3SecCorr    = (TH3D*)fcorr->Get("h3SecCorr");
 290
 291     // retrieve 3D thnsparse correction matrices for tracking efficiency and secondaries (note different binning!) --- this is used!!!!
 292     THnSparse* hSTrackEff = (THnSparse*) fcorr->Get("hSTrackEff");
 293     THnSparse* hSTrackCorr = (THnSparse*) fcorr->Get("hSTrackCorr");
 294     THnSparse* hSSecCont = (THnSparse*) fcorr->Get("hSSecCont");
 295     THnSparse* hSSecCorr = (THnSparse*) fcorr->Get("hSSecCorr");
 296
 297
 298     // retrieve correction matrices for tracking efficiency (note different binning!)
 299     TH1D* h1TrackCorr_pt  = (TH1D*)fcorr->Get("h1TrackCorr_pt");
 300     TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
 301
 302     // retrieve correction matrices for secondaries (note different binning!)
 303     TH1D* h1SecCorr_pt  = (TH1D*)fcorr->Get("h1SecCorr_pt");
 304     TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");
 305
 306     // create corrected spectra (as clone of raw data)
 307     TH1D* h1Corrected_pt  = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
 308     TH1D* h1Corrected_eta  = h1RecTrackHist2_eta->Clone("h1Corrected_eta");
 309
 310     // secondaries correction for pt spectrum
 311     for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
 312         Double_t pt = h1Corrected_pt->GetBinCenter(i);
 313         Double_t val = h1Corrected_pt->GetBinContent(i);
 314         Double_t err = h1Corrected_pt->GetBinError(i);
 315         if (pt >= 50) { pt = 49.5; }  // above 50 GeV corr matr have low statistics
 316         Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
 317         Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
 318         Double_t cval = val*secCorr;
 319         Double_t cerr = TMath::Sqrt(val*val*secCorrErr*secCorrErr + err*err*secCorr*secCorr);
 320         h1Corrected_pt->SetBinContent(i,cval);
 321         h1Corrected_pt->SetBinError(i,cerr);
 322     }
 323
 324     // tracking efficiency correction pt spectrum
 325     for (int i=1; i <= h1Corrected_pt->GetNbinsX() ;i++) {
 326         Double_t pt = h1Corrected_pt->GetBinCenter(i);
 327         Double_t val = h1Corrected_pt->GetBinContent(i);
 328         Double_t err = h1Corrected_pt->GetBinError(i);
 329         if (pt >= 50) { pt = 49.5; } // above 50 GeV corr matr have low statistics
 330         Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
 331         Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
 332
 333         Double_t cval = corrTrackMatch * val * effCorr;
 334         Double_t cerr = corrTrackMatch * TMath::Sqrt(val*val*effCorrErr*effCorrErr + err*err*effCorr*effCorr);
 335         h1Corrected_pt->SetBinContent(i,cval);
 336         h1Corrected_pt->SetBinError(i,cerr);
 337     }
 338
 339
 340
 341    // efficiency and contamination correction for thnsparse
 342    for (Long64_t j = 0; j < hSRecTrack->GetNbins(); j++) {
 343        Int_t tc[4];
 344        Double_t tval    = hSRecTrack->GetBinContent(j,tc);
 345        Double_t terr    = hSRecTrack->GetBinError(j);
 346        Double_t tzv     = hSRecTrack->GetAxis(0)->GetBinCenter(tc[0]);
 347        Double_t tpt     = hSRecTrack->GetAxis(1)->GetBinCenter(tc[1]);
 348        Double_t teta    = hSRecTrack->GetAxis(2)->GetBinCenter(tc[2]);
 349        Double_t tmultMB = hSRecTrack->GetAxis(3)->GetBinCenter(tc[3]);
 350        if (tzv >= 10.0) { tzv = 9.9; }
 351        if (teta >= 0.8) {teta = 0.79;}
 352        if (tzv < -10.0) { tzv = -10.; }
 353        if (teta < -0.8) { teta = -0.8; }
 354        if (tpt >= 50.) { tpt = 49.5; } // above 50 GeV corr matr have low statistics
 355        if (tpt < 0.15) { tpt = 0.175; } // also below 0.15
 356        Double_t xvals[3];
 357        xvals[0] = tzv;
 358        xvals[1] = tpt;
 359        xvals[2] = teta;
 360        Double_t effCorr    = hSTrackCorr->GetBinContent(hSTrackCorr->GetBin(xvals,kFALSE));
 361        Double_t effCorrErr = hSTrackCorr->GetBinError(hSTrackCorr->GetBin(xvals,kFALSE));
 362        Double_t secCorr    = hSSecCont->GetBinContent(hSSecCont->GetBin(xvals,kFALSE));
 363        Double_t secCorrErr = hSSecCont->GetBinError(hSSecCont->GetBin(xvals,kFALSE));
 364 //        Double_t effCorr    = h3TrackCorr->GetBinContent(h3TrackCorr->FindBin(tzv,tpt,teta));
 365 //        Double_t effCorrErr = h3TrackCorr->GetBinError(h3TrackCorr->FindBin(tzv,tpt,teta));
 366 //        Double_t secCorr    = 1. - h3SecCont->GetBinContent(h3SecCont->FindBin(tzv,tpt,teta));
 367 //        Double_t secCorrErr = h3SecCont->GetBinError(h3SecCont->FindBin(tzv,tpt,teta));
 368        if (effCorr < 1.) { cout << "bin empty, use efficiency 1!" << endl; effCorr=1.; }
 369        if (secCorr < 0.001) { cout << "bin empty, use contamination 0!" << endl; }
 370        if (effCorrErr < 1e-9) { cout << "eff error empty!" << endl; }
 371        if (secCorrErr < 1e-9) { cout << "cont error empty!" << endl;}
 372        Double_t corr    = corrTrackMatch * effCorr*(1.-secCorr);
 373        //Double_t corrErr = corrTrackMatch * TMath::Sqrt(effCorr*secCorrErr + secCorr*effCorrErr);
 374        Double_t ctval = tval*corr;
 375        Double_t cterr = terr*corr; // errors are not correlated
 376        //Double_t cterr = TMath::Sqrt(tval*tval*corrErr*corrErr + terr*terr*corr*corr); // errors are not correlated
 377        hSRecTrack->SetBinContent(j,ctval);
 378        hSRecTrack->SetBinError(j,cterr);
 379     }
 380
 381     TH1D* effCorrPt = hSRecTrack->Projection(1)->Clone("effCorrPt");
 382     TH1D* effCorrPtErr2 = hSRecTrack->Projection(1)->Clone("effCorrPtErr2");
 383     TH1D* effCorrPtNorm = hSRecTrack->Projection(1)->Clone("effCorrPtNorm");
 384     effCorrPt->Reset();
 385     effCorrPtErr2->Reset();
 386 //     effCorrPtNorm->Reset();
 387     /*
 388 for (Long64_t j = 0; j < hSRecTrackAllMult->GetNbins(); j++) {
 389        Int_t tc[3];
 390        Double_t tval    = hSRecTrackAllMult->GetBinContent(j,tc);
 391        Double_t tzv     = hSRecTrackAllMult->GetAxis(0)->GetBinCenter(tc[0]);
 392        Double_t tpt     = hSRecTrackAllMult->GetAxis(1)->GetBinCenter(tc[1]);
 393        Double_t teta    = hSRecTrackAllMult->GetAxis(2)->GetBinCenter(tc[2]);
 394        if (tpt >= 50.) { tpt = 49.5; } // above 50 GeV corr matr have low statistics
 395        Double_t effCorr    = h3TrackCorr->GetBinContent(h3TrackCorr->FindBin(tzv,tpt,teta));
 396        Double_t effCorrErr = h3TrackCorr->GetBinError(h3TrackCorr->FindBin(tzv,tpt,teta));
 397        Double_t secCorr    = 1. - h3SecCont->GetBinContent(h3SecCont->FindBin(tzv,tpt,teta));
 398        Double_t secCorrErr = h3SecCont->GetBinError(h3SecCont->FindBin(tzv,tpt,teta));
 399        if (effCorr < 0.1) { cout << "bin empty, use efficiency 1!" << endl; effCorr=1.; effCorrErr=0.; }
 400        if (secCorr < 0.1) { cout << "bin empty, use contamination 0!" << endl; secCorr=1.; secCorrErr=0.; }
 401        Double_t corr    = effCorr*secCorr;
 402        Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated
 403        effCorrPt->Fill(tpt,tval*corr);
 404        effCorrPtErr2->Fill(tpt,tval*tval*corrErr*corrErr);
 405        effCorrPtNorm->Fill(tpt,tval);
 406     }
 407
 408     effCorrPt->Divide(effCorrPtNorm);
 409     for (Int_t i = 4; i <= effCorrPtErr2->GetNbinsX(); i++) { effCorrPtErr2->SetBinContent(i,TMath::Sqrt(effCorrPtErr2->GetBinContent(i)));  }
 410     effCorrPtErr2->Divide(effCorrPtNorm);
 411     cout << effCorrPt->GetBinContent(4) << endl;
 412     cout << effCorrPtErr2->GetBinContent(4) << endl;
 413     cout << effCorrPtNorm->GetBinContent(4) << endl;
 414
 415
 416
 417     TH1D* RecTrackPtCorrected = hSRecTrack->Projection(1)->Clone("RecTrackPtCorrected");
 418     RecTrackPtCorrected->Reset();
 419
 420
 421     for (Int_t i = 4; i <= RecTrackPtCorrected->GetNbinsX(); i++) {
 422        if (0 == effCorrPtErr2->GetBinContent(i)) continue;
 423        Double_t val     = effCorrPt->GetBinContent(i);
 424        cout << (val*val/effCorrPtNorm->GetBinContent(i)) << " " << effCorrPtErr2->GetBinContent(i) << endl;
 425        Double_t err2  = (val*val/effCorrPtNorm->GetBinContent(i)) + effCorrPtErr2->GetBinContent(i);
 426        RecTrackPtCorrected->SetBinContent(i,val);
 427        RecTrackPtCorrected->SetBinError(i,TMath::Sqrt(err2));
 428     }
 429     */
 430
 431 for (Long64_t j = 0; j < hSRecTrackAllMult->GetNbins(); j++) {
 432        Int_t tc[3];
 433        Double_t tval    = hSRecTrackAllMult->GetBinContent(j,tc);
 434        Double_t terr    = TMath::Sqrt(tval); //hSRecTrackAllMult->GetBinError(j);
 435        Double_t tzv     = hSRecTrackAllMult->GetAxis(0)->GetBinCenter(tc[0]);
 436        Double_t tpt     = hSRecTrackAllMult->GetAxis(1)->GetBinCenter(tc[1]);
 437        Double_t teta    = hSRecTrackAllMult->GetAxis(2)->GetBinCenter(tc[2]);
 438        Double_t tmultMB = hSRecTrack->GetAxis(3)->GetBinCenter(tc[3]);
 439        if (tzv >= 10.0) { tzv = 9.9; }
 440        if (teta >= 0.8) {teta = 0.79;}
 441        if (tzv < -10.0) { tzv = -10.; }
 442        if (teta < -0.8) { teta = -0.8; }
 443        if (tpt >= 50.) { tpt = 49.5; } // above 50 GeV corr matr have low statistics
 444        if (tpt < 0.15) { tpt = 0.175; } // also below 0.15
 445       Double_t xvals[3];
 446        xvals[0] = tzv;
 447        xvals[1] = tpt;
 448        xvals[2] = teta;
 449        Double_t effCorr    = hSTrackCorr->GetBinContent(hSTrackCorr->GetBin(xvals,kFALSE));
 450        Double_t effCorrErr = hSTrackCorr->GetBinError(hSTrackCorr->GetBin(xvals,kFALSE));
 451        Double_t secCont    = hSSecCont->GetBinContent(hSSecCont->GetBin(xvals,kFALSE));
 452        Double_t secContErr = hSSecCont->GetBinError(hSSecCont->GetBin(xvals,kFALSE));
 453 //        Double_t effCorr    = h3TrackCorr->GetBinContent(h3TrackCorr->FindBin(tzv,tpt,teta));
 454 //        Double_t effCorrErr = h3TrackCorr->GetBinError(h3TrackCorr->FindBin(tzv,tpt,teta));
 455 //        Double_t secCorr    = 1. - h3SecCont->GetBinContent(h3SecCont->FindBin(tzv,tpt,teta));
 456 //        Double_t secCorrErr = h3SecCont->GetBinError(h3SecCont->FindBin(tzv,tpt,teta));
 457        if (effCorr < 1.) { cout << "bin empty, use efficiency 1!" << endl; effCorr=1.; }
 458        if (secCont < 0.001) { cout << "bin empty, use contamination 0!" << endl; }
 459        if (effCorrErr < 1e-9) { cout << "eff error empty!" << endl; }
 460        if (secContErr < 1e-9) { cout << "cont error empty!" << endl;}
 461        Double_t secCorr = (1.-secCont);
 462        Double_t corr    = corrTrackMatch * (effCorr*secCorr);
 463        Double_t corrErr = corrTrackMatch *  TMath::Sqrt(effCorr*secCorrErr + secCorr*effCorrErr);
 464        Double_t ctval = tval*corr;
 465        Double_t cterr = terr*corr; // errors are not correlated
 466        //Double_t cterr = TMath::Sqrt(tval*tval*corrErr*corrErr + terr*terr*corr*corr); // errors are not correlated
 467        hSRecTrackAllMult->SetBinContent(j,ctval);
 468        hSRecTrackAllMult->SetBinError(j,cterr);
 469     }
 470
 471     // create final spectrum in multiplicity bins
 472     TH2D* dNdPtMult = hSRecTrack->Projection(3,1)->Clone("dNdPtMult");
 473
 474     // create final spectra (as clone of corrected data)
 475     TH1D* dNdPt = hSRecTrackAllMult->Projection(1)->Clone("dNdPt");
 476 /*
 477     // errors on pt spectrum from corrections
 478     for (int i=1; i <= dNdPt->GetNbinsX() ; i++) {
 479         Double_t pt = dNdPt->GetBinCenter(i);
 480         Double_t val = dNdPt->GetBinContent(i);
 481         Double_t err = dNdPt->GetBinError(i);
 482         if (pt >= 50) { pt = 49.5; }  // above 50 GeV corr matr have low statistics
 483         Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
 484         Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
 485         Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
 486         Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
 487         Double_t corr    = effCorr*secCorr;
 488         Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated
 489         Double_t cval = val*secCorr;
 490         Double_t cerr = TMath::Sqrt(val*val*corrErr*corrErr + err*err*corr*corr);
 491         dNdPt->SetBinError(i,cerr);
 492     }
 493
 494     // errors on pt spectrum from corrections
 495     for (int i=1; i <= dNdPtMult->GetNbinsX() ; i++) {
 496     for (int k=1; k <= dNdPtMult->GetNbinsY() ; k++) {
 497         Double_t pt = dNdPt->GetBinCenter(i);
 498         Double_t val = dNdPtMult->GetBinContent(i,k);
 499         Double_t err = dNdPtMult->GetBinError(i,k);
 500         if (pt >= 50) { pt = 49.5; }  // above 50 GeV corr matr have low statistics
 501         Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
 502         Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
 503         Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
 504         Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
 505         Double_t corr    = effCorr*secCorr;
 506         Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated
 507         Double_t cval = val*secCorr;
 508         Double_t cerr = TMath::Sqrt(val*val*corrErr*corrErr + err*err*corr*corr);
 509         dNdPtMult->SetBinError(i,k,cerr);
 510     }
 511     }
 512 */
 513
 514     // for eta the correction is simpler because of same binning
 515     h1Corrected_eta->Multiply(h1SecCorr_eta);
 516     h1Corrected_eta->Multiply(h1TrackCorr_eta);
 517
 518     Hists->Add(h1Corrected_pt);
 519     Hists->Add(h1Corrected_eta);
 520
 521
 522     TH1D* dNdPt_simple = (TH1D*) h1Corrected_pt->Clone("dNdPt_simple");
 523     TH1D* dNdPt_simple_nores  = (TH1D*) h1Corrected_pt->Clone("dNdPt_simple_nores");
 524
 525
 526     TH1D* dNdEta  = hSRecTrackAllMult->Projection(2)->Clone("dNdEta");
 527     TH1D* dNdEta_simple  = h1Corrected_eta->Clone("dNdEta_simple");
 528
 529     // also uncorrected spectra (as clone of raw data)
 530     TH1D* dNdPt_raw   = (TH1D*) h1RecTrackHist2_pt->Clone("dNdPt_raw");
 531     TH1D* dNdEta_raw  = (TH1D*) h1RecTrackHist2_eta->Clone("dNdEta_raw");
 532
 533     // also uncorrected spectra (as clone of raw data)
 534     TH1D* dNdPt_nores          = (TH1D*) dNdPt->Clone("dNdPt_nores");
 535     //TH1D* dNdEta_nores  = h1Corrected_eta->Clone("dNdEta_nores");
 536
 537
 538
 539 //TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,50.);
 540 //TF1 *fcent  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,50.);
 541 TFile* fptcorr = TFile::Open("ptcorr_140511.root");
 542
 543 TF1 *fun = 0;
 544 TString id = TString(idstring);
 545 if ( fname.Contains("900GeV") ) {
 546   TF1 *fun  = (TF1*) fptcorr->Get("ptcorr_900");
 547   cout << "++++++++++++++++++++++++++++++++++++++" << endl;
 548   cout << "900GeV pt-resolution correction used!" << endl;
 549   cout << "++++++++++++++++++++++++++++++++++++++" << endl;
 550 } elseif ( fname.Contains("2.76TeV") ) {
 551   TF1 *fun  = (TF1*) fptcorr->Get("ptcorr_2760");
 552   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 553   cout << "2.76TeV pt-resolution correction used!" << endl;
 554   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 555 } elseif ( fname.Contains("7TeV") ) {
 556   TF1 *fun  = (TF1*) fptcorr->Get("ptcorr_7000");
 557   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 558   cout << "7TeV pt-resolution correction used!" << endl;
 559   cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
 560 } else {
 561   fun = TF1("none","1.",10.,50.);
 562   cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
 563   cout << "ERROR! NO pt-resolution correction used!" << endl;
 564   cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
 565 }
 566
 567     // normalization and finalization
 568     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 569    for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
 570         Double_t pt = dNdPt->GetBinCenter(i);
 571         Double_t width = dNdPt->GetBinWidth(i);
 572         Double_t val = dNdPt->GetBinContent(i);
 573         Double_t err = dNdPt->GetBinError(i);
 574         Double_t corrPtResol = 1.0;
 575         corrPtResol = fun->Eval(pt);
 576         if (pt < 10.) { corrPtResol = 1.0; }
 577         if (corrPtResol > 1.0 ) { corrPtResol = 1.0; }
 578         Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 579         Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 580 //         Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * CorrEvent* pt);
 581 //         Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * CorrEvent* pt);
 582         dNdPt->SetBinContent(i,cval);
 583         dNdPt->SetBinError(i,cerr);
 584     }
 585
 586     // normalization and finalization without resolution correction
 587     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 588    for (int i=1; i <= dNdPt_nores->GetNbinsX() ;i++) {
 589         Double_t pt = dNdPt_nores->GetBinCenter(i);
 590         Double_t width = dNdPt_nores->GetBinWidth(i);
 591         Double_t val = dNdPt_nores->GetBinContent(i);
 592         Double_t err = dNdPt_nores->GetBinError(i);
 593         Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 594         Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 595         dNdPt_nores->SetBinContent(i,cval);
 596         dNdPt_nores->SetBinError(i,cerr);
 597     }
 598
 599     // normalization and finalization without resolution correction
 600     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 601    for (int i=1; i <= dNdPt_simple_nores->GetNbinsX() ;i++) {
 602         Double_t pt = dNdPt_simple_nores->GetBinCenter(i);
 603         Double_t width = dNdPt_simple_nores->GetBinWidth(i);
 604         Double_t val = dNdPt_simple_nores->GetBinContent(i);
 605         Double_t err = dNdPt_simple_nores->GetBinError(i);
 606         Double_t cval = (val)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 607         Double_t cerr = (err)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 608         dNdPt_simple_nores->SetBinContent(i,cval);
 609         dNdPt_simple_nores->SetBinError(i,cerr);
 610     }
 611
 612     // for dndeta again simpler
 613     dNdEta->Scale(1,"width");
 614     dNdEta->Scale(1./InelasticEvents);
 615     dNdEta_simple->Scale(1,"width");
 616     dNdEta_simple->Scale(1./InelasticEvents);
 617
 618
 619     // normalization and finalization
 620     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 621    for (int i=1; i <= dNdPt_simple->GetNbinsX() ;i++) {
 622         Double_t pt = dNdPt_simple->GetBinCenter(i);
 623         Double_t width = dNdPt_simple->GetBinWidth(i);
 624         Double_t val = dNdPt_simple->GetBinContent(i);
 625         Double_t err = dNdPt_simple->GetBinError(i);
 626         Double_t corrPtResol = 1.0;
 627         corrPtResol = fun->Eval(pt);
 628         if (pt < 10.) { corrPtResol = 1.0; }
 629         if (corrPtResol > 1.0 ) { corrPtResol = 1.0; }
 630         Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 631         Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 632         dNdPt_simple->SetBinContent(i,cval);
 633         dNdPt_simple->SetBinError(i,cerr);
 634     }
 635     /*
 636     TH1D* dNdPt2 = (TH1D*) dNdPt->Clone("dNdPt2");
 637     TH1D* dNdPt2_simple = (TH1D*) dNdPt_simple->Clone("dNdPt2_simple");
 638     dNdPt2->Scale(InelasticEvents / (ReconstructedEvents * CorrEvent));
 639     dNdPt2_simple->Scale(InelasticEvents / (ReconstructedEvents * CorrEvent));
 640     Hists->Add(dNdPt2);
 641     Hists->Add(dNdPt2_simple);
 642     */
 643
 644     // normalization and finalization
 645     // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
 646    for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
 647         Double_t pt = dNdPt_raw->GetBinCenter(i);
 648         Double_t width = dNdPt_raw->GetBinWidth(i);
 649         Double_t val = dNdPt_raw->GetBinContent(i);
 650         Double_t err = dNdPt_raw->GetBinError(i);
 651         Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 652         Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * InelasticEvents * pt);
 653         dNdPt_raw->SetBinContent(i,cval);
 654         dNdPt_raw->SetBinError(i,cerr);
 655     }
 656     // for dndeta again simpler
 657     dNdEta_raw->Scale(1,"width");
 658     dNdEta_raw->Scale(1./InelasticEvents);
 659
 660     dNdEta->SetMarkerStyle(21);
 661     dNdPt->SetMarkerStyle(21);
 662     dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 663     dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 664
 665     dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 666     dNdPt_nores->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
 667     dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
 668
 669     Hists->Add(dNdEta_raw);
 670     Hists->Add(dNdPt_raw);
 671     Hists->Add(dNdEta);
 672     Hists->Add(dNdEta_simple);
 673     Hists->Add(dNdPt);
 674     Hists->Add(dNdPt_nores);
 675     Hists->Add(dNdPt_simple);
 676     Hists->Add(dNdPt_simple_nores);
 677     Hists->Add(hSRecTrack);
 678     Hists->Add(hSRecTrackAllMult);
 679
 680
 681    // plot pictures and save to gifdir
 682     for (i=0; i < Hists->LastIndex(); i++) {
 683         TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
 684         if (gifdir && ctmp) {
 685             TString gif(gifdir);
 686             gif += '/';
 687             gif += ctmp->GetName();
 688             gif += ".gif";
 689             ctmp->SaveAs(gif.Data(),"gif");
 690             delete ctmp;
 691         }
 692     }
 693
 694     // save all correction matrices and control histograms to file
 695     if (!outfile) { return; }
 696     TFile *out = TFile::Open(outfile,"RECREATE");
 697     Hists->Write();
 698     out->Close();
 699
 700     return ReconstructedEvents;
 701
 702 }
 703
 704 //___________________________________________________________________________
 705 TCanvas* PlotHist(TObject* hobj, const char* label=0)
 706 {
 707     TH1* h = dynamic_cast<TH1*>(hobj);
 708     if (!h) return 0;
 709     if (h->GetDimension() > 2) return 0;
 710     h->SetStats(0);
 711     if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); }
 712     TString t(label);
 713     if (label) t += "_";
 714     t += h->GetName();
 715     h->SetTitle(t.Data());
 716     TCanvas* c = new TCanvas(t.Data(),t.Data());
 717     if (h->GetDimension() >= 1) {
 718         TString xlabel(h->GetXaxis()->GetTitle());
 719         if (xlabel.Contains("Pt")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 720         if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }
 721     }
 722     if (h->GetDimension() == 2) {
 723         TString ylabel(h->GetYaxis()->GetTitle());
 724         if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 725         if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
 726         h->Draw("COLZ");
 727     }
 728     if (h->GetDimension() == 1) {
 729         h->Draw();
 730     }
 731     return c;
 732
 733 }