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

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

// tmp setting
//      const char* datafile = "/lustre/alice/train/V006.PbPb/2011-03-15_0009.5917/mergedPeriods/PbPb/2.76ATeV/LHC10h.pass1/mknichel_dNdPtPbPb_TPCITS_VZERO1.root";     
//      const char* datatask = "mknichel_dNdPtPbPb_TPCITS_VZERO";
//      const char* corrfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/corrMatr_LHC10h8_TPCITS.root";     
//      const char* idstring = "c0";
//      const char* outfile = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/finalSpectra_LHC10h.pass1_TPCITS.root";     
//      const char* gifdir = "/u/mknichel/alice/dNdPt_PbPb/2011-03-15/LHC10h.pass1";
 
    
Int_t c_first = 1;
Int_t c_last = 11;

TString id = TString(idstring);
if ( id.Contains("c0-5") ) { c_first = c_last = 1; }
if ( id.Contains("c5-10") ) { c_first = c_last = 2; }
if ( id.Contains("c10-20") ) { c_first = c_last = 3; }
if ( id.Contains("c20-30") ) { c_first = c_last = 4; }
if ( id.Contains("c30-40") ) { c_first = c_last = 5; }
if ( id.Contains("c40-50") ) { c_first = c_last = 6; }
if ( id.Contains("c50-60") ) { c_first = c_last = 7; }
if ( id.Contains("c60-70") ) { c_first = c_last = 8; }
if ( id.Contains("c70-80") ) { c_first = c_last = 9; }
if ( id.Contains("c80-90") ) { c_first = c_last = 10; }
if ( id.Contains("c90-100") ) { c_first = c_last = 11; }
    
    // 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 saved
    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	
    fRecEventHist2->GetAxis(2)->SetRange(c_first,c_last); // select centrality    
    TH2D* h2RecEvent2All = (TH2D*) fRecEventHist2->Projection(0,1)->Clone("h2RecEvent2All");
    fRecEventHist2->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
    TH2D* h2RecEvent2 = (TH2D*) 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(3)->SetRange(c_first,c_last); // select centrality    
    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 >= 50) { pt = 49.5; }  // above 50 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 effCorr    = h1TrackCorr_pt->GetBinContent(h1TrackCorr_pt->FindBin(pt));
        Double_t effCorrErr = h1TrackCorr_pt->GetBinError(h1TrackCorr_pt->FindBin(pt));        
        Double_t corr    = effCorr*secCorr;
        Double_t corrErr = effCorr*secCorrErr + secCorr*effCorrErr; // errors are correlated        
        Double_t cval = val*corr;
        Double_t cerr = TMath::Sqrt(val*val*corrErr*corrErr + err*err*corr*corr);        
        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");
    
    TH1D* dNdPt_nores   = h1Corrected_pt->Clone("dNdPt_nores");    
       

//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.);
TFile* fptcorr = TFile::Open("ptcorrPbPb_150511.root");

TF1 * fun = 0;

  TF1 *fcent  = new TF1("fcent","1.01074e+00-1.98127e-03*x-1.19903e-04*x*x",5.,50.);    
  TF1 *fperi  = new TF1("fperi","1.00343-0.000608425*x-6.7038e-05*x*x",5.,50.);
Double_t downscale = 1.;
            if (c_first != c_last) {
             cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
             cout << "WARNING: pt resolution correction error!" << endl;
             cout << " (works only for single centraliy bins) " << endl;
             cout << "++++++++++++++++++++++++++++++++++++++++" << endl;
        }
        if (c_first == 1) {
		fun = (TF1*) fptcorr->Get("ptcorr_c0");
		downscale = 0.8;
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "0-5% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	
        } else if (c_first == 2) { 
        	downscale = 0.8;
		fun =(TF1*) fptcorr->Get("ptcorr_c5");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "5-10% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	        
        } else if (c_first == 3) { 
        	downscale = 0.8;
		fun =(TF1*) fptcorr->Get("ptcorr_c10");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "10-20% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	        
        } else if (c_first == 4) { 
        	downscale = 0.9;
		fun =(TF1*) fptcorr->Get("ptcorr_c20");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "20-30% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	        
        } else if (c_first == 5) { 
		fun =(TF1*) fptcorr->Get("ptcorr_c30");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "30-40% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	        
        } else if (c_first == 6) { 
		fun =(TF1*) fptcorr->Get("ptcorr_c40");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "40-50% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	        
        } else if (c_first == 7) { 
		fun =(TF1*) fptcorr->Get("ptcorr_c50");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "50-60% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	                
	} else if (c_first == 8) { 
		fun =(TF1*) fptcorr->Get("ptcorr_c60");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;
 		cout << "60-70% central pt-resolution correction used!" << endl;
 		cout << "+++++++++++++++++++++++++++++++++++++++++++" << endl;   	                   
        } else {
        	fun =(TF1*) fptcorr->Get("ptcorr_c70");
  		cout << "+++++++++++++++++++++++++++++++++++++++++++++" << endl;  		
 		cout << "70-80% central 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 = 1.-((1.-fun->Eval(pt))*downscale );
        if (pt < 10.) { 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 without resolution correction
    // 1/N_evt 1/(2 pi pt) 1/width 1/etarange
   for (int i=1; i <= dNdPt_nores->GetNbinsX() ;i++) {
        Double_t pt = dNdPt_nores->GetBinCenter(i);
        Double_t width = dNdPt_nores->GetBinWidth(i);
        Double_t val = dNdPt_nores->GetBinContent(i);
        Double_t err = dNdPt_nores->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_nores->SetBinContent(i,cval);
        dNdPt_nores->SetBinError(i,cerr);
    }    
    
    // 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}) (GeV/c)^{-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}) (GeV/c)^{-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);
    Hists->Add(dNdPt_nores);
    
    
   // 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 , 100.); }
        if (xlabel.Contains("p_{T}")) { c->SetLogx();  c->SetLogy();  h->GetXaxis()->SetRangeUser(0.1 , 100.); }        
    }
    if (h->GetDimension() == 2) {  
        TString ylabel(h->GetYaxis()->GetTitle());
        if (ylabel.Contains("Pt")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
        if (ylabel.Contains("p_{T}")) { c->SetLogy(); h->GetYaxis()->SetRangeUser(0.1 , 100.); }
        h->Draw("COLZ");
    }        
    if (h->GetDimension() == 1) {
        h->Draw();
    }
    return c;

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