TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / 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 }