--- /dev/null
+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/PWG1 -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.so");
+ gSystem->Load("libANALYSISalice.so");
+ gSystem->Load("libTENDER.so");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libPWG0base.so");
+ gSystem->Load("libPWG0dep");
+ gSystem->Load("libPWG0selectors.so");
+
+
+ // 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;
+
+}
--- /dev/null
+Double_t GenerateCorrMatr(const char* mcfile, const char* mctask, const char* idstring ,const char* outfile, const char* gifdir = 0)
+//Double_t GenerateCorrMatr()
+{
+
+//tmp setting
+// const char* mcfile = "/lustre/alice/train/V005.PbPb_MC/2010-11-29_0108.4251/mergedPeriods/MC/PbPb/LHC10h8/mknichel_dNdPtPbPb_gt_v0_c0_syst4.root";
+// const char* mctask = "mknichel_dNdPtPbPb_gt_v0_c0_syst4";
+// const char* idstring = "gt_v0_c0_syst4";
+// const char* outfile = "/u/mknichel/alice/dNdPt/2010-11-29_0318/corrMatr_gt_v0_c0_syst4.root";
+// const char* gifdir = "/u/mknichel/alice/dNdPt/2010-11-29";
+
+ // settings vor zVertex cut (event and track level)
+ Double_t zVert = 10.0;
+
+ // setting on eta cut (track level)
+ Double_t eta = 0.8;
+
+ // strangeness scaling factor (for secondaries from strange decays)
+ Double_t sscale = 1.4;
+
+ //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/PWG1 -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.so");
+ gSystem->Load("libANALYSISalice.so");
+ gSystem->Load("libTENDER.so");
+ gSystem->Load("libCORRFW.so");
+ gSystem->Load("libPWG0base.so");
+ gSystem->Load("libPWG0dep");
+ gSystem->Load("libPWG0selectors.so");
+
+
+
+
+
+ // make plots nicer
+ gROOT->SetStyle("Plain");
+ gStyle->SetPalette(1);
+
+ // array for all correction matrices
+ TObjArray* CorrMatr = new TObjArray();
+
+ // array for all control histograms
+ TObjArray* ContHist = new TObjArray();
+
+
+ // load mc information
+ TFile* fmc = TFile::Open(mcfile,"READ");
+ if (!fmc) return -1;
+ TList* lmc = dynamic_cast<TList*>(fmc->Get(mctask));
+ if (!lmc) return -1;
+ AlidNdPtAnalysisPbPb *obj = dynamic_cast<AlidNdPtAnalysisPbPb*>(lmc->FindObject("dNdPtAnalysisPbPb"));
+ if (!obj) return -1;
+
+ //Event statistics
+ THnSparse *fRecEventMatrix = obj->GetRecEventMatrix(); //all reconstructed events
+ TH2D* h2RecEventAll = fRecEventMatrix->Projection(0,1)->Clone("h2RecEventAll");
+ fRecEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ TH2D* h2RecEvent = fRecEventMatrix->Projection(0,1)->Clone("h2RecEvent");
+ Double_t MCReconstructedEvents = h2RecEvent->Integral();
+ Double_t MCReconstructedEventsAll = h2RecEventAll->Integral();
+ ContHist->Add(h2RecEvent);
+ ContHist->Add(h2RecEventAll);
+
+ THnSparse *fTriggerEventMatrix = obj->GetTriggerEventMatrix(); //all triggered events
+ TH2D* h2TriggerEventAll = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEventAll");
+ fTriggerEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ TH2D* h2TriggerEvent = fTriggerEventMatrix->Projection(0,1)->Clone("h2TriggerEvent");
+ Double_t MCTriggeredEvents = h2TriggerEvent->Integral();
+ Double_t MCTriggeredEventsAll = h2TriggerEventAll->Integral();
+ ContHist->Add(h2TriggerEvent);
+ ContHist->Add(h2TriggerEventAll);
+
+ THnSparse *fGenEventMatrix = obj->GetGenEventMatrix(); //all generated events
+ TH2D* h2GenEventAll = fGenEventMatrix->Projection(0,1)->Clone("h2GenEventAll");
+ fGenEventMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ TH2D* h2GenEvent = fGenEventMatrix->Projection(0,1)->Clone("h2GenEvent");
+ Double_t MCGeneratedEvents = h2GenEvent->Integral();
+ Double_t MCGeneratedEventsAll = h2GenEventAll->Integral();
+ ContHist->Add(h2RecEvent);
+ ContHist->Add(h2RecEventAll);
+
+ printf("=== generated MC events for correction matrices %lf ===\n",MCGeneratedEvents);
+ printf("=== triggered MC events for correction matrices %lf ===\n",MCTriggeredEvents);
+ printf("=== recontructed MC events for correction matrices %lf ===\n",MCReconstructedEvents);
+ printf("\n");
+ printf("=== cut on the zVertex +- %lf ===\n",zVert);
+ printf("=== generated MC events (before zVertex cut) %lf ===\n",MCGeneratedEventsAll);
+ printf("=== triggered MC events (before zVertex cut) %lf ===\n",MCTriggeredEventsAll);
+ printf("=== recontructed MC events (before zVertex cut) %lf ===\n",MCReconstructedEventsAll);
+
+ TH1D* h1MCGeneratedEvents = new TH1D("h1MCGeneratedEvents","h1MCGeneratedEvents",1,0,1);
+ TH1D* h1MCTriggeredEvents = new TH1D("h1MCTriggeredEvents","h1MCTriggeredEvents",1,0,1);
+ TH1D* h1MCReconstructedEvents = new TH1D("h1MCReconstructedEvents","h1MCReconstructedEvents",1,0,1);
+ TH1D* h1MCGeneratedEventsAll = new TH1D("h1MCGeneratedEventsAll","h1MCGeneratedEventsAll",1,0,1);
+ TH1D* h1MCTriggeredEventsAll = new TH1D("h1MCTriggeredEventsAll","h1MCTriggeredEventsAll",1,0,1);
+ TH1D* h1MCReconstructedEventsAll = new TH1D("h1MCReconstructedEventsAll","h1MCReconstructedEventsAll",1,0,1);
+
+ h1MCGeneratedEvents->Fill(0,MCGeneratedEvents);
+ h1MCGeneratedEvents->SetEntries(MCGeneratedEvents);
+ h1MCTriggeredEvents->Fill(0,MCTriggeredEvents);
+ h1MCTriggeredEvents->SetEntries(MCTriggeredEvents);
+ h1MCReconstructedEvents->Fill(0,MCReconstructedEvents);
+ h1MCReconstructedEvents->SetEntries(MCReconstructedEvents);
+ h1MCGeneratedEventsAll->Fill(0,MCGeneratedEventsAll);
+ h1MCGeneratedEventsAll->SetEntries(MCGeneratedEventsAll);
+ h1MCTriggeredEventsAll->Fill(0,MCTriggeredEventsAll);
+ h1MCTriggeredEventsAll->SetEntries(MCTriggeredEventsAll);
+ h1MCReconstructedEventsAll->Fill(0,MCReconstructedEventsAll);
+ h1MCReconstructedEventsAll->SetEntries(MCReconstructedEventsAll);
+
+ ContHist->Add(h1MCGeneratedEvents);
+ ContHist->Add(h1MCTriggeredEvents);
+ ContHist->Add(h1MCReconstructedEvents);
+ ContHist->Add(h1MCGeneratedEventsAll);
+ ContHist->Add(h1MCTriggeredEventsAll);
+ ContHist->Add(h1MCReconstructedEventsAll);
+
+ // efficienfy and correction matrices for tigger and vertex efficiency
+ TH2D* h2EventTriggerEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2GenEventAll,"h2EventTriggerEffAll");
+ TH2D* h2EventTriggerCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2TriggerEventAll,"h2EventTriggerCorrAll");
+ TH2D* h2EventTriggerEff = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2GenEvent,"h2EventTriggerEff");
+ TH2D* h2EventTriggerCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2TriggerEvent,"h2EventTriggerCorr");
+
+ TH2D* h2EventRecEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2TriggerEventAll,"h2EventRecEffAll");
+ TH2D* h2EventRecCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEventAll,h2RecEventAll,"h2EventRecCorrAll");
+ TH2D* h2EventRecEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2TriggerEvent,"h2EventRecEff");
+ TH2D* h2EventRecCorr = AlidNdPtHelper::GenerateCorrMatrix(h2TriggerEvent,h2RecEvent,"h2EventRecCorr");
+
+ TH2D* h2EventEffAll = AlidNdPtHelper::GenerateCorrMatrix(h2RecEventAll,h2GenEventAll,"h2EventEffAll");
+ TH2D* h2EventCorrAll = AlidNdPtHelper::GenerateCorrMatrix(h2GenEventAll,h2RecEventAll,"h2EventCorrAll");
+ TH2D* h2EventEff = AlidNdPtHelper::GenerateCorrMatrix(h2RecEvent,h2GenEvent,"h2EventEff");
+ TH2D* h2EventCorr = AlidNdPtHelper::GenerateCorrMatrix(h2GenEvent,h2RecEvent,"h2EventCorr");
+
+ CorrMatr->Add(h2EventTriggerEffAll);
+ CorrMatr->Add(h2EventTriggerCorrAll);
+ CorrMatr->Add(h2EventTriggerEff);
+ CorrMatr->Add(h2EventTriggerCorr);
+ CorrMatr->Add(h2EventRecEffAll);
+ CorrMatr->Add(h2EventRecCorrAll);
+ CorrMatr->Add(h2EventRecEff);
+ CorrMatr->Add(h2EventRecCorr);
+ CorrMatr->Add(h2EventEffAll);
+ CorrMatr->Add(h2EventCorrAll);
+ CorrMatr->Add(h2EventEff);
+ CorrMatr->Add(h2EventCorr);
+
+
+
+ // all recontructed
+ THnSparse *fRecTrackMatrix = obj->GetRecTrackMatrix();
+ fRecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fRecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+ TH3D* h3RecTrack = fRecTrackMatrix->Projection(0,1,2)->Clone("h3RecTrack");
+ TH2D* h2RecTrack_zv_pt = fRecTrackMatrix->Projection(0,1)->Clone("h2RecTrack_zv_pt");
+ TH2D* h2RecTrack_zv_eta = fRecTrackMatrix->Projection(0,2)->Clone("h2RecTrack_zv_eta");
+ TH2D* h2RecTrack_pt_eta = fRecTrackMatrix->Projection(1,2)->Clone("h2RecTrack_pt_eta");
+ TH1D* h1RecTrack_zv = fRecTrackMatrix->Projection(0)->Clone("h1RecTrack_zv");
+ TH1D* h1RecTrack_pt = fRecTrackMatrix->Projection(1)->Clone("h1RecTrack_pt");
+ TH1D* h1RecTrack_eta = fRecTrackMatrix->Projection(2)->Clone("h1RecTrack_eta");
+ Double_t MCReconstructedTracks = h3RecTrack->Integral();
+
+ ContHist->Add(h3RecTrack);
+ ContHist->Add(h2RecTrack_zv_pt);
+ ContHist->Add(h2RecTrack_zv_eta);
+ ContHist->Add(h2RecTrack_pt_eta);
+ ContHist->Add(h1RecTrack_zv);
+ ContHist->Add(h1RecTrack_pt);
+ ContHist->Add(h1RecTrack_eta);
+
+ // recontructed primary tracks
+ THnSparse *fRecPrimTrackMatrix = obj->GetRecPrimTrackMatrix();
+ THnSparse *fRecTrackMatrixScaled = fRecPrimTrackMatrix->Clone("fRecTrackMatrixScaled"); //used later for secondaries scaling
+ fRecPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fRecPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+ TH3D* h3RecPrimTrack = fRecPrimTrackMatrix->Projection(0,1,2)->Clone("h3RecPrimTrack");
+ TH2D* h2RecPrimTrack_zv_pt = fRecPrimTrackMatrix->Projection(0,1)->Clone("h2RecPrimTrack_zv_pt");
+ TH2D* h2RecPrimTrack_zv_eta = fRecPrimTrackMatrix->Projection(0,2)->Clone("h2RecPrimTrack_zv_eta");
+ TH2D* h2RecPrimTrack_pt_eta = fRecPrimTrackMatrix->Projection(1,2)->Clone("h2RecPrimTrack_pt_eta");
+ TH1D* h1RecPrimTrack_zv = fRecPrimTrackMatrix->Projection(0)->Clone("h1RecPrimTrack_zv");
+ TH1D* h1RecPrimTrack_pt = fRecPrimTrackMatrix->Projection(1)->Clone("h1RecPrimTrack_pt");
+ TH1D* h1RecPrimTrack_eta = fRecPrimTrackMatrix->Projection(2)->Clone("h1RecPrimTrack_eta");
+ Double_t MCReconstructedPrimTracks = h3RecPrimTrack->Integral();
+
+ ContHist->Add(h3RecPrimTrack);
+ ContHist->Add(h2RecPrimTrack_zv_pt);
+ ContHist->Add(h2RecPrimTrack_zv_eta);
+ ContHist->Add(h2RecPrimTrack_pt_eta);
+ ContHist->Add(h1RecPrimTrack_zv);
+ ContHist->Add(h1RecPrimTrack_pt);
+ ContHist->Add(h1RecPrimTrack_eta);
+
+ // recontructed secondary tracks
+ THnSparse *fRecSecTrackMatrix = obj->GetRecSecTrackMatrix();
+ THnSparse *fRecSecTrackMatrixScaled = fRecSecTrackMatrix->Clone("fRecSecTrackMatrixScaled"); //used later for secondaries scaling
+ fRecSecTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fRecSecTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+ TH3D* h3RecSecTrack = fRecSecTrackMatrix->Projection(0,1,2)->Clone("h3RecSecTrack");
+ TH2D* h2RecSecTrack_zv_pt = fRecSecTrackMatrix->Projection(0,1)->Clone("h2RecSecTrack_zv_pt");
+ TH2D* h2RecSecTrack_zv_eta = fRecSecTrackMatrix->Projection(0,2)->Clone("h2RecSecTrack_zv_eta");
+ TH2D* h2RecSecTrack_pt_eta = fRecSecTrackMatrix->Projection(1,2)->Clone("h2RecSecTrack_pt_eta");
+ TH1D* h1RecSecTrack_zv = fRecSecTrackMatrix->Projection(0)->Clone("h1RecSecTrack_zv");
+ TH1D* h1RecSecTrack_pt = fRecSecTrackMatrix->Projection(1)->Clone("h1RecSecTrack_pt");
+ TH1D* h1RecSecTrack_eta = fRecSecTrackMatrix->Projection(2)->Clone("h1RecSecTrack_eta");
+ Double_t MCReconstructedSecTracks = h3RecSecTrack->Integral();
+
+ ContHist->Add(h3RecSecTrack);
+ ContHist->Add(h2RecSecTrack_zv_pt);
+ ContHist->Add(h2RecSecTrack_zv_eta);
+ ContHist->Add(h2RecSecTrack_pt_eta);
+ ContHist->Add(h1RecSecTrack_zv);
+ ContHist->Add(h1RecSecTrack_pt);
+ ContHist->Add(h1RecSecTrack_eta);
+
+ // generated primary tracks
+ THnSparse *fGenPrimTrackMatrix = obj->GetGenPrimTrackMatrix();
+ fGenPrimTrackMatrix->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fGenPrimTrackMatrix->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+ TH3D* h3GenPrimTrack = fGenPrimTrackMatrix->Projection(0,1,2)->Clone("h3GenPrimTrack");
+ TH2D* h2GenPrimTrack_zv_Pt = fGenPrimTrackMatrix->Projection(0,1)->Clone("h2GenPrimTrack_zv_pt");
+ TH2D* h2GenPrimTrack_zv_eta = fGenPrimTrackMatrix->Projection(0,2)->Clone("h2GenPrimTrack_zv_eta");
+ TH2D* h2GenPrimTrack_pt_eta = fGenPrimTrackMatrix->Projection(1,2)->Clone("h2GenPrimTrack_pt_eta");
+ TH1D* h1GenPrimTrack_zv = fGenPrimTrackMatrix->Projection(0)->Clone("h1GenPrimTrack_zv");
+ TH1D* h1GenPrimTrack_pt = fGenPrimTrackMatrix->Projection(1)->Clone("h1GenPrimTrack_pt");
+ TH1D* h1GenPrimTrack_eta = fGenPrimTrackMatrix->Projection(2)->Clone("h1GenPrimTrack_eta");
+ Double_t MCGeneratedPrimTracks = h3GenPrimTrack->Integral();
+
+ ContHist->Add(h3GenPrimTrack);
+ ContHist->Add(h2GenPrimTrack_zv_pt);
+ ContHist->Add(h2GenPrimTrack_zv_eta);
+ ContHist->Add(h2GenPrimTrack_pt_eta);
+ ContHist->Add(h1GenPrimTrack_zv);
+ ContHist->Add(h1GenPrimTrack_pt);
+ ContHist->Add(h1GenPrimTrack_eta);
+ printf("\n");
+ printf("==============================================================\n");
+ printf("=== recontructed MC tracks %lf ===\n",MCReconstructedTracks);
+ printf("=== recontructed MC secondary tracks %lf ===\n",MCReconstructedSecTracks);
+ printf("=== recontructed MC primary tracks %lf ===\n",MCReconstructedPrimTracks);
+ printf("=== generated MC primary track %lf ===\n",MCGeneratedPrimTracks);
+ printf("==============================================================\n");
+ printf("\n");
+
+// THnSparse *fSparseTriggerTrackEvent = obj->GetTriggerTrackEventMatrix();//Tracks from triggered events
+// THnSparse *fSparseVtxTrackEvent = obj->GetRecTrackEventMatrix();//Tracks from events with rec. vtx
+// THnSparse *fSparseGenTrackEvent = obj->GetGenTrackEventMatrix();//generated TrackEvent matrix
+
+ // tracking efficiencies + corrections
+ TH3D* h3TrackEff = AlidNdPtHelper::GenerateCorrMatrix(h3RecPrimTrack,h3GenPrimTrack,"h3TrackEff");
+ TH3D* h3TrackCorr = AlidNdPtHelper::GenerateCorrMatrix(h3GenPrimTrack,h3RecPrimTrack,"h3TrackCorr");
+
+ TH2D* h2TrackEff_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_pt,h2GenPrimTrack_zv_pt,"h2TrackEff_zv_pt");
+ TH2D* h2TrackCorr_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_pt,h2RecPrimTrack_zv_pt,"h2TrackCorr_zv_pt");
+ TH2D* h2TrackEff_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_zv_eta,h2GenPrimTrack_zv_eta,"h2TrackEff_zv_eta");
+ TH2D* h2TrackCorr_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_zv_eta,h2RecPrimTrack_zv_eta,"h2TrackCorr_zv_eta");
+ TH2D* h2TrackEff_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecPrimTrack_pt_eta,h2GenPrimTrack_pt_eta,"h2TrackEff_pt_eta");
+ TH2D* h2TrackCorr_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2GenPrimTrack_pt_eta,h2RecPrimTrack_pt_eta,"h2TrackCorr_pt_eta");
+
+ TH1D* h1TrackEff_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_zv,h1GenPrimTrack_zv,"h1TrackEff_zv");
+ TH1D* h1TrackCorr_zv = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_zv,h1RecPrimTrack_zv,"h1TrackCorr_zv");
+ TH1D* h1TrackEff_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_pt,h1GenPrimTrack_pt,"h1TrackEff_pt");
+ TH1D* h1TrackCorr_pt = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_pt,h1RecPrimTrack_pt,"h1TrackCorr_pt");
+ TH1D* h1TrackEff_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecPrimTrack_eta,h1GenPrimTrack_eta,"h1TrackEff_eta");
+ TH1D* h1TrackCorr_eta = AlidNdPtHelper::GenerateCorrMatrix(h1GenPrimTrack_eta,h1RecPrimTrack_eta,"h1TrackCorr_eta");
+
+ CorrMatr->Add(h3TrackEff);
+ CorrMatr->Add(h3TrackCorr);
+ CorrMatr->Add(h2TrackEff_zv_pt);
+ CorrMatr->Add(h2TrackCorr_zv_pt);
+ CorrMatr->Add(h2TrackEff_zv_eta);
+ CorrMatr->Add(h2TrackCorr_zv_eta);
+ CorrMatr->Add(h2TrackEff_pt_eta);
+ CorrMatr->Add(h2TrackCorr_pt_eta);
+ CorrMatr->Add(h1TrackEff_zv);
+ CorrMatr->Add(h1TrackCorr_zv);
+ CorrMatr->Add(h1TrackEff_pt);
+ CorrMatr->Add(h1TrackCorr_pt);
+ CorrMatr->Add(h1TrackEff_eta);
+ CorrMatr->Add(h1TrackCorr_eta);
+
+ // scale the secondaries before calculating correction matrices
+ for (Long64_t i = 0; i < fRecSecTrackMatrixScaled->GetNbins(); i++) {
+ Int_t c[3];
+ Double_t val = fRecSecTrackMatrixScaled->GetBinContent(i,c);
+ Double_t err = fRecSecTrackMatrixScaled->GetBinError(i);
+ Double_t pt = fRecSecTrackMatrixScaled->GetAxis(1)->GetBinCenter(c[1]);
+ Double_t scale = GetStrangenessCorrFactorPbPb(pt, sscale);
+ fRecSecTrackMatrixScaled->SetBinContent(c,val*scale);
+ fRecSecTrackMatrixScaled->SetBinError(c,err*scale);
+ }
+
+ // for correct determination of secondaries contamination, also the total total tracks have to be scaled
+ // this is done by taking primaries and adding the scaled secondaries
+ fRecTrackMatrixScaled->Add(fRecSecTrackMatrixScaled);
+
+ fRecSecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fRecSecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+
+ TH3D* h3RecSecTrackScaled = fRecSecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecSecTrackScaled");
+ TH2D* h2RecSecTrackScaled_zv_pt = fRecSecTrackMatrixScaled->Projection(0,1)->Clone("h2RecSecTrackScaled_zv_pt");
+ TH2D* h2RecSecTrackScaled_zv_eta = fRecSecTrackMatrixScaled->Projection(0,2)->Clone("h2RecSecTrackScaled_zv_eta");
+ TH2D* h2RecSecTrackScaled_pt_eta = fRecSecTrackMatrixScaled->Projection(1,2)->Clone("h2RecSecTrackScaled_pt_eta");
+ TH1D* h1RecSecTrackScaled_zv = fRecSecTrackMatrixScaled->Projection(0)->Clone("h1RecSecTrackScaled_zv");
+ TH1D* h1RecSecTrackScaled_pt = fRecSecTrackMatrixScaled->Projection(1)->Clone("h1RecSecTrackScaled_pt");
+ TH1D* h1RecSecTrackScaled_eta = fRecSecTrackMatrixScaled->Projection(2)->Clone("h1RecSecTrackScaled_eta");
+
+ ContHist->Add(h3RecSecTrackScaled);
+ ContHist->Add(h2RecSecTrackScaled_zv_pt);
+ ContHist->Add(h2RecSecTrackScaled_zv_eta);
+ ContHist->Add(h2RecSecTrackScaled_pt_eta);
+ ContHist->Add(h1RecSecTrackScaled_zv);
+ ContHist->Add(h1RecSecTrackScaled_pt);
+ ContHist->Add(h1RecSecTrackScaled_eta);
+
+ fRecTrackMatrixScaled->GetAxis(0)->SetRangeUser(-zVert, zVert-0.01);//zVer
+ fRecTrackMatrixScaled->GetAxis(2)->SetRangeUser(-eta, eta-0.01);//eta
+
+ TH3D* h3RecTrackScaled = fRecTrackMatrixScaled->Projection(0,1,2)->Clone("h3RecTrackScaled");
+ TH2D* h2RecTrackScaled_zv_pt = fRecTrackMatrixScaled->Projection(0,1)->Clone("h2RecTrackScaled_zv_pt");
+ TH2D* h2RecTrackScaled_zv_eta = fRecTrackMatrixScaled->Projection(0,2)->Clone("h2RecTrackScaled_zv_eta");
+ TH2D* h2RecTrackScaled_pt_eta = fRecTrackMatrixScaled->Projection(1,2)->Clone("h2RecTrackScaled_pt_eta");
+ TH1D* h1RecTrackScaled_zv = fRecTrackMatrixScaled->Projection(0)->Clone("h1RecTrackScaled_zv");
+ TH1D* h1RecTrackScaled_pt = fRecTrackMatrixScaled->Projection(1)->Clone("h1RecTrackScaled_pt");
+ TH1D* h1RecTrackScaled_eta = fRecTrackMatrixScaled->Projection(2)->Clone("h1RecTrackScaled_eta");
+
+ ContHist->Add(h3RecTrackScaled);
+ ContHist->Add(h2RecTrackScaled_zv_pt);
+ ContHist->Add(h2RecTrackScaled_zv_eta);
+ ContHist->Add(h2RecTrackScaled_pt_eta);
+ ContHist->Add(h1RecTrackScaled_zv);
+ ContHist->Add(h1RecTrackScaled_pt);
+ ContHist->Add(h1RecTrackScaled_eta);
+
+ // create histograms for secondaries contamination and correction
+
+ TH3D* h3SecCont = AlidNdPtHelper::GenerateCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCont");
+ TH3D* h3SecCorr = AlidNdPtHelper::GenerateContCorrMatrix(h3RecSecTrackScaled,h3RecTrackScaled,"h3SecCorr");
+ TH2D* h2SecCont_zv_pt = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCont_zv_pt");
+ TH2D* h2SecCorr_zv_pt = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_pt,h2RecTrackScaled_zv_pt,"h2SecCorr_zv_pt");
+ TH2D* h2SecCont_zv_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCont_zv_eta");
+ TH2D* h2SecCorr_zv_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_zv_eta,h2RecTrackScaled_zv_eta,"h2SecCorr_zv_eta");
+ TH2D* h2SecCont_pt_eta = AlidNdPtHelper::GenerateCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCont_pt_eta");
+ TH2D* h2SecCorr_pt_eta = AlidNdPtHelper::GenerateContCorrMatrix(h2RecSecTrackScaled_pt_eta,h2RecTrackScaled_pt_eta,"h2SecCorr_pt_eta");
+ TH1D* h1SecCont_zv = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCont_zv");
+ TH1D* h1SecCorr_zv = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_zv,h1RecTrackScaled_zv,"h1SecCorr_zv");
+ TH1D* h1SecCont_pt = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCont_pt");
+ TH1D* h1SecCorr_pt = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_pt,h1RecTrackScaled_pt,"h1SecCorr_pt");
+ TH1D* h1SecCont_eta = AlidNdPtHelper::GenerateCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCont_eta");
+ TH1D* h1SecCorr_eta = AlidNdPtHelper::GenerateContCorrMatrix(h1RecSecTrackScaled_eta,h1RecTrackScaled_eta,"h1SecCorr_eta");
+
+ CorrMatr->Add(h3SecCont);
+ CorrMatr->Add(h3SecCorr);
+ CorrMatr->Add(h2SecCont_zv_pt);
+ CorrMatr->Add(h2SecCorr_zv_pt);
+ CorrMatr->Add(h2SecCont_zv_eta);
+ CorrMatr->Add(h2SecCorr_zv_eta);
+ CorrMatr->Add(h2SecCont_pt_eta);
+ CorrMatr->Add(h2SecCorr_pt_eta);
+ CorrMatr->Add(h1SecCont_zv);
+ CorrMatr->Add(h1SecCorr_zv);
+ CorrMatr->Add(h1SecCont_pt);
+ CorrMatr->Add(h1SecCorr_pt);
+ CorrMatr->Add(h1SecCont_eta);
+ CorrMatr->Add(h1SecCorr_eta);
+
+ // plot pictures and save to gifdir
+ for (i=0; i < CorrMatr->LastIndex(); i++) {
+ TCanvas* ctmp = PlotHist(CorrMatr->At(i),idstring);
+ if (gifdir && ctmp) {
+ TString gif(gifdir);
+ gif += '/';
+ gif += ctmp->GetName();
+ gif += ".gif";
+ ctmp->SaveAs(gif.Data(),"gif");
+ delete ctmp;
+ }
+ }
+ for (i=0; i < ContHist->LastIndex(); i++) {
+ TCanvas* ctmp = PlotHist(ContHist->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");
+ CorrMatr->Write();
+ ContHist->Write();
+ out->Close();
+
+ return MCReconstructedEvents;
+
+}
+
+
+//_____________________________________________________________________________
+Double_t GetStrangenessCorrFactorPbPb(Double_t pt, Double_t s)
+{
+ // data driven correction factor for secondaries (PbPb)
+
+ if (pt <= 0.25) return 1.0;
+ if (pt <= 0.5) return GetLinearInterpolationValue(0.25,1.0,0.5,0.60+0.40*s, pt);
+ if (pt <= 1.0) return GetLinearInterpolationValue(0.5,0.60+0.40*s,1.0,0.53+0.47*s, pt);
+ if (pt <= 2.0) return GetLinearInterpolationValue(1.0,0.53+0.47*s,2.0,0.44+0.56*s, pt);
+ if (pt <= 5.0) return GetLinearInterpolationValue(2.0,0.44+0.56*s,5.0,0.33+0.67*s, pt);
+ return 0.33+0.67*s;
+}
+
+
+//___________________________________________________________________________
+Double_t GetLinearInterpolationValue(const Double_t x1,const Double_t y1,const Double_t x2,const Double_t y2, const Double_t pt)
+{
+ //
+ // linear interpolation
+ //
+ return ((y2-y1)/(x2-x1))*pt+(y2-(((y2-y1)/(x2-x1))*x2));
+}
+
+//___________________________________________________________________________
+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(); h->GetXaxis()->SetRangeUser(0.1 , 50.); }
+ }
+ if (h->GetDimension() == 1) {
+ 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;
+
+}
+
+Int_t CheckLoadLibrary(const char* library)
+{
+ // checks if a library is already loaded, if not loads the library
+
+ if (strlen(gSystem->GetLibraries(Form("%s.so", library), "", kFALSE)) > 0)
+ return 1;
+
+ return gSystem->Load(library);
+}