[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / ChargedHadrons / dNdPt / macros / ApplyCorrections.C

Double_t ApplyCorrections(const char* datafile, const char* datatask, const char* corrfile, const char* idstring ,const char* outfile, const char* gifdir = 0)
//void ApplyCorrections()
{

// tmp setting
//     const char* datafile = "/lustre/alice/train/V005.PbPb/2010-11-28_0258.4230/mergedPeriods/data/PbPb/LHC10h.pass1/mknichel_dNdPtPbPb_gt_v0_c0.root";
//     const char* datatask = "mknichel_dNdPtPbPb_gt_v0_c0";
//     const char* corrfile = "/u/mknichel/alice/dNdPt/2010-11-28/corrMatr_gt_v0_c0.root";
//     const char* idstring = "gt_v0_c0";
//     const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-28/finalSpectra_v0_c0.root";
//     const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-28";
    
    // settings vor zVertex cut (event and track level)
    Double_t zVert = 10.0;
    
    // setting on eta cut (track level)
    Double_t eta = 0.8;
    
    //load required libraries
    //load required libraries    
    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");
    
  gSystem->Load("libCore");
  gSystem->Load("libPhysics");
  gSystem->Load("libMinuit");
  gSystem->Load("libGui");
  gSystem->Load("libXMLParser");

  gSystem->Load("libGeom");
  gSystem->Load("libVMC");

  gSystem->Load("libNet");

  gSystem->Load("libSTEERBase");
  gSystem->Load("libESD");
  gSystem->Load("libCDB");
  gSystem->Load("libRAWDatabase");
  gSystem->Load("libRAWDatarec");
  gSystem->Load("libANALYSIS");

    
    
    gSystem->Load("libANALYSIS");
    gSystem->Load("libANALYSISalice");
    gSystem->Load("libTender");
    gSystem->Load("libCORRFW");
    gSystem->Load("libPWG0base");
    gSystem->Load("libPWG0dep"); 
    gSystem->Load("libPWG0selectors");
    

    // make plots nicer
    gROOT->SetStyle("Plain");
    gStyle->SetPalette(1);
    
    // array for all histograms to be saves
    TObjArray* Hists = new TObjArray();
    
    // open file with correction matrices
    TFile *fcorr = TFile::Open(corrfile,"READ");
    if (!fcorr) return -2;
    
    // load data
    TFile* fdata = TFile::Open(datafile,"READ");
    if (!fdata) return -1;
    TList* ldata = dynamic_cast<TList*>(fdata->Get(datatask));
    if (!ldata) return -1;
    AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(ldata->FindObject("dNdPtAnalysisPbPb"));
    if (!obj) return -1;
    
    //Event statistics
    THnSparse *fRecEventHist2 = obj->GetRecEventHist2(); //reconstructed events	
    TH2D* h2RecEvent2All = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
    fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
    TH2D* h2RecEvent2 = fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2");
    
    Double_t ReconstructedEvents = h2RecEvent2->Integral();
    Double_t ReconstructedEventsAll = h2RecEvent2All->Integral();
    
    Hists->Add(h2RecEvent2);
    Hists->Add(h2RecEvent2All);
  
    printf("=== Number of events from DATA                      %lf ===\n",ReconstructedEvents);
    printf("=== Number of events from DATA (before zVertex cut) %lf ===\n",ReconstructedEventsAll);        
  
    TH1D* h1ReconstructedEvents = new TH1D("h1ReconstructedEvents","h1ReconstructedEvents",1,0,1);
    TH1D* h1ReconstructedEventsAll = new TH1D("h1ReconstructedEventsAll","h1ReconstructedEventsAll",1,0,1);
    
    h1ReconstructedEvents->Fill(0,ReconstructedEvents);
    h1ReconstructedEvents->SetEntries(ReconstructedEvents);
    h1ReconstructedEventsAll->Fill(0,ReconstructedEventsAll);
    h1ReconstructedEventsAll->SetEntries(ReconstructedEventsAll);
        
    Hists->Add(h1ReconstructedEvents);
    Hists->Add(h1ReconstructedEventsAll);

     // retrieve tracks
    THnSparse* fRecTrackHist2 = obj->GetRecTrackHist2(2); //after all cuts (2)
    fRecTrackHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01); //zVertex
    fRecTrackHist2->GetAxis(2)->SetRangeUser(-eta, eta-0.01); // eta   
   
    TH1D* h1RecTrackHist2_zv = fRecTrackHist2->Projection(0)->Clone("h1RecTrackHist2_zv");
    TH1D* h1RecTrackHist2_pt = fRecTrackHist2->Projection(1)->Clone("h1RecTrackHist2_pt");
    TH1D* h1RecTrackHist2_eta = fRecTrackHist2->Projection(2)->Clone("h1RecTrackHist2_eta");
   
    Hists->Add(h1RecTrackHist2_zv);
    Hists->Add(h1RecTrackHist2_pt);
    Hists->Add(h1RecTrackHist2_eta);

    // retrieve correction matrices for tracking efficiency (note different binning!)
    TH1D* h1TrackCorr_pt  = (TH1D*)fcorr->Get("h1TrackCorr_pt");  
    TH1D* h1TrackCorr_eta = (TH1D*)fcorr->Get("h1TrackCorr_eta");
        
    // retrieve correction matrices for secondaries (note different binning!)
    TH1D* h1SecCorr_pt  = (TH1D*)fcorr->Get("h1SecCorr_pt");  
    TH1D* h1SecCorr_eta = (TH1D*)fcorr->Get("h1SecCorr_eta");

    // create corrected spectra (as clone of raw data)
    TH1D* h1Corrected_pt  = h1RecTrackHist2_pt->Clone("h1Corrected_pt");
    TH1D* h1Corrected_eta  = h1RecTrackHist2_eta->Clone("h1Corrected_eta");

    // secondaries correction for pt spectrum
    for (int i=1; i <= h1Corrected_pt->GetNbinsX() ; i++) {
        Double_t pt = h1Corrected_pt->GetBinCenter(i);
        Double_t val = h1Corrected_pt->GetBinContent(i);
        Double_t err = h1Corrected_pt->GetBinError(i);
        if (pt >= 19) { pt = 19; }  // above 20 GeV corr matr have low statistics
        Double_t secCorr    = h1SecCorr_pt->GetBinContent(h1SecCorr_pt->FindBin(pt));
        Double_t secCorrErr = h1SecCorr_pt->GetBinError(h1SecCorr_pt->FindBin(pt));
        Double_t cval = val*secCorr;
        Double_t cerr = TMath::Sqrt(val*val*secCorrErr*secCorrErr + err*err*secCorr*secCorr);
        h1Corrected_pt->SetBinContent(i,cval);
        h1Corrected_pt->SetBinError(i,cerr);
    }

    // tracking efficiency correction pt spectrum
    for (int i=1; i <= h1Corrected_pt->GetNbinsX() ;i++) {
        Double_t pt = h1Corrected_pt->GetBinCenter(i);
        Double_t val = h1Corrected_pt->GetBinContent(i);
        Double_t err = h1Corrected_pt->GetBinError(i);
        if (pt >= 19) { pt = 19; } // above 20 GeV corr matr have low statistics
        Double_t effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
        Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));
        
        Double_t cval = val*effCorr;
        Double_t cerr = TMath::Sqrt(val*val*effCorrErr*effCorrErr + err*err*effCorr*effCorr);
        h1Corrected_pt->SetBinContent(i,cval);
        h1Corrected_pt->SetBinError(i,cerr);
    }
    
    // for eta the correction is simpler because of same binning
    h1Corrected_eta->Multiply(h1SecCorr_eta);
    h1Corrected_eta->Multiply(h1TrackCorr_eta);
    
    Hists->Add(h1Corrected_pt);
    Hists->Add(h1Corrected_eta);

    // create final spectra (as clone of corrected data)
    TH1D* dNdPt   = h1Corrected_pt->Clone("dNdPt");
    TH1D* dNdEta  = h1Corrected_eta->Clone("dNdEta");
       
    // also uncorrected spectra (as clone of raw data)
    TH1D* dNdPt_raw   = h1RecTrackHist2_pt->Clone("dNdPt_raw");
    TH1D* dNdEta_raw  = h1RecTrackHist2_eta->Clone("dNdEta_raw");
       

//TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
//TF1 *fcent  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);
TString id = TString(idstring);
if ( id.Contains("c0") || id.Contains("c5") ) {
  TF1 *fun  = new TF1("cent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,40.);    
  cout << "++++++++++++++++++++++++++++++++++++++" << endl;
  cout << "CENTRAL pt-resolution correction used!" << endl;
  cout << "++++++++++++++++++++++++++++++++++++++" << endl;
} else {
TF1 *fun  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,40.);
  cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
  cout << "PERIPHERAL pt-resolution correction used!" << endl;
  cout << "+++++++++++++++++++++++++++++++++++++++++" << endl;
}

    // normalization and finalization
    // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
   for (int i=1; i <= dNdPt->GetNbinsX() ;i++) {
        Double_t pt = dNdPt->GetBinCenter(i);
        Double_t width = dNdPt->GetBinWidth(i);
        Double_t val = dNdPt->GetBinContent(i);
        Double_t err = dNdPt->GetBinError(i);
        Double_t corrPtResol = 1.0;
        corrPtResol = fun->Eval(pt);
        if (pt < 5.) { corrPtResol = 1.0; }
        Double_t cval = (val * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
        Double_t cerr = (err * corrPtResol)/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
        dNdPt->SetBinContent(i,cval);
        dNdPt->SetBinError(i,cerr);
    }
    // for dndeta again simpler
    dNdEta->Scale(1,"width");
    dNdEta->Scale(1./ReconstructedEvents);
    
    // normalization and finalization
    // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
   for (int i=1; i <= dNdPt_raw->GetNbinsX() ;i++) {
        Double_t pt = dNdPt_raw->GetBinCenter(i);
        Double_t width = dNdPt_raw->GetBinWidth(i);
        Double_t val = dNdPt_raw->GetBinContent(i);
        Double_t err = dNdPt_raw->GetBinError(i);
        Double_t cval = val/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
        Double_t cerr = err/(width * 2.0 * TMath::Pi() * 1.6 * ReconstructedEvents * pt);
        dNdPt_raw->SetBinContent(i,cval);
        dNdPt_raw->SetBinError(i,cerr);
    }
    // for dndeta again simpler
    dNdEta_raw->Scale(1,"width");
    dNdEta_raw->Scale(1./ReconstructedEvents);
     
    dNdEta->SetMarkerStyle(21);
    dNdPt->SetMarkerStyle(21);
    dNdPt->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
    dNdEta->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
    
    dNdPt_raw->SetTitle("; p_{T} (GeV/c) ; 1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{ch})/(d#eta dp_{T})^{-2}");
    dNdEta_raw->SetTitle("; #eta ; 1/N_{evt} (d^{2}N_{ch})/(d#eta)");
    
    Hists->Add(dNdEta_raw);
    Hists->Add(dNdPt_raw);
    Hists->Add(dNdEta);
    Hists->Add(dNdPt);
    
    
   // plot pictures and save to gifdir
    for (i=0; i < Hists->LastIndex(); i++) {    
        TCanvas* ctmp = PlotHist(Hists->At(i),idstring);
        if (gifdir && ctmp) {
            TString gif(gifdir);
            gif += '/';
            gif += ctmp->GetName();
            gif += ".gif";
            ctmp->SaveAs(gif.Data(),"gif");     
            delete ctmp;
        }
    }  

    // save all correction matrices and control histograms to file
    if (!outfile) { return; }
    TFile *out = TFile::Open(outfile,"RECREATE");
    Hists->Write();
    out->Close();    

    return ReconstructedEvents;

}

//___________________________________________________________________________
TCanvas* PlotHist(TObject* hobj, const char* label=0)
{
    TH1* h = dynamic_cast<TH1*>(hobj);
    if (!h) return 0;
    if (h->GetDimension() > 2) return 0;
    h->SetStats(0);
    if ( TString(h->GetName()).Contains("Events")) { h->SetStats(1); } 
    TString t(label);
    if (label) t += "_";
    t += h->GetName();
    h->SetTitle(t.Data());
    TCanvas* c = new TCanvas(t.Data(),t.Data());
    if (h->GetDimension() >= 1) {
        TString xlabel(h->GetXaxis()->GetTitle());
        if (xlabel.Contains("Pt")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }
        if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 50.); }        
    }
    if (h->GetDimension() == 2) {  
        TString ylabel(h->GetYaxis()->GetTitle());
        if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
        if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 50.); }
        h->Draw("COLZ");
    }        
    if (h->GetDimension() == 1) {
        h->Draw();
    }
    return c;

}
Commit	Line	Data
7c96ea29	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
2bfe5463	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");
7c96ea29	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
4070f709	43	gSystem->Load("libANALYSIS");
	44	gSystem->Load("libANALYSISalice");
	45	gSystem->Load("libTender");
	46	gSystem->Load("libCORRFW");
	47	gSystem->Load("libPWG0base");
7c96ea29	48	gSystem->Load("libPWG0dep");
4070f709	49	gSystem->Load("libPWG0selectors");
7c96ea29	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(valvalsecCorrErrsecCorrErr + errerrsecCorrsecCorr);
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(valvaleffCorrErreffCorrErr + errerreffCorreffCorr);
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.000608425x-6.7038e-05xx",5.,40.);
168	//TF1 fcent = new TF1("cent","1.01074e+00-1.98127e-03x-1.19903e-04xx",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-03x-1.19903e-04xx",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.000608425x-6.7038e-05xx",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	}