From: rpreghen Date: Fri, 11 May 2012 12:32:42 +0000 (+0000) Subject: TPC rdEdx analysis code, macros and more (P.Christiansen) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=4ebdd20e07ab796e7c73d4a1a08b8162d6e1307b;p=u%2Fmrichter%2FAliRoot.git TPC rdEdx analysis code, macros and more (P.Christiansen) --- diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/README.txt b/PWGLF/SPECTRA/IdentifiedHighPt/README.txt new file mode 100644 index 00000000000..25adb1ee008 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/README.txt @@ -0,0 +1,79 @@ +The idea is to document the different steps needed to make spectra. + +STEP 1: GRID CODE +================= + +In the directory grid you will find the code used to produce the trees. For +info on how to use have a look at runAAF.C + +STEP 2: EXTRACT TREES +===================== + +The file produced on the grid contains the trees in the list. This means that +one cannot direcrtly chain tghem. Therefore we use the code in extract_code to +make new files with the trees only. + +Example: +./merge.sh aortizve/Trees_LHC10b_Pass3/files/ +and +./merge.sh aortizve/Trees_LHC10b_Pass3/files/ HighPtDeDxV0 + +STEP 3: COMPILE LIBRARY +======================= + +cd lib +make clean +make + + +STEP 4: DETERMINE RATIOS +======================== + +This is the biggest step. + +mkdir ratios_7tevb +cd ratios_7tevb + +First we need to create the text files we want to analyze. +Example: +find /home/pchristi/work/analysis/7tev/ | grep HighPtDeDx_Tree | grep new | grep 117059 > 7tev_b_test.dat +find /home/pchristi/work/analysis/7tev/ | grep HighPtDeDx_Tree | grep new > 7tev_b.dat + +find /home/pchristi/work/analysis/7tev/ | grep HighPtDeDxV0_Tree | grep new | grep 117059 > 7tev_b_test_v0.dat +find /home/pchristi/work/analysis/7tev/ | grep HighPtDeDxV0_Tree | grep new > 7tev_b_v0.dat + +ln -s ../macros/run_code.C . +ln -s ../macros/calibrate_de_dx.C . + +cp ../macros/drawText.C . +Edit the text here. This macro is used to tag the pictures. + +Follow the example in the macro run_code.C + +Step 1-5 is about determining the dE/dx calibrations and the input data to the +fits in pT. + +Now it is time for extracting the uncorrected ratios. + +ln -s ../macros/fit_yields_final.C . + +This is documented in fit_yields_final.C + + +There is still things missing: +- Option to generate tree when generating the data. +- The code to estimate systematic errors +- Efficiency code +- Corrected fractions +- Spectra and RAA code + +But that will come soon - hopefully next week. + + + + +To zip: + +ls -1 README.txt lib/Makefile lib/*.cxx lib/*.h macros/*.C > tozip.txt +tar -hcvzf analysis.tgz -T tozip.txt + diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/corrected_fraction.C b/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/corrected_fraction.C new file mode 100644 index 00000000000..ef070b3efb7 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/corrected_fraction.C @@ -0,0 +1,195 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "my_tools.C" + + +/* + Info: + * Efficiencies are currently set by hand!!!!! + * Systematic errors needs an update also. + * Need to accomodate k+p at soem point! + + To run: + + gSystem->AddIncludePath("-I../macros") + gSystem->AddIncludePath("-I../lib") + gROOT->SetMacroPath(".:../macros") + .L my_tools.C++ + .L corrected_fraction.C+ + + + CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_pi.root", 1, "PP", 0, 10.0); + CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_k.root", 2, "PP", 0, 10.0); + CorrectedFraction("../ratios_7tevb/fit_yields_results/7tev_b.root", "fraction_7tev_b_p.root", 3, "PP", 0, 10.0); + + + And in the shell one can add files like: + + hadd fraction_7tev_b_all.root fraction_7tev_b_pi.root fraction_7tev_b_k.root fraction_7tev_b_p.root + + */ + +void CorrectedFraction(const Char_t* inFileName, const Char_t* outFileName, Int_t pid, + const Char_t* endName="PP", Int_t centBin=0, + Double_t ptMax=20.0, const Char_t* outDir="plots"); +TH1D* GetSystErrorHist(TH1D* hRatio, Int_t centBin, TF1* electronFraction); + + +//________________________________________________________________________________________ +void CorrectedFraction(const Char_t* inFileName, const Char_t* outFileName, Int_t pid, + const Char_t* endName, Int_t centBin, + Double_t ptMax, const Char_t* outDir) +{ + gStyle->SetOptStat(0); + + + TFile* fileData = FindFileFresh(inFileName); + + if(!fileData) + return; + + TF1* fElectronFraction = (TF1*)fileData->Get("fElectronFraction"); + cout << "Electron fraction found!" << endl; + TH1D* hPionRatio = 0; + if(pid==1) { + hPionRatio = (TH1D*)fileData->Get("hPionRatio"); + CutHistogram(hPionRatio, 2.0, ptMax); + hPionRatio->SetMarkerStyle(24); + } else if(pid==2) { + hPionRatio = (TH1D*)fileData->Get("hKaonRatio"); + CutHistogram(hPionRatio, 3.0, ptMax); + hPionRatio->SetMarkerStyle(25); + } else if(pid==3) { + hPionRatio = (TH1D*)fileData->Get("hProtonRatio"); + CutHistogram(hPionRatio, 3.0, ptMax); + hPionRatio->SetMarkerStyle(25); + } + + // Global variable + TH1D* hSystFraction = 0; + if(pid==1) + hSystFraction = GetSystErrorHist(hPionRatio, centBin, fElectronFraction); + else + hSystFraction = GetSystErrorHist(hPionRatio, centBin, 0); + hSystFraction->SetLineColor(1); + hSystFraction->SetMarkerStyle(1); + hSystFraction->SetFillStyle(1001); + hSystFraction->SetFillColor(kGray); + + TH1D* hPionFraction = (TH1D*)hPionRatio->Clone("hPionFraction"); + if(pid==1 && !fElectronFraction) { + TF1 f1("f1", "1.0", 0.0, 50.0); + cout << "NO ELECTRON FRACTION!!!" << endl; + hPionFraction->Add(&f1, -0.01); // correct for muons and electrons + CutHistogram(hPionFraction, 3.0, ptMax); + hSystFraction->Add(&f1, -0.01); // correct for muons and electrons + CutHistogram(hSystFraction, 3.0, ptMax); + } + + if(pid==1 || pid ==3) { + + hPionFraction->Scale(0.94); // correct for efficiency ratio + hSystFraction->Scale(0.94); // correct for efficiency ratio + } else { + + TF1* effRatioK = new TF1("effRatioK", "exp(-[1]*x)+[0]", 0.0, 50.0); + effRatioK->SetParameters(9.82065e-01, 1.28157e+00); + hPionFraction->Multiply(effRatioK); // correct for efficiency ratio + hSystFraction->Multiply(effRatioK); // correct for efficiency ratio + } + + if(pid==1) + hPionFraction->SetMarkerStyle(20); + else + hPionFraction->SetMarkerStyle(20); + TLatex latex; + latex.SetNDC(); + latex.SetTextSize(0.05); + + TCanvas* cRatios = new TCanvas("cRatios", "Particle fractions"); + cRatios->Clear(); + hSystFraction->GetXaxis()->SetRangeUser(0, ptMax); + hSystFraction->GetYaxis()->SetRangeUser(0, 1); + hSystFraction->Draw("E5"); + hPionRatio->Draw("SAME"); + hPionFraction->Draw("SAME"); + latex.DrawLatex(0.6, 0.95, Form("%s", endName)); + + CreateDir(outDir); + if(pid==1) { + hPionFraction->SetName(Form("hPionFraction_%s", endName)); + hSystFraction->SetName(Form("hPionFractionSyst_%s", endName)); + cRatios->SaveAs(Form("%s/fractions_%s_pions.gif", outDir, endName)); + } else if (pid==2) { + hPionFraction->SetName(Form("hKaonFraction_%s", endName)); + hSystFraction->SetName(Form("hKaonFractionSyst_%s", endName)); + cRatios->SaveAs(Form("%s/fractions_%s_kaons.gif", outDir, endName)); + } else if (pid==3) { + hPionFraction->SetName(Form("hProtonFraction_%s", endName)); + hSystFraction->SetName(Form("hProtonFractionSyst_%s", endName)); + cRatios->SaveAs(Form("%s/fractions_%s_protons.gif", outDir, endName)); + } + + TFile* outFile = new TFile(outFileName, "RECREATE"); + hPionFraction->Write(); + hSystFraction->Write(); + outFile->Close(); +} + +//__________________________________________________________________________________ +TH1D* GetSystErrorHist(TH1D* hRatio, Int_t centBin, TF1* electronFraction) +{ + TFile* fileSyst = FindFile("syst_vs_pt.root"); + + TF1* fSyst = (TF1*)fileSyst->Get(Form("systHigh%d", centBin)); + fSyst->SetRange(2.0, 20.0); + fSyst->Print(); + + Double_t syst_error_mc = 0.03; // pp + if (centBin>0) + syst_error_mc = 0.05; // Pb+Pb + + Double_t syst_error_correction = 0.01; // pp + if (centBin>0) + syst_error_correction = 0.03; // Pb+Pb + + TH1D* hSystError = (TH1D*)hRatio->Clone("hPionFractionSyst"); + hSystError->Multiply(fSyst); + + const Int_t nBins = hSystError->GetNbinsX(); + for(Int_t bin = 1; bin <= nBins; bin++) { + + Double_t value = hRatio->GetBinContent(bin); + Double_t stat_error = hSystError->GetBinContent(bin); + + if(value==0) + continue; + + Double_t syst_error = stat_error*stat_error; + syst_error += value*value*syst_error_mc*syst_error_mc; + syst_error += value*value*syst_error_correction*syst_error_correction; + + + if (electronFraction) { + + Double_t systEandMu = electronFraction->Eval(hRatio->GetXaxis()->GetBinCenter(bin)); + syst_error += systEandMu*systEandMu; + } + + syst_error = TMath::Sqrt(syst_error); + hSystError->SetBinContent(bin, value); + hSystError->SetBinError(bin, syst_error); + } + + return hSystError; +} + diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/syst_vs_pt.root b/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/syst_vs_pt.root new file mode 100644 index 00000000000..bbfe6dab393 Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/corrected_fraction/syst_vs_pt.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/createEfficiency.C b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/createEfficiency.C new file mode 100644 index 00000000000..1b48ea74cb9 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/createEfficiency.C @@ -0,0 +1,697 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "DebugClasses.C" +#include "my_tools.C" + +#include + +using namespace std; + +/* + To run code: + ============ + + Info: + * I did not recheck this code. For now I would just use the efficiency values I have. + * The code could be made nicer. Esepcially some plots could be put in folders. + + Use AliRoot because of AliXRDPROOFtoolkit: + gSystem->AddIncludePath("-I$ALICE_ROOT/TPC") + gSystem->AddIncludePath("-I../lib") + gSystem->AddIncludePath("-I../grid") + gSystem->AddIncludePath("-I../macros") + gROOT->SetMacroPath(".:../macros:../grid:../lib/") + .L my_tools.C+ + .L DebugClasses.C+ + .L createEfficiency.C+ + + Examples of visualization: + DrawEfficiency("lhc10d_eff_pythia.root", "eff.root") + + // This is the correction I did for the low pT guys + DrawCorrection("lhc10d_eff_pythia.root", "lhc10d_eff_phojet.root") + + + + CreateEff("lhc10d_mc_all.dat", 0, "lhc10d_eff_all.root") + + +*/ + +void DrawEfficiency(const Char_t* fileName, const Char_t* effFileName); +void DrawCorrection(const Char_t* fileNamePYTHIA, const Char_t* fileNamePHOJET); +TH1D* HistInvert(TH1D* hist); + + +void CreateEff(const Char_t* mcFileName, Int_t maxEvents, const Char_t* mcOutFileName, + Float_t centLow=-20, Float_t centHigh=-5) +{ + gStyle->SetOptStat(0); + // + // Create output + // + TFile* outFile = new TFile(mcOutFileName, "RECREATE"); + + const Int_t nPid = 7; + TH1D* hMcIn[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcOut[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcSec[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcEff[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcInNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcOutNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcSecNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcEffNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcInPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcOutPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcSecPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hMcEffPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + + Int_t color[nPid] = {1, 2, 3, 4, 5, 1, 1}; + + const Int_t nPtBins = 68; + Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, + 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, + 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 , + 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 , + 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0, + 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, + 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 }; + + for(Int_t pid = 0; pid < nPid; pid++) { + + hMcIn[pid] = new TH1D(Form("hIn%d", pid), Form("Efficiency (pid %d)", pid), + nPtBins, xBins); + hMcInNeg[pid] = new TH1D(Form("hInNeg%d", pid), Form("Efficiency (pid %d, q < 0)", pid), + nPtBins, xBins); + hMcInPos[pid] = new TH1D(Form("hInPos%d", pid), Form("Efficiency (pid %d, q < 0)", pid), + nPtBins, xBins); + + hMcIn[pid]->Sumw2(); + hMcIn[pid]->SetMarkerStyle(29); + hMcIn[pid]->SetMarkerColor(color[pid]); + hMcInNeg[pid]->Sumw2(); + hMcInNeg[pid]->SetMarkerStyle(24); + hMcInNeg[pid]->SetMarkerColor(color[pid]); + hMcInPos[pid]->Sumw2(); + hMcInPos[pid]->SetMarkerStyle(20); + hMcInPos[pid]->SetMarkerColor(color[pid]); + + hMcOut[pid] = new TH1D(Form("hMcOut%d", pid), Form("MC out (pid %d)", pid), + nPtBins, xBins); + hMcOutNeg[pid] = new TH1D(Form("hMcOutNeg%d", pid), Form("MC out (pid %d, q < 0)", pid), + nPtBins, xBins); + hMcOutPos[pid] = new TH1D(Form("hMcOutPos%d", pid), Form("MC out (pid %d, q < 0)", pid), + nPtBins, xBins); + + hMcOut[pid]->Sumw2(); + hMcOut[pid]->SetMarkerStyle(29); + hMcOut[pid]->SetMarkerColor(color[pid]); + hMcOutNeg[pid]->Sumw2(); + hMcOutNeg[pid]->SetMarkerStyle(24); + hMcOutNeg[pid]->SetMarkerColor(color[pid]); + hMcOutPos[pid]->Sumw2(); + hMcOutPos[pid]->SetMarkerStyle(20); + hMcOutPos[pid]->SetMarkerColor(color[pid]); + + hMcSec[pid] = new TH1D(Form("hSec%d", pid), Form("Secondaries (pid %d)", pid), + nPtBins, xBins); + hMcSecNeg[pid] = new TH1D(Form("hSecNeg%d", pid), Form("Secondaries (pid %d, q < 0)", pid), + nPtBins, xBins); + hMcSecPos[pid] = new TH1D(Form("hSecPos%d", pid), Form("Secondaries (pid %d, q < 0)", pid), + nPtBins, xBins); + + hMcSec[pid]->Sumw2(); + hMcSec[pid]->SetMarkerStyle(29); + hMcSec[pid]->SetMarkerColor(color[pid]); + hMcSecNeg[pid]->Sumw2(); + hMcSecNeg[pid]->SetMarkerStyle(24); + hMcSecNeg[pid]->SetMarkerColor(color[pid]); + hMcSecPos[pid]->Sumw2(); + hMcSecPos[pid]->SetMarkerStyle(20); + hMcSecPos[pid]->SetMarkerColor(color[pid]); + + hMcEff[pid] = new TH1D(Form("hEff%d", pid), Form("Efficiency (pid %d)", pid), + nPtBins, xBins); + hMcEffNeg[pid] = new TH1D(Form("hEffNeg%d", pid), Form("Efficiency (pid %d, q < 0)", pid), + nPtBins, xBins); + hMcEffPos[pid] = new TH1D(Form("hEffPos%d", pid), Form("Efficiency (pid %d, q < 0)", pid), + nPtBins, xBins); + + hMcEff[pid]->Sumw2(); + hMcEff[pid]->SetMarkerStyle(29); + hMcEff[pid]->SetMarkerColor(color[pid]); + hMcEffNeg[pid]->Sumw2(); + hMcEffNeg[pid]->SetMarkerStyle(24); + hMcEffNeg[pid]->SetMarkerColor(color[pid]); + hMcEffPos[pid]->Sumw2(); + hMcEffPos[pid]->SetMarkerStyle(20); + hMcEffPos[pid]->SetMarkerColor(color[pid]); + } + + TTree* Tree = 0; + if(strstr(mcFileName, ".dat")) { + + AliXRDPROOFtoolkit tool; + TChain* chain = tool.MakeChain(mcFileName,"tree", 0, 1000); + chain->Lookup(); + Tree = chain; + } else { + TFile* mcFile = FindFileFresh(mcFileName); + if(!mcFile) + return; + + Tree = (TTree*)mcFile->Get("tree"); + } + if(!Tree) + return; + + + DeDxEvent* event = 0; + TClonesArray* trackArray = 0; + TClonesArray* mcTrackArray = 0; + Tree->SetBranchAddress("event", &event); + Tree->SetBranchAddress("track" , &trackArray); + Tree->SetBranchAddress("trackMC" , &mcTrackArray); + Int_t nEvents = Tree->GetEntries(); + cout << "Number of events: " << nEvents << endl; + + if(maxEvents>0 && maxEvents < nEvents) { + + nEvents = maxEvents; + cout << "N events was reduced to: " << maxEvents << endl; + } + + Int_t currentRun = 0; + + for(Int_t n = 0; n < nEvents; n++) { + + Tree->GetEntry(n); + + if((n+1)%1000000==0) + cout << "Event: " << n+1 << "/" << nEvents << endl; + + if(event->run != currentRun) { + + cout << "New run: " << event->run << endl; + currentRun = event->run; + } + + if(event->cent < centLow || event->cent > centHigh) + continue; + + const Int_t nMcTracks = mcTrackArray->GetEntries(); + + for(Int_t i = 0; i < nMcTracks; i++) { + + DeDxTrackMC* trackMC = (DeDxTrackMC*)mcTrackArray->At(i); + + // if(TMath::Abs(trackMC->pdgMC)==3312 || TMath::Abs(trackMC->pdgMC)==3334) + // continue; // Xi-! + + hMcIn[0]->Fill(trackMC->ptMC); + hMcIn[trackMC->pidMC]->Fill(trackMC->ptMC); + + if(trackMC->qMC < 0) { + + hMcInNeg[0]->Fill(trackMC->ptMC); + hMcInNeg[trackMC->pidMC]->Fill(trackMC->ptMC); + } else { + + hMcInPos[0]->Fill(trackMC->ptMC); + hMcInPos[trackMC->pidMC]->Fill(trackMC->ptMC); + } + } + + const Int_t nTracks = trackArray->GetEntries(); + + for(Int_t i = 0; i < nTracks; i++) { + + DeDxTrack* track = (DeDxTrack*)trackArray->At(i); + + if(!(track->filter&1)) + continue; + + // if(TMath::Abs(track->mother)==3312 || TMath::Abs(track->mother)==3334) + // continue; // Xi+- or Omega+-! + + hMcOut[0]->Fill(track->pt); + hMcOut[track->pid]->Fill(track->pt); + + if(track->q < 0) { + + hMcOutNeg[0]->Fill(track->pt); + hMcOutNeg[track->pid]->Fill(track->pt); + } else { + + hMcOutPos[0]->Fill(track->pt); + hMcOutPos[track->pid]->Fill(track->pt); + } + + if(track->primary==0) { + hMcSec[0]->Fill(track->pt); + hMcSec[track->pid]->Fill(track->pt); + + if(track->q < 0) { + + hMcSecNeg[0]->Fill(track->pt); + hMcSecNeg[track->pid]->Fill(track->pt); + } else { + + hMcSecPos[0]->Fill(track->pt); + hMcSecPos[track->pid]->Fill(track->pt); + } + } + } + } + + TH1D* hMcInPiKP = (TH1D*)hMcIn[1]->Clone("hMcInPiKP"); + hMcInPiKP->Add(hMcIn[2]); + hMcInPiKP->Add(hMcIn[3]); + hMcInPiKP->Divide(hMcIn[0]); + + for(Int_t pid = 0; pid < nPid; pid++) { + + hMcSec[pid]->Divide(hMcSec[pid], hMcOut[pid]); + hMcSecNeg[pid]->Divide(hMcSecNeg[pid], hMcOutNeg[pid]); + hMcSecPos[pid]->Divide(hMcSecPos[pid], hMcOutPos[pid]); + + hMcEff[pid] ->Divide(hMcOut[pid], hMcIn[pid]); + hMcEffNeg[pid]->Divide(hMcOutNeg[pid], hMcInNeg[pid]); + hMcEffPos[pid]->Divide(hMcOutPos[pid], hMcInPos[pid]); + } + + outFile->Write(); + outFile->Close(); +} + +//_______________________________________________________________________ +void DrawEfficiency(const Char_t* fileName, const Char_t* effFileName) +{ + TFile* file = FindFileFresh(fileName); + if(!file) + return; + + gStyle->SetOptStat(0); + + const Double_t ptmin = 3.01; + const Double_t ptmax = 19.999; + // const Double_t ptmax = 49.999; + + const Double_t effmin = 0.6; + const Double_t effmax = 0.9; + + // const Double_t effmin = 0.0; + // const Double_t effmax = 1.5; + + const Double_t secmin = 0.0; + const Double_t secmax = 0.1; + + const Int_t nPid = 7; + + TH1D* hEff[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hEffNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hEffPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hSec[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hSecNeg[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hSecPos[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hOut[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + + for(Int_t pid = 0; pid < nPid; pid++) { + + hEff[pid] = (TH1D*)file->Get(Form("hEff%d", pid)); + + hEffNeg[pid] = (TH1D*)file->Get(Form("hEffNeg%d", pid)); + + hEffPos[pid] = (TH1D*)file->Get(Form("hEffPos%d", pid)); + + hSec[pid] = (TH1D*)file->Get(Form("hSec%d", pid)); + + hSecNeg[pid] = (TH1D*)file->Get(Form("hSecNeg%d", pid)); + + hSecPos[pid] = (TH1D*)file->Get(Form("hSecPos%d", pid)); + + hOut[pid] = (TH1D*)file->Get(Form("hMcOut%d", pid)); + + // hEff[pid]->Rebin(4); + // hEffNeg[pid]->Rebin(4); + // hEffPos[pid]->Rebin(4); + + // hSec[pid]->Rebin(4); + // hSecNeg[pid]->Rebin(4); + // hSecPos[pid]->Rebin(4); + + // hOut[pid]->Rebin(4); + + // hEff[pid]->Scale(0.25); + // hEffNeg[pid]->Scale(0.25); + // hEffPos[pid]->Scale(0.25); + + // hSec[pid]->Scale(0.25); + // hSecNeg[pid]->Scale(0.25); + // hSecPos[pid]->Scale(0.25); + + // hOut[pid]->Scale(0.25); + + hEff[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hEff[pid]->GetYaxis()->SetRangeUser(effmin, effmax); + hEffNeg[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hEffNeg[pid]->GetYaxis()->SetRangeUser(effmin, effmax); + hEffPos[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hEffPos[pid]->GetYaxis()->SetRangeUser(effmin, effmax); + hSec[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hSec[pid]->GetYaxis()->SetRangeUser(secmin, secmax); + hSecNeg[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hSecNeg[pid]->GetYaxis()->SetRangeUser(secmin, secmax); + hSecPos[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hSecPos[pid]->GetYaxis()->SetRangeUser(secmin, secmax); + hOut[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + } + + TCanvas* cChEff = new TCanvas("cChEff", "All efficiency", 800, 300); + cChEff->Clear(); + cChEff->Divide(2, 1); + cChEff->cd(1); + // TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0); + TF1* chEff = new TF1("chEff", "[0]*(1-[1]/x)", 0.0, 50.0); + chEff->SetLineColor(1); + chEff->SetParameters(0.7, 0.5); + hEff[0]->Fit(chEff, "", "", ptmin, ptmax); + + cChEff->cd(2); + hEffPos[0]->Draw(); + hEffNeg[0]->Draw("SAME"); + chEff->DrawCopy("SAME"); + cChEff->SaveAs("eff_charged.gif"); + + TCanvas* cPionEff = new TCanvas("cPionEff", "Pion efficiency", 800, 300); + cPionEff->Clear(); + cPionEff->Divide(2, 1); + cPionEff->cd(1); + // TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0); + TF1* pionEff = new TF1("pionEff", "[0]*(1-[1]/x)", 0.0, 50.0); + pionEff->SetLineColor(2); + pionEff->SetParameters(0.7, 0.5); + hEff[1]->Fit(pionEff, "", "", ptmin, ptmax); + + cPionEff->cd(2); + hEffPos[1]->Draw(); + hEffNeg[1]->Draw("SAME"); + pionEff->DrawCopy("SAME"); + cPionEff->SaveAs("eff_pion.gif"); + + TCanvas* cKaonEff = new TCanvas("cKaonEff", "Kaon efficiency", 800, 300); + cKaonEff->Clear(); + cKaonEff->Divide(2, 1); + cKaonEff->cd(1); + // TF1* kaonEff = new TF1("kaonEff", "pol0", 0.0, 50.0); + TF1* kaonEff = new TF1("kaonEff", "[0]*(1-[1]/x)", 0.0, 50.0); + kaonEff->SetLineColor(3); + kaonEff->SetParameters(0.7, 0.5); + hEff[2]->Fit(kaonEff, "", "", ptmin, ptmax); + + cKaonEff->cd(2); + hEffPos[2]->Draw(); + hEffNeg[2]->Draw("SAME"); + kaonEff->DrawCopy("SAME"); + cKaonEff->SaveAs("eff_kaon.gif"); + + TCanvas* cProtonEff = new TCanvas("cProtonEff", "Proton efficiency", 800, 300); + cProtonEff->Clear(); + cProtonEff->Divide(2, 1); + cProtonEff->cd(1); + // TF1* protonEff = new TF1("protonEff", "pol0", 0.0, 50.0); + TF1* protonEff = new TF1("protonEff", "[0]*(1-[1]/x)", 0.0, 50.0); + protonEff->SetLineColor(4); + protonEff->SetParameters(0.7, 0.5); + hEff[3]->Fit(protonEff, "", "", ptmin, ptmax); + + cProtonEff->cd(2); + hEffPos[3]->Draw(); + hEffNeg[3]->Draw("SAME"); + protonEff->DrawCopy("SAME"); + cProtonEff->SaveAs("eff_proton.gif"); + + TCanvas* cAllEff = new TCanvas("cAllEff", "All efficiency", 400, 300); + cAllEff->Clear(); + + TH1D* hDummy = (TH1D*)hEff[1]->Clone("hDummy"); + hDummy->Reset(); + hDummy->SetTitle("Efficiency vs p_{T}; p_{T} [GeV/c]; Efficiency"); + hDummy->Draw(); + + chEff->DrawCopy("SAME"); + pionEff->DrawCopy("SAME"); + kaonEff->DrawCopy("SAME"); + protonEff->DrawCopy("SAME"); + cAllEff->SaveAs("eff_all.gif"); + + // TCanvas* cPionSec = new TCanvas("cPionSec", "Pion sec", 800, 300); + // cPionSec->Clear(); + // cPionSec->Divide(2, 1); + // cPionSec->cd(1); + // TF1* pionSec = new TF1("pionSec", "pol0", 0.0, 50.0); + // hSec[1]->Fit(pionSec, "", "", ptmin, ptmax); + + // cPionSec->cd(2); + // hSecPos[1]->Draw(); + // hSecNeg[1]->Draw("SAME"); + // pionSec->Draw("SAME"); + // cPionSec->SaveAs("sec_pion.gif"); + + // TCanvas* cKaonSec = new TCanvas("cKaonSec", "Kaon sec", 800, 300); + // cKaonSec->Clear(); + // cKaonSec->Divide(2, 1); + // cKaonSec->cd(1); + // TF1* kaonSec = new TF1("kaonSec", "pol0", 0.0, 50.0); + // hSec[2]->Fit(kaonSec, "", "", ptmin, ptmax); + + // cKaonSec->cd(2); + // hSecPos[2]->Draw(); + // hSecNeg[2]->Draw("SAME"); + // kaonSec->Draw("SAME"); + // cKaonSec->SaveAs("sec_kaon.gif"); + + // TCanvas* cProtonSec = new TCanvas("cProtonSec", "Proton sec", 800, 300); + // cProtonSec->Clear(); + // cProtonSec->Divide(2, 1); + // cProtonSec->cd(1); + // TF1* protonSec = new TF1("protonSec", "pol0", 0.0, 50.0); + // hSec[3]->Fit(protonSec, "", "", ptmin, ptmax); + + // cProtonSec->cd(2); + // hSecPos[3]->Draw(); + // hSecNeg[3]->Draw("SAME"); + // protonSec->Draw("SAME"); + // cProtonSec->SaveAs("sec_proton.gif"); + + + TCanvas* cEffRatioPi = new TCanvas("cEffRatioPi", "eff pi / eff all", 400, 300); + cEffRatioPi->Clear(); + // TF1* pionEff = new TF1("pionEff", "pol0", 0.0, 50.0); + TF1* effRatioPi = new TF1("effRatioPi", "pol0", 0.0, 50.0); + effRatioPi->SetLineColor(6); + effRatioPi->SetParameters(0.7, 0.5); + TH1D* hEffRatioPi = (TH1D*)hEff[1]->Clone("hEffRatioPi"); + hEffRatioPi->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{pi}"); + hEffRatioPi->Divide(hEff[0], hEff[1], 1, 1, "B"); + // hEffRatioPi->Divide(hEff[1], hEff[0]); + hEffRatioPi->GetXaxis()->SetRangeUser(0.0, 19.99); + hEffRatioPi->GetYaxis()->SetRangeUser(0.8, 1.1); + hEffRatioPi->SetStats(kTRUE); + hEffRatioPi->Fit(effRatioPi, "", "", ptmin, ptmax); + cEffRatioPi->SaveAs("eff_ratio_pi.gif"); + cEffRatioPi->SaveAs("eff_ratio_pi.pdf"); + + TCanvas* cEffRatioK = new TCanvas("cEffRatioK", "eff K / eff all", 400, 300); + cEffRatioK->Clear(); + TF1* effRatioK = new TF1("effRatioK", "exp(-[1]*x)+[0]", 0.0, 50.0); + effRatioK->SetParameters(1.0, 1.0); + effRatioK->SetLineColor(6); + TH1D* hEffChargedRebinned = (TH1D*)hEff[0]->Clone("hEffChargedRebinned"); + hEffChargedRebinned->Rebin(2); + TH1D* hEffRatioK = (TH1D*)hEff[2]->Clone("hEffRatioK"); + hEffRatioK->Rebin(2); + hEffRatioK->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{K}"); + hEffRatioK->Divide(hEffChargedRebinned, hEffRatioK, 1, 1, "B"); + // hEffRatioK->Divide(hEff[1], hEff[0]); + hEffRatioK->GetXaxis()->SetRangeUser(0.0, 19.99); + hEffRatioK->GetYaxis()->SetRangeUser(0.8, 1.1); + hEffRatioK->SetStats(kTRUE); + hEffRatioK->Fit(effRatioK, "", "", ptmin, ptmax); + cEffRatioK->SaveAs("eff_ratio_K.gif"); + cEffRatioK->SaveAs("eff_ratio_K.pdf"); + + TCanvas* cEffRatioP = new TCanvas("cEffRatioP", "eff p / eff all", 400, 300); + cEffRatioP->Clear(); + TF1* effRatioP = new TF1("effRatioP", "pol0", 0.0, 50.0); + effRatioP->SetLineColor(6); + effRatioP->SetParameters(0.7, 0.5); + // TH1D* hEffChargedRebinned = (TH1D*)hEff[0]->Clone("hEffChargedRebinned"); + // hEffChargedRebinned->Rebin(2); + TH1D* hEffRatioP = (TH1D*)hEff[3]->Clone("hEffRatioP"); + hEffRatioP->Rebin(2); + hEffRatioP->SetTitle("; p_{T} [GeV/c]; #epsilon_{ch}/#epsilon_{p}"); + hEffRatioP->Divide(hEffChargedRebinned, hEffRatioP, 1, 1, "B"); + // hEffRatioP->Divide(hEff[1], hEff[0]); + hEffRatioP->GetXaxis()->SetRangeUser(0.0, 19.99); + hEffRatioP->GetYaxis()->SetRangeUser(0.8, 1.1); + hEffRatioP->SetStats(kTRUE); + hEffRatioP->Fit(effRatioP, "", "", ptmin, ptmax); + cEffRatioP->SaveAs("eff_ratio_p.gif"); + cEffRatioP->SaveAs("eff_ratio_p.pdf"); + + + // TCanvas* cMuonRatio = new TCanvas("cMuonRatio", "muon / all", 400, 300); + // cMuonRatio->Clear(); + // // TF1* pionMuon = new TF1("pionMuon", "pol0", 0.0, 50.0); + // TF1* muonRatio = new TF1("muonRatio", "pol0", 0.0, 50.0); + // muonRatio->SetLineColor(1); + // muonRatio->SetParameter(0, 0.01); + // TH1D* hMuonRatio = (TH1D*)hOut[5]->Clone("hMuonRatio"); + // hMuonRatio->SetTitle("; p_{T} [GeV/c]; muon/all"); + // hMuonRatio->Divide(hOut[5], hOut[0], 1, 1, "B"); + // hMuonRatio->GetYaxis()->SetRangeUser(0.0, 0.05); + // hMuonRatio->SetStats(kTRUE); + // hMuonRatio->Fit(muonRatio, "", "", ptmin, ptmax); + // cMuonRatio->SaveAs("muon_ratio.gif"); + // cMuonRatio->SaveAs("muon_ratio.pdf"); + + // TCanvas* cElectronRatio = new TCanvas("cElectronRatio", "electron / all", 400, 300); + // cElectronRatio->Clear(); + // // TF1* pionElectron = new TF1("pionElectron", "pol0", 0.0, 50.0); + // TF1* electronRatio = new TF1("electronRatio", "pol0", 0.0, 50.0); + // electronRatio->SetLineColor(1); + // electronRatio->SetParameter(0, 0.01); + // TH1D* hElectronRatio = (TH1D*)hOut[4]->Clone("hElectronRatio"); + // hElectronRatio->SetTitle("; p_{T} [GeV/c]; electron/all"); + // hElectronRatio->Divide(hOut[4], hOut[0], 1, 1, "B"); + // hElectronRatio->GetYaxis()->SetRangeUser(0.0, 0.05); + // hElectronRatio->SetStats(kTRUE); + // hElectronRatio->SetMarkerColor(6); + // hElectronRatio->Fit(electronRatio, "", "", ptmin, ptmax); + // cElectronRatio->SaveAs("electron_ratio.gif"); + // cElectronRatio->SaveAs("electron_ratio.pdf"); + + + TFile* effFile = new TFile(effFileName, "RECREATE"); + chEff->Write(); + pionEff->Write(); + kaonEff->Write(); + protonEff->Write(); + effFile->Close(); +} + + +//_______________________________________________________________________ +void DrawCorrection(const Char_t* fileNamePYTHIA, const Char_t* fileNamePHOJET) +{ + TFile* filePYTHIA = FindFileFresh(fileNamePYTHIA); + if(!filePYTHIA) + return; + + TFile* filePHOJET = FindFileFresh(fileNamePHOJET); + if(!filePHOJET) + return; + + gStyle->SetOptStat(0); + + const Double_t ptmin = 0.0; + const Double_t ptmax = 4.999; + // const Double_t ptmax = 49.999; + + const Double_t effmin = 0.95; + const Double_t effmax = 1.12; + + const Int_t nPid = 7; + + TH1D* hInPYTHIA[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + TH1D* hInPHOJET[nPid] = {0, 0, 0, 0, 0, 0, 0 }; + + for(Int_t pid = 0; pid < nPid; pid++) { + + hInPYTHIA[pid] = (TH1D*)filePYTHIA->Get(Form("hIn%d", pid)); + hInPYTHIA[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hInPYTHIA[pid]->GetYaxis()->SetRangeUser(effmin, effmax); + + hInPHOJET[pid] = (TH1D*)filePHOJET->Get(Form("hIn%d", pid)); + hInPHOJET[pid]->GetXaxis()->SetRangeUser(ptmin, ptmax); + hInPHOJET[pid]->GetYaxis()->SetRangeUser(effmin, effmax); + } + + hInPYTHIA[1]->Add(hInPYTHIA[2]); + hInPYTHIA[1]->Add(hInPYTHIA[3]); + hInPYTHIA[0]->Divide(hInPYTHIA[1], hInPYTHIA[0], 1, 1, "B"); + + TH1D* histPYTHIA = HistInvert(hInPYTHIA[0]); + histPYTHIA->SetLineColor(2); + histPYTHIA->SetMarkerColor(2); + + hInPHOJET[1]->Add(hInPHOJET[2]); + hInPHOJET[1]->Add(hInPHOJET[3]); + hInPHOJET[0]->Divide(hInPHOJET[1], hInPHOJET[0], 1, 1, "B"); + + TH1D* histPHOJET = HistInvert(hInPHOJET[0]); + histPHOJET->SetLineColor(4); + histPHOJET->SetMarkerColor(4); + + TCanvas* cChCorrection = new TCanvas("cChCorrection", "All efficiency", 400, 300); + cChCorrection->Clear(); + + cChCorrection->SetGridy(); + + histPYTHIA->GetYaxis()->SetRangeUser(effmin, effmax); + histPYTHIA->SetTitle("N_{ch}/(#pi + K + p) for primaries at generator level; p_{T} [GeV/c]; Ratio"); + histPYTHIA->Draw(); + + histPHOJET->Draw("SAME"); + + TLegend* legend = new TLegend(0.55, 0.22, 0.79, 0.42); + legend->SetBorderSize(0); + legend->SetFillColor(0); + legend->AddEntry(histPYTHIA, "PYTHIA", "P"); + legend->AddEntry(histPHOJET, "PHOJET", "P"); + legend->Draw(); + + cChCorrection->SaveAs("charged_over_piKp.gif"); + + TFile* file = new TFile("correction.root", "RECREATE"); + histPYTHIA->SetName("histPYTHIA"); + histPHOJET->SetName("histPHOJET"); + histPYTHIA->Write(); + histPHOJET->Write(); + file->Close(); +} + +//_______________________________________________________________________ +TH1D* HistInvert(TH1D* hist) +{ + TH1D* histNew = new TH1D(*hist); + histNew->Reset(); + histNew->SetName(Form("%s_inv", hist->GetName())); + histNew->SetTitle(Form("%s (inverted)", hist->GetTitle())); + + const Int_t nBins = hist->GetNbinsX(); + + for(Int_t i = 1; i <= nBins; i++) { + + if(hist->GetBinContent(i) == 0) + continue; + + histNew->SetBinContent(i, 1.0/hist->GetBinContent(i)); + histNew->SetBinError(i, hist->GetBinError(i)/hist->GetBinContent(i)/hist->GetBinContent(i)); + } + + return histNew; +} diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_phojet.root b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_phojet.root new file mode 100644 index 00000000000..94cf9cb53aa Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_phojet.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_pythia.root b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_pythia.root new file mode 100644 index 00000000000..7ee4a33a13d Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/efficiency/lhc10d_eff_pythia.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_2.76TeV.root b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_2.76TeV.root new file mode 100644 index 00000000000..c93dc0cd9e4 Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_2.76TeV.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_7TeV.root b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_7TeV.root new file mode 100644 index 00000000000..3633cd49069 Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_7TeV.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_900GeV.root b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_900GeV.root new file mode 100644 index 00000000000..465ec8cca66 Binary files /dev/null and b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/charged_spectra/dNdPt_pp_900GeV.root differ diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/convertToMyBins.C b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/convertToMyBins.C new file mode 100644 index 00000000000..9499d88b0d1 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/convertToMyBins.C @@ -0,0 +1,343 @@ +#include +#include +#include +#include + +#include "my_tools.C" + +#include + +using namespace std; + +/* + Info: + + Method to convert charged pT histograms to have the same histogram binning + as our histograms have. They in general go to 100 GeV/c in pT while we only + go to 50. In that range the binning is EXACTLY the same (and the code checks + it). + + In charged spectra there is NO NORMALIZATION uncertainty. We add it by hand + to the syst error. + + * Need at some point to get the final spectra + * Need at some point to get the final normalization uncertainty + + + To run: + + gSystem->AddIncludePath("-I../macros") + gSystem->AddIncludePath("-I../lib") + gROOT->SetMacroPath(".:../macros") + .L my_tools.C+ + .L convertToMyBins.C+ + + + ConvertPPfile("charged_spectra/dNdPt_pp_7TeV.root", "NOPID_pp_7TeV.root", 0.034); + ConvertPPfile("charged_spectra/dNdPt_pp_2.76TeV.root", "NOPID_pp_2.76TeV.root", 0.034); + ConvertPPfile("charged_spectra/dNdPt_pp_900GeV.root", "NOPID_pp_900GeV.root", 0.034); + + + + + Older methods that can be useful for AA: + The method convertToMyBinsAA is for pp and AA spectra. + The method convertToMyBinsRaa is for RAA. + The method convertToMyBinsOld is for AA data when the centralities had to be merged. + + */ + + + +TH1D* Convert(TH1D* histIn, Double_t addSyst=0); +TH1D* ConvertGraph(TGraphErrors* graphIn); +void convertToMyBinsAA(); +void convertToMyBinsRaa(const Char_t* fileName = "PbPb_RAA_2760GeV_200511.root"); +void convertToMyBinsOld(const Char_t* fileName = "dndpt_all_2011-05-15.root"); + + +void ConvertPPfile(const Char_t* inFileName, const Char_t* outFileName, Double_t normSyst=0.0) +{ + TFile* fileDataCharged = FindFileFresh(inFileName); + + if(!fileDataCharged) + return; + + TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat"); + TH1D* hPtCharged = Convert(hPtHelp); + hPtCharged->SetName("hDnDptCharged_PP"); + hPtCharged->SetMarkerStyle(29); + + hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst"); + TH1D* hPtChargedSyst = Convert(hPtHelp, normSyst); + hPtChargedSyst->SetName("hDnDptChargedSyst_PP"); + hPtChargedSyst->SetMarkerStyle(29); + + TFile* fileOut = new TFile(outFileName, "RECREATE"); + hPtCharged->Write(); + hPtChargedSyst->Write(); + fileOut->Close(); +} + +//________________________________________________________________________________________ +TH1D* Convert(TH1D* histIn, Double_t addSyst) +{ + const Int_t nPtBins = 68; + Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, + 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, + 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 , + 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 , + 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0, + 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, + 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 }; + + + TH1D* hist = new TH1D("hist", "", nPtBins, xBins); + hist->SetTitle(histIn->GetTitle()); + hist->GetXaxis()->SetTitle(histIn->GetXaxis()->GetTitle()); + hist->GetYaxis()->SetTitle(histIn->GetYaxis()->GetTitle()); + + const Double_t deltapt = 0.0001; + + for(Int_t bin = 1; bin <= nPtBins; bin++) { + + // check bin size + if(TMath::Abs(hist->GetXaxis()->GetBinLowEdge(bin) - + histIn->GetXaxis()->GetBinLowEdge(bin)) > deltapt) { + cout << "pt edge low does not agree!!!!!!!!!!!!!!!" << endl; + return 0; + } + if(TMath::Abs(hist->GetXaxis()->GetBinUpEdge(bin) - + histIn->GetXaxis()->GetBinUpEdge(bin)) > deltapt) { + cout << "pt edge high does not agree!!!!!!!!!!!!!!!" << endl; + return 0; + } + if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) - + histIn->GetXaxis()->GetBinCenter(bin)) > deltapt) { + cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl; + return 0; + } + + + hist->SetBinContent(bin, histIn->GetBinContent(bin)); + + Double_t error = histIn->GetBinError(bin); + if(addSyst) { + + Double_t extra_error = addSyst*histIn->GetBinContent(bin); + error = TMath::Sqrt(error*error + extra_error*extra_error); + } + hist->SetBinError(bin, error); + } + + return hist; +} + +//________________________________________________________________________________________ +TH1D* ConvertGraph(TGraphErrors* graphIn) +{ + const Int_t nPtBins = 68; + Double_t xBins[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, + 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, + 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 , + 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 , + 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0, + 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0, + 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 }; + + + TH1D* hist = new TH1D("hist", "", nPtBins, xBins); + // hist->SetTitle(graphIn->GetTitle()); + // hist->GetXaxis()->SetTitle(graphIn->GetXaxis()->GetTitle()); + // hist->GetYaxis()->SetTitle(graphIn->GetYaxis()->GetTitle()); + + const Double_t deltapt = 0.0001; + + // Jacek does not fill the first 3 bins of his graph? + for(Int_t bin = 4; bin <= nPtBins; bin++) { + + // check bin size + if(TMath::Abs(hist->GetXaxis()->GetBinCenter(bin) - + graphIn->GetX()[bin-1]) > deltapt) { + cout << "pt center does not agree!!!!!!!!!!!!!!!" << endl; + return 0; + } + + + hist->SetBinContent(bin, graphIn->GetY()[bin-1]); + hist->SetBinError(bin, graphIn->GetEY()[bin-1]); + } + + return hist; +} + +//________________________________________________________________________________________ +void convertToMyBinsAA() +{ + const Int_t nCent = 7; + const Int_t cent[nCent+1] = {-5, 0, 5, 10, 20, 40, 60, 80}; + const Int_t colors[nCent] = {kBlack, kRed+1, kOrange+7, kOrange, kGreen+2, kCyan+2, kBlue+1}; + + TH1D* hPtCharged[nCent] = {0, 0, 0, 0, 0, 0, 0}; + TH1D* hPtChargedSyst[nCent] = {0, 0, 0, 0, 0, 0, 0}; + + for (Int_t i = 0; i < nCent; i++) { + + if(cent[i]==-5) { + + TFile* fileDataCharged = FindFileFresh("dNdPt_pp_2.76TeV.root"); + + if(!fileDataCharged) + return; + + TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat"); + hPtCharged[i] = Convert(hPtHelp); + hPtCharged[i]->SetName("hDnDptCharged_PP"); + hPtCharged[i]->SetMarkerColor(colors[i]); + hPtCharged[i]->SetMarkerStyle(29); + hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst"); + hPtChargedSyst[i] = Convert(hPtHelp, 0.034); // 3.4 % norm error + /* + Dear Peter, + + for 2.76TeV the systematic error for MB->INEL is 3.40%,, we decided to + drop the uncertainty from VdM scan since we don't really use it... + + Michael + */ + hPtChargedSyst[i]->SetName("hDnDptChargedSyst_PP"); + hPtChargedSyst[i]->SetMarkerColor(colors[i]); + hPtChargedSyst[i]->SetMarkerStyle(29); + } else { + + TFile* fileDataCharged = FindFileFresh(Form("dNdPt_PbPb_2.76ATeV_centrality_%d-%d.root", cent[i], cent[i+1])); + + if(!fileDataCharged) + return; + + TH1D* hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_stat"); + hPtCharged[i] = Convert(hPtHelp); + hPtCharged[i]->SetName(Form("hDnDptCharged_%d_%d", cent[i], cent[i+1])); + hPtCharged[i]->SetMarkerColor(colors[i]); + hPtCharged[i]->SetMarkerStyle(29); + hPtHelp = (TH1D*)fileDataCharged->Get("dNdPt_syst"); + hPtChargedSyst[i] = Convert(hPtHelp); + hPtChargedSyst[i]->SetName(Form("hDnDptChargedSyst_%d_%d", cent[i], cent[i+1])); + hPtChargedSyst[i]->SetMarkerColor(colors[i]); + hPtChargedSyst[i]->SetMarkerStyle(29); + } + } + + TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE"); + for (Int_t i = 0; i < nCent; i++) { + hPtCharged[i]->Write(); + hPtChargedSyst[i]->Write(); + } + fileOut->Close(); +} + +//________________________________________________________________________________________ +void convertToMyBinsRaa(const Char_t* fileName) +{ + TFile* fileData = FindFileFresh(fileName); + + if(!fileData) + return; + + TGraphErrors* hRaa_0_5 = (TGraphErrors*)fileData->Get("raa_c0_5_stat"); + TH1D* hRaaCharged_0_5 = ConvertGraph(hRaa_0_5); + hRaaCharged_0_5->SetName("hRaaCharged_0_5"); + + TGraphErrors* hRaa_5_10 = (TGraphErrors*)fileData->Get("raa_c5_10_stat"); + TH1D* hRaaCharged_5_10 = ConvertGraph(hRaa_5_10); + hRaaCharged_5_10->SetName("hRaaCharged_5_10"); + + TGraphErrors* hRaa_10_20 = (TGraphErrors*)fileData->Get("raa_c10_20_stat"); + TH1D* hRaaCharged_10_20 = ConvertGraph(hRaa_10_20); + hRaaCharged_10_20->SetName("hRaaCharged_10_20"); + + TGraphErrors* hRaa_20_40 = (TGraphErrors*)fileData->Get("raa_c20_40_stat"); + TH1D* hRaaCharged_20_40 = ConvertGraph(hRaa_20_40); + hRaaCharged_20_40->SetName("hRaaCharged_20_40"); + + TGraphErrors* hRaa_40_60 = (TGraphErrors*)fileData->Get("raa_c40_60_stat"); + TH1D* hRaaCharged_40_60 = ConvertGraph(hRaa_40_60); + hRaaCharged_40_60->SetName("hRaaCharged_40_60"); + + TGraphErrors* hRaa_60_80 = (TGraphErrors*)fileData->Get("raa_c60_80_stat"); + TH1D* hRaaCharged_60_80 = ConvertGraph(hRaa_60_80); + hRaaCharged_60_80->SetName("hRaaCharged_60_80"); + + TGraphErrors* hRaa_0_20 = (TGraphErrors*)fileData->Get("raa_c0_20_stat"); + TH1D* hRaaCharged_0_20 = ConvertGraph(hRaa_0_20); + hRaaCharged_0_20->SetName("hRaaCharged_0_20"); + + TGraphErrors* hRaa_40_80 = (TGraphErrors*)fileData->Get("raa_c40_80_stat"); + TH1D* hRaaCharged_40_80 = ConvertGraph(hRaa_40_80); + hRaaCharged_40_80->SetName("hRaaCharged_40_80"); + + + TFile* fileOut = new TFile("Raa_Charged.root", "RECREATE"); + hRaaCharged_0_5->Write(); + hRaaCharged_5_10->Write(); + hRaaCharged_10_20->Write(); + hRaaCharged_20_40->Write(); + hRaaCharged_40_60->Write(); + hRaaCharged_60_80->Write(); + hRaaCharged_0_20->Write(); + hRaaCharged_40_80->Write(); + fileOut->Close(); +} + +//________________________________________________________________________________________ + +void convertToMyBinsOld(const Char_t* fileName) +{ + TFile* fileData = FindFileFresh(fileName); + + if(!fileData) + return; + + TH1D* hDnDPt_0_5 = (TH1D*)fileData->Get("dNdPt_PbPb_c0_5"); + TH1D* hDnDptCharged_0_5 = Convert(hDnDPt_0_5); + hDnDptCharged_0_5->SetName("hDnDptCharged_0_5"); + + TH1D* hDnDPt_5_10 = (TH1D*)fileData->Get("dNdPt_PbPb_c5_10"); + TH1D* hDnDptCharged_5_10 = Convert(hDnDPt_5_10); + hDnDptCharged_5_10->SetName("hDnDptCharged_5_10"); + + TH1D* hDnDPt_10_20 = (TH1D*)fileData->Get("dNdPt_PbPb_c10_20"); + TH1D* hDnDptCharged_10_20 = Convert(hDnDPt_10_20); + hDnDptCharged_10_20->SetName("hDnDptCharged_10_20"); + + TH1D* hDnDPt_20_40 = (TH1D*)fileData->Get("dNdPt_PbPb_c20_30"); + TH1D* hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c30_40"); + hDnDPt_20_40->Add(hHelp); + hDnDPt_20_40->Scale(1.0/2.0); + TH1D* hDnDptCharged_20_40 = Convert(hDnDPt_20_40); + hDnDptCharged_20_40->SetName("hDnDptCharged_20_40"); + + TH1D* hDnDPt_40_60 = (TH1D*)fileData->Get("dNdPt_PbPb_c40_50"); + hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c50_60"); + hDnDPt_40_60->Add(hHelp); + hDnDPt_40_60->Scale(1.0/2.0); + TH1D* hDnDptCharged_40_60 = Convert(hDnDPt_40_60); + hDnDptCharged_40_60->SetName("hDnDptCharged_40_60"); + + TH1D* hDnDPt_60_80 = (TH1D*)fileData->Get("dNdPt_PbPb_c60_70"); + hHelp = (TH1D*)fileData->Get("dNdPt_PbPb_c70_80"); + hDnDPt_60_80->Add(hHelp); + hDnDPt_60_80->Scale(1.0/2.0); + TH1D* hDnDptCharged_60_80 = Convert(hDnDPt_60_80); + hDnDptCharged_60_80->SetName("hDnDptCharged_60_80"); + + + TFile* fileOut = new TFile("dNdPt_Charged.root", "RECREATE"); + hDnDptCharged_0_5->Write(); + hDnDptCharged_5_10->Write(); + hDnDptCharged_10_20->Write(); + hDnDptCharged_20_40->Write(); + hDnDptCharged_40_60->Write(); + hDnDptCharged_60_80->Write(); + fileOut->Close(); +} diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/makeSpectra.C b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/makeSpectra.C new file mode 100644 index 00000000000..a9f29cfa811 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/final_spectra/makeSpectra.C @@ -0,0 +1,112 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "my_tools.C" + + +/* + Info: + Multiplies the fraction with the charged spectrum. + Applies a rapidity correction in case of K and p. + + To run: + + gSystem->AddIncludePath("-I../macros") + gSystem->AddIncludePath("-I../lib") + gROOT->SetMacroPath(".:../macros") + .L my_tools.C+ + .L makeSpectra.C+ + + MakeSpectra("../corrected_fraction/fraction_7tev_b_all.root", "NOPID_pp_7TeV.root", "final_spectra_7tev_b.root"); + + + */ + +Double_t rap_correction(Double_t* x, Double_t* par); + +void MakeSpectra(const Char_t* fileName, + const Char_t* fileNameCharged, + const Char_t* outFileName, + const Char_t* endName="PP") +{ + TFile* fileData = FindFileFresh(fileName); + + if(!fileData) + return; + + TFile* fileDataCharged = FindFileFresh(fileNameCharged); + + if(!fileDataCharged) + return; + + const Int_t nPid = 3; + const Char_t* pidNames[nPid] = {"Pion", "Kaon", "Proton"}; + const Char_t* titleNames[nPid] = {"#pi^{+} + #pi^{-}", "K^{+} + K^{-}", "p + #bar{p}"}; + const Double_t mass[nPid] = {0.140, 0.494, 0.938}; + + TH1D* hPtSpectrum[nPid] = {0, 0, 0}; + TH1D* hPtSpectrumSyst[nPid] = {0, 0, 0}; + + TH1D* hPtCharged = (TH1D*)fileDataCharged->Get(Form("hDnDptCharged_%s", endName)); + TH1D* hPtChargedSyst = (TH1D*)fileDataCharged->Get(Form("hDnDptChargedSyst_%s", endName)); + + TF1* fRap = new TF1("fRap", rap_correction, 0.0, 50.0, 2); + + for(Int_t pid = 0; pid < nPid; pid++) { + + // Pions + TH1D* hFraction = (TH1D*)fileData->Get(Form("h%sFraction_%s", pidNames[pid], endName)); + TH1D* hFractionSyst = (TH1D*)fileData->Get(Form("h%sFractionSyst_%s", pidNames[pid], endName)); + + hPtSpectrum[pid] = (TH1D*)hFraction->Clone(Form("h%sSpectrum_%s", pidNames[pid], endName)); + hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid])); + hPtSpectrum[pid]->Multiply(hPtCharged); + + hPtSpectrumSyst[pid] = (TH1D*)hFractionSyst->Clone(Form("h%sSpectrumSyst_%s", pidNames[pid], endName)); + hPtSpectrumSyst[pid]->SetMarkerStyle(1); + hPtSpectrum[pid]->GetYaxis()->SetTitle(Form("1/N_{evt} 1/(2#pi p_{T}) (d^{2}N_{%s})/(dy dp_{T}) (GeV/c)^{-2}", titleNames[pid])); + hPtSpectrumSyst[pid]->Multiply(hPtChargedSyst); + hPtSpectrumSyst[pid]->SetFillStyle(1001); + hPtSpectrumSyst[pid]->SetFillColor(kGray); + + // + // Rapidity correction + // + fRap->SetParameters(0.8, mass[pid]); + hPtSpectrum[pid]->Divide(fRap); + hPtSpectrumSyst[pid]->Divide(fRap); + } + + TFile* fileOut = new TFile(outFileName, "RECREATE"); + for(Int_t pid = 0; pid < nPid; pid++) { + hPtSpectrumSyst[pid]->Write(); + hPtSpectrum[pid]->Write(); + } + fileOut->Close(); +} + +//______________________________________________________________________________ +Double_t rap_correction(Double_t* x, Double_t* par) +{ + Double_t pt = x[0]; + + Double_t eta = par[0]; + Double_t mass = par[1]; + + const Double_t mt = TMath::Sqrt(pt*pt + mass*mass); + + const Double_t rap = TMath::ASinH(pt/mt*TMath::SinH(eta)); + + return rap/eta; +} diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDx.C b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDx.C new file mode 100644 index 00000000000..ae28a2ebe02 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDx.C @@ -0,0 +1,185 @@ +AliAnalysisTask* AddTask(Bool_t AnalysisMC, const Char_t* taskname, Int_t typerun, UInt_t kTriggerInt[], Float_t minc[], Float_t maxc[] ) +{ + // Creates a pid task and adds it to the analysis manager + + // Get the pointer to the existing analysis manager via the static + //access method + //========================================================================= + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTaskHighPtDeDx", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the + // analysis manager The availability of MC handler can also be + // checked here. + // ========================================================================= + if (!mgr->GetInputEventHandler()) { + Error("AddTaskHighPtDeDx", "This task requires an input event handler"); + return NULL; + } + + + + // + // Add track filters, with Golden Cuts + // + AliAnalysisFilter* trackFilterGolden = new AliAnalysisFilter("trackFilter"); + AliESDtrackCuts* esdTrackCutsGolden = + AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); + trackFilterGolden->AddCuts(esdTrackCutsGolden); + + + //old cuts without golden cut + AliAnalysisFilter* trackFilter0 = new AliAnalysisFilter("trackFilter"); + Bool_t clusterCut=0; + Bool_t selPrimaries=kTRUE; + AliESDtrackCuts* esdTrackCutsL0 = new AliESDtrackCuts; + // TPC + if(clusterCut == 0) esdTrackCutsL0->SetMinNClustersTPC(70); + else if (clusterCut == 1) { + esdTrackCutsL0->SetMinNCrossedRowsTPC(70); + esdTrackCutsL0->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); + } + else { + AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut)); + esdTrackCutsL0->SetMinNClustersTPC(70); + } + esdTrackCutsL0->SetMaxChi2PerClusterTPC(4); + esdTrackCutsL0->SetAcceptKinkDaughters(kFALSE); + esdTrackCutsL0->SetRequireTPCRefit(kTRUE); + // ITS + esdTrackCutsL0->SetRequireITSRefit(kTRUE); + esdTrackCutsL0->SetClusterRequirementITS(AliESDtrackCuts::kSPD, + AliESDtrackCuts::kAny); + if(selPrimaries) { + // 7*(0.0026+0.0050/pt^1.01) + esdTrackCutsL0->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); + } + esdTrackCutsL0->SetMaxDCAToVertexZ(2); + esdTrackCutsL0->SetDCAToVertex2D(kFALSE); + esdTrackCutsL0->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL0->SetMaxChi2PerClusterITS(1e10); + esdTrackCutsL0->SetMaxChi2TPCConstrainedGlobal(1e10); + + + trackFilter0->AddCuts(esdTrackCutsL0); + + + + AliAnalysisFilter* trackFilterTPC = new AliAnalysisFilter("trackFilterTPC"); + AliESDtrackCuts* esdTrackCutsTPC = + AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); + trackFilterTPC->AddCuts(esdTrackCutsTPC); + + + + // Create the task and configure it + //======================================================================== + if(typerun==2){//heavy ion analysis + + + AliAnalysisTaskHighPtDeDx* taskHighPtDeDx[6]; + for( Int_t i=0; i<6; ++i ){ + taskHighPtDeDx[i]=0; + Char_t TaskName[256]={0}; + sprintf(TaskName,"%s_%1.0f_%1.0f",taskname,minc[i],maxc[i]); + + taskHighPtDeDx[i] = new AliAnalysisTaskHighPtDeDx(TaskName); + TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + taskHighPtDeDx[i]->SetAnalysisType(type); + taskHighPtDeDx[i]->SetAnalysisMC(AnalysisMC); + taskHighPtDeDx[i]->SetAnalysisPbPb(kTRUE); + taskHighPtDeDx[i]->SetProduceTPCBranch(kFALSE); + taskHighPtDeDx[i]->SetDebugLevel(0); + taskHighPtDeDx[i]->SetEtaCut(0.8); + taskHighPtDeDx[i]->SetVtxCut(10.0); + taskHighPtDeDx[i]->SetMinPt(0.0); // default 2.0 + taskHighPtDeDx[i]->SetLowPtFraction(0.01); // keep 1% of tracks below min pt + taskHighPtDeDx[i]->SetTreeOption(1); + taskHighPtDeDx[i]->SetTrigger1(kTriggerInt[0]); + taskHighPtDeDx[i]->SetTrigger2(kTriggerInt[1]); + taskHighPtDeDx[i]->SetMinCent(minc[i]); + taskHighPtDeDx[i]->SetMaxCent(maxc[i]); + //Set Filtesr + taskHighPtDeDx[i]->SetTrackFilterGolden(trackFilterGolden); + taskHighPtDeDx[i]->SetTrackFilter(trackFilter0); + taskHighPtDeDx[i]->SetTrackFilterTPC(trackFilterTPC); + + mgr->AddTask(taskHighPtDeDx[i]); + + } + + // Create ONLY the output containers for the data produced by the + // task. Get and connect other common input/output containers via + // the manager as below + //======================================================================= + AliAnalysisDataContainer *cout_hist[6]; + for( Int_t i=0; i<6; ++i ){ + + cout_hist[i]=0; + Char_t outFileName[256]={0}; + sprintf(outFileName,"%s_Tree_%1.0f_%1.0f.root",taskname,minc[i],maxc[i]); + //AliAnalysisDataContainer *cout_hist = 0; + + cout_hist[i] = mgr->CreateContainer(Form("output_%1.0f_%1.0f",minc[i],maxc[i]), TList::Class(), AliAnalysisManager::kOutputContainer, outFileName); + mgr->ConnectInput (taskHighPtDeDx[i], 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(taskHighPtDeDx[i], 1, cout_hist[i]); + + } + + // Return task pointer at the end + return taskHighPtDeDx[0]; + } + if(typerun==3){//pp analysis + + + AliAnalysisTaskHighPtDeDx* taskHighPtDeDx = new AliAnalysisTaskHighPtDeDx("taskHighPtDeDxpp"); + TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + taskHighPtDeDx->SetAnalysisType(type); + taskHighPtDeDx->SetAnalysisMC(AnalysisMC); + taskHighPtDeDx->SetAnalysisPbPb(kFALSE); + taskHighPtDeDx->SetProduceTPCBranch(kFALSE); + taskHighPtDeDx->SetDebugLevel(0); + taskHighPtDeDx->SetEtaCut(0.8); + taskHighPtDeDx->SetVtxCut(10.0); + taskHighPtDeDx->SetMinPt(0.0); // default 2.0 + taskHighPtDeDx->SetLowPtFraction(0.01); // keep 1% of tracks below min pt + taskHighPtDeDx->SetTreeOption(1); + taskHighPtDeDx->SetTrigger1(kTriggerInt[0]); + taskHighPtDeDx->SetTrigger2(kTriggerInt[1]); + //Set Filtesr + taskHighPtDeDx->SetTrackFilterGolden(trackFilterGolden); + taskHighPtDeDx->SetTrackFilter(trackFilter0); + taskHighPtDeDx->SetTrackFilterTPC(trackFilterTPC); + + mgr->AddTask(taskHighPtDeDx); + + // Create ONLY the output containers for the data produced by the + // task. Get and connect other common input/output containers via + // the manager as below + //======================================================================= + + AliAnalysisDataContainer *cout_histdedx; + + cout_histdedx=0; + Char_t outFileName[256]={0}; + sprintf(outFileName,"%s_Tree.root",taskname); + //AliAnalysisDataContainer *cout_hist = 0; + + cout_histdedx = mgr->CreateContainer("outputdedx", TList::Class(), AliAnalysisManager::kOutputContainer, outFileName); + mgr->ConnectInput (taskHighPtDeDx, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(taskHighPtDeDx, 1, cout_histdedx); + + // Return task pointer at the end + return taskHighPtDeDx; + + + + + + } + + +} diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDxV0.C b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDxV0.C new file mode 100644 index 00000000000..14ff6c190b9 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AddTaskHighPtDeDxV0.C @@ -0,0 +1,190 @@ +AliAnalysisTask* AddTask(Bool_t AnalysisMC, const Char_t* taskname, Int_t typerun, UInt_t kTriggerInt[], Float_t minc[], Float_t maxc[] ) +{ + // Creates a pid task and adds it to the analysis manager + + // Get the pointer to the existing analysis manager via the static + //access method + //========================================================================= + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTaskHighPtDeDx", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the + // analysis manager The availability of MC handler can also be + // checked here. + // ========================================================================= + if (!mgr->GetInputEventHandler()) { + Error("AddTaskHighPtDeDx", "This task requires an input event handler"); + return NULL; + } + + // + // Add track filters, with Golden Cuts + // + AliAnalysisFilter* trackFilterGolden = new AliAnalysisFilter("trackFilter"); + AliESDtrackCuts* esdTrackCutsGolden = + AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); + trackFilterGolden->AddCuts(esdTrackCutsGolden); + + + //old cuts without golden cut + AliAnalysisFilter* trackFilter0 = new AliAnalysisFilter("trackFilter"); + Bool_t clusterCut=0; + Bool_t selPrimaries=kTRUE; + AliESDtrackCuts* esdTrackCutsL0 = new AliESDtrackCuts; + // TPC + if(clusterCut == 0) esdTrackCutsL0->SetMinNClustersTPC(70); + else if (clusterCut == 1) { + esdTrackCutsL0->SetMinNCrossedRowsTPC(70); + esdTrackCutsL0->SetMinRatioCrossedRowsOverFindableClustersTPC(0.8); + } + else { + AliWarningClass(Form("Wrong value of the clusterCut parameter (%d), using cut on Nclusters",clusterCut)); + esdTrackCutsL0->SetMinNClustersTPC(70); + } + esdTrackCutsL0->SetMaxChi2PerClusterTPC(4); + esdTrackCutsL0->SetAcceptKinkDaughters(kFALSE); + esdTrackCutsL0->SetRequireTPCRefit(kTRUE); + // ITS + esdTrackCutsL0->SetRequireITSRefit(kTRUE); + esdTrackCutsL0->SetClusterRequirementITS(AliESDtrackCuts::kSPD, + AliESDtrackCuts::kAny); + if(selPrimaries) { + // 7*(0.0026+0.0050/pt^1.01) + esdTrackCutsL0->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01"); + } + esdTrackCutsL0->SetMaxDCAToVertexZ(2); + esdTrackCutsL0->SetDCAToVertex2D(kFALSE); + esdTrackCutsL0->SetRequireSigmaToVertex(kFALSE); + esdTrackCutsL0->SetMaxChi2PerClusterITS(1e10); + esdTrackCutsL0->SetMaxChi2TPCConstrainedGlobal(1e10); + + + trackFilter0->AddCuts(esdTrackCutsL0); + + + + AliAnalysisFilter* trackFilterTPC = new AliAnalysisFilter("trackFilterTPC"); + AliESDtrackCuts* esdTrackCutsTPC = + AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); + trackFilterTPC->AddCuts(esdTrackCutsTPC); + + + // Create the task and configure it + //======================================================================== + if(typerun==2){//heavy ion analysis + + AliAnalysisTaskHighPtDeDxV0* taskHighPtDeDxV0[6]; + for( Int_t i=0; i<6; ++i ){ + taskHighPtDeDxV0[i]=0; + Char_t TaskName[256]={0}; + sprintf(TaskName,"%s_%1.0f_%1.0f",taskname,minc[i],maxc[i]); + + taskHighPtDeDxV0[i] = new AliAnalysisTaskHighPtDeDxV0(TaskName); + TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + + taskHighPtDeDxV0[i]->SetAnalysisType(type); + taskHighPtDeDxV0[i]->SetAnalysisMC(AnalysisMC); + taskHighPtDeDxV0[i]->SetAnalysisPbPb(kTRUE); + taskHighPtDeDxV0[i]->SetProduceTPCBranch(kFALSE); + taskHighPtDeDxV0[i]->SetDebugLevel(0); + taskHighPtDeDxV0[i]->SetEtaCut(0.8); + taskHighPtDeDxV0[i]->SetVtxCut(10.0); + taskHighPtDeDxV0[i]->SetMinPt(0.0); // def: 2.0 + taskHighPtDeDxV0[i]->SetMassCut(0.1); // def: 0.03 + taskHighPtDeDxV0[i]->SetTreeOption(1); + taskHighPtDeDxV0[i]->SetRequireRecV0(kFALSE); // def: kTRUE + taskHighPtDeDxV0[i]->SetStoreMcIn(kTRUE); // def: kFALSE + taskHighPtDeDxV0[i]->SetTrigger1(kTriggerInt[0]); + taskHighPtDeDxV0[i]->SetTrigger2(kTriggerInt[1]); + taskHighPtDeDxV0[i]->SetMinCent(minc[i]); + taskHighPtDeDxV0[i]->SetMaxCent(maxc[i]); + //Set Filtesr + taskHighPtDeDxV0[i]->SetTrackFilterGolden(trackFilterGolden); + taskHighPtDeDxV0[i]->SetTrackFilter(trackFilter0); + taskHighPtDeDxV0[i]->SetTrackFilterTPC(trackFilterTPC); + + mgr->AddTask(taskHighPtDeDxV0[i]); + + + } + + // Create ONLY the output containers for the data produced by the + // task. Get and connect other common input/output containers via + // the manager as below + //======================================================================= + AliAnalysisDataContainer *cout_hist[6]; + for( Int_t i=0; i<6; ++i ){ + + cout_hist[i]=0; + Char_t outFileName[256]={0}; + sprintf(outFileName,"%s_Tree_%1.0f_%1.0f.root",taskname,minc[i],maxc[i]); + //AliAnalysisDataContainer *cout_hist = 0; + + cout_hist[i] = mgr->CreateContainer(Form("output_%1.0f_%1.0f",minc[i],maxc[i]), TList::Class(), AliAnalysisManager::kOutputContainer, outFileName); + mgr->ConnectInput (taskHighPtDeDxV0[i], 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(taskHighPtDeDxV0[i], 1, cout_hist[i]); + + } + + // Return task pointer at the end + return taskHighPtDeDxV0[0]; + } + + if(typerun==3){//pp analysis + + AliAnalysisTaskHighPtDeDxV0* taskHighPtDeDxV0 = new AliAnalysisTaskHighPtDeDxV0("taskHighPtDeDxV0pp"); + TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD" + + taskHighPtDeDxV0->SetAnalysisType(type); + taskHighPtDeDxV0->SetAnalysisMC(AnalysisMC); + taskHighPtDeDxV0->SetAnalysisPbPb(kFALSE); + taskHighPtDeDxV0->SetProduceTPCBranch(kFALSE); + taskHighPtDeDxV0->SetDebugLevel(0); + taskHighPtDeDxV0->SetEtaCut(0.8); + taskHighPtDeDxV0->SetVtxCut(10.0); + taskHighPtDeDxV0->SetMinPt(0.0); // def: 2.0 + taskHighPtDeDxV0->SetMassCut(0.03); // def: 0.03 + taskHighPtDeDxV0->SetTreeOption(1); + taskHighPtDeDxV0->SetRequireRecV0(kTRUE); // def: kTRUE + taskHighPtDeDxV0->SetStoreMcIn(kFALSE); // def: kFALSE + taskHighPtDeDxV0->SetTrigger1(kTriggerInt[0]); + taskHighPtDeDxV0->SetTrigger2(kTriggerInt[1]); + //Set Filtesr + taskHighPtDeDxV0->SetTrackFilterGolden(trackFilterGolden); + taskHighPtDeDxV0->SetTrackFilter(trackFilter0); + taskHighPtDeDxV0->SetTrackFilterTPC(trackFilterTPC); + + mgr->AddTask(taskHighPtDeDxV0); + + + + // Create ONLY the output containers for the data produced by the + // task. Get and connect other common input/output containers via + // the manager as below + //======================================================================= + AliAnalysisDataContainer *cout_histdedxv0; + + cout_histdedxv0=0; + Char_t outFileName[256]={0}; + sprintf(outFileName,"%s_Tree.root",taskname); + //AliAnalysisDataContainer *cout_hist = 0; + + cout_histdedxv0 = mgr->CreateContainer("output", TList::Class(), AliAnalysisManager::kOutputContainer, outFileName); + mgr->ConnectInput (taskHighPtDeDxV0, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(taskHighPtDeDxV0, 1, cout_histdedxv0); + + // Return task pointer at the end + return taskHighPtDeDxV0; + + + } + + + + + + +} diff --git a/PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDx.cxx b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDx.cxx new file mode 100644 index 00000000000..096f51ead54 --- /dev/null +++ b/PWGLF/SPECTRA/IdentifiedHighPt/grid/AliAnalysisTaskHighPtDeDx.cxx @@ -0,0 +1,1705 @@ +#include "AliAnalysisTaskHighPtDeDx.h" + +// ROOT includes +#include +#include +#include +#include +#include +#include + +// AliRoot includes +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include "AliCentrality.h" + +#include +#include +#include + + +// STL includes +#include +using namespace std; + + +// +// Responsible: +// Alexandru Dobrin (Wayne State) +// Peter Christiansen (Lund) +// + +/* +To do: + +Separate the code into two + +*/ + + + + +ClassImp(AliAnalysisTaskHighPtDeDx) + +const Double_t AliAnalysisTaskHighPtDeDx::fgkClight = 2.99792458e-2; + +//_____________________________________________________________________________ +AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(): + AliAnalysisTaskSE(), + fESD(0x0), + fAOD(0x0), + fMC(0x0), + fMCStack(0x0), + fMCArray(0x0), + fTrackFilter(0x0), + fTrackFilterGolden(0x0), + fTrackFilterTPC(0x0), + fAnalysisType("ESD"), + fAnalysisMC(kFALSE), + fAnalysisPbPb(kFALSE), + fTPCBranch(kFALSE), + ftrigBit1(0x0), + ftrigBit2(0x0), + fRandom(0x0), + fEvent(0x0), + fTrackArrayGlobalPar(0x0), + fTrackArrayTPCPar(0x0), + fTrackArrayMC(0x0), + fVZEROArray(0x0), + fVtxCut(10.0), + fEtaCut(0.9), + fMinPt(0.1), + fMinCent(0.0), + fMaxCent(100.0), + fLowPtFraction(0.01), + fTreeOption(0), + fMcProcessType(-999), + fTriggeredEventMB(-999), + fVtxStatus(-999), + fZvtx(-999), + fZvtxMC(-999), + fRun(-999), + fEventId(-999), + fListOfObjects(0), + fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0), + fTree(0x0), + fn1(0), + fn2(0), + fcent(0) + + +{ + // Default constructor (should not be used) + + // fRandom = new TRandom(0); // 0 means random seed +} + +//______________________________________________________________________________ +AliAnalysisTaskHighPtDeDx::AliAnalysisTaskHighPtDeDx(const char *name): + AliAnalysisTaskSE(name), + fESD(0x0), + fAOD(0x0), + fMC(0x0), + fMCStack(0x0), + fMCArray(0x0), + fTrackFilter(0x0), + fTrackFilterGolden(0x0), + fTrackFilterTPC(0x0), + fAnalysisType("ESD"), + fAnalysisMC(kFALSE), + fAnalysisPbPb(kFALSE), + fTPCBranch(kFALSE), + ftrigBit1(0x0), + ftrigBit2(0x0), + fRandom(0x0), + fEvent(0x0), + fTrackArrayGlobalPar(0x0), + fTrackArrayMC(0x0), + fVtxCut(10.0), + fEtaCut(0.9), + fMinPt(0.1), + fLowPtFraction(0.01), + fTreeOption(0), + fMcProcessType(-999), + fTriggeredEventMB(-999), + fVtxStatus(-999), + fZvtx(-999), + fZvtxMC(-999), + fRun(-999), + fEventId(-999), + fListOfObjects(0), + fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0), + fTree(0x0), + fn1(0), + fn2(0), + fcent(0) + +{ + + if(fTree)fTree=0; + // Output slot #1 writes into a TList + DefineOutput(1, TList::Class()); +} + +//_____________________________________________________________________________ +AliAnalysisTaskHighPtDeDx::~AliAnalysisTaskHighPtDeDx() +{ + // Destructor + // histograms are in the output list and deleted when the output + // list is deleted by the TSelector dtor + if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { + delete fListOfObjects; + fListOfObjects = 0; + } + if (fRandom) delete fRandom; + fRandom=0; + + // //for proof running; I cannot create tree do to memory limitations -> work with THnSparse + // if (fListOfObjects && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fOutputList; + + + +} + +//______________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::UserCreateOutputObjects() +{ + // This method is called once per worker node + // Here we define the output: histograms and debug tree if requested + // We also create the random generator here so it might get different seeds... + fRandom = new TRandom(0); // 0 means random seed + + OpenFile(1); + fListOfObjects = new TList(); + fListOfObjects->SetOwner(); + + // + // Histograms + // + fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3); + fListOfObjects->Add(fEvents); + + fn1=new TH1F("fn1","fn1",5001,-1,5000); + fListOfObjects->Add(fn1); + + fn2=new TH1F("fn2","fn2",5001,-1,5000); + fListOfObjects->Add(fn2); + + fcent=new TH1F("fcent","fcent",104,-2,102); + fListOfObjects->Add(fcent); + + fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5); + fListOfObjects->Add(fVtx); + + fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30); + fListOfObjects->Add(fVtxBeforeCuts); + + fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30); + fListOfObjects->Add(fVtxAfterCuts); + + if (fAnalysisMC) { + fVtxMC = new TH1I("fVtxMC","Vtx info - no trigger cut (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5); + fListOfObjects->Add(fVtxMC); + } + + if (fTreeOption) { + + fTree = new TTree("tree","Event data"); + fEvent = new DeDxEvent(); + fTree->Branch("event", &fEvent); + + fTrackArrayGlobalPar = new TClonesArray("DeDxTrack", 1000); + fTree->Bronch("trackGlobalPar", "TClonesArray", &fTrackArrayGlobalPar); + if(fTPCBranch){ + fTrackArrayTPCPar = new TClonesArray("DeDxTrack", 1000); + fTree->Bronch("trackTPCPar", "TClonesArray", &fTrackArrayTPCPar); + } + + fVZEROArray = new TClonesArray("VZEROCell", 1000); + fTree->Bronch("cellVZERO", "TClonesArray", &fVZEROArray); + + if (fAnalysisMC) { + fTrackArrayMC = new TClonesArray("DeDxTrackMC", 1000); + fTree->Bronch("trackMC", "TClonesArray", &fTrackArrayMC); + } + + fTree->SetDirectory(0); + + fListOfObjects->Add(fTree); + + } + + // Post output data. + PostData(1, fListOfObjects); +} + +//______________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::UserExec(Option_t *) +{ + // Main loop + + // + // First we make sure that we have valid input(s)! + // + AliVEvent *event = InputEvent(); + if (!event) { + Error("UserExec", "Could not retrieve event"); + return; + } + + + if (fAnalysisType == "ESD"){ + fESD = dynamic_cast(event); + if(!fESD){ + Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + } else { + fAOD = dynamic_cast(event); + if(!fAOD){ + Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + } + + + + if (fAnalysisMC) { + + if (fAnalysisType == "ESD"){ + fMC = dynamic_cast(MCEvent()); + if(!fMC){ + Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + + fMCStack = fMC->Stack(); + + if(!fMCStack){ + Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + } else { // AOD + + fMC = dynamic_cast(MCEvent()); + if(fMC) + fMC->Dump(); + + fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles"); + if(!fMCArray){ + Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__); + this->Dump(); + return; + } + } + } + + + // Get trigger decision + fTriggeredEventMB = 0; //init + + + if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())) + ->IsEventSelected() & ftrigBit1 ){ + fn1->Fill(1); + fTriggeredEventMB = 1; //event triggered as minimum bias + } + if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())) + ->IsEventSelected() & ftrigBit2 ){ + // From AliVEvent: + // kINT7 = BIT(1), // V0AND trigger, offline V0 selection + fTriggeredEventMB += 2; + fn2->Fill(1); + } + + + // Get process type for MC + fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD + + // real data that are not triggered we skip + if(!fAnalysisMC && !fTriggeredEventMB) + return; + + if (fAnalysisMC) { + + if (fAnalysisType == "ESD"){ + + AliHeader* headerMC = fMC->Header(); + if (headerMC) { + + AliGenEventHeader* genHeader = headerMC->GenEventHeader(); + TArrayF vtxMC(3); // primary vertex MC + vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy + if (genHeader) { + genHeader->PrimaryVertex(vtxMC); + } + fZvtxMC = vtxMC[2]; + + // PYTHIA: + AliGenPythiaEventHeader* pythiaGenHeader = + dynamic_cast(headerMC->GenEventHeader()); + if (pythiaGenHeader) { //works only for pythia + fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType()); + } + // PHOJET: + AliGenDPMjetEventHeader* dpmJetGenHeader = + dynamic_cast(headerMC->GenEventHeader()); + if (dpmJetGenHeader) { + fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType()); + } + } + } else { // AOD + + AliAODMCHeader* mcHeader = dynamic_cast(fAOD->FindListObject("mcHeader")); + if(mcHeader) { + fZvtxMC = mcHeader->GetVtxZ(); + + if(strstr(mcHeader->GetGeneratorName(), "Pythia")) { + fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType()); + } else { + fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType()); + } + } + } + } + + // There are 3 cases + // Vertex: NO - status -1 + // Vertex: YES : outside cut - status 0 + // : inside cut - status 1 + // We have to be careful how we normalize because we probably want to + // normalize to: + // Nevents=(No vertex + outside + inside)/(out + in)*in + + + if (fAnalysisType == "ESD") + fZvtx = GetVertex(fESD); + else // AOD + fZvtx = GetVertex(fAOD); + + fVtxStatus = -999; + + if(fZvtx<-990) { + + fVtxStatus = -1; + if(fTriggeredEventMB) + fVtx->Fill(0); + if(fAnalysisMC) + fVtxMC->Fill(0); + } else { + + if(fTriggeredEventMB) + fVtx->Fill(1); + if(fAnalysisMC) + fVtxMC->Fill(1); + fVtxBeforeCuts->Fill(fZvtx); + fVtxStatus = 0; + if (TMath::Abs(fZvtx) < fVtxCut) { + fVtxAfterCuts->Fill(fZvtx); + fVtxStatus = 1; + } + } + + // store MC event data no matter what + if(fAnalysisMC) { + + if (fAnalysisType == "ESD") { + ProcessMCTruthESD(); + } else { // AOD + ProcessMCTruthAOD(); + } + } + + + + Float_t centrality = -10; + // only analyze triggered events + if(fTriggeredEventMB) { + + if (fAnalysisType == "ESD"){ + if(fAnalysisPbPb){ + AliCentrality *centObject = fESD->GetCentrality(); + centrality = centObject->GetCentralityPercentile("V0M"); + if((centrality>fMaxCent)||(centralityFill(centrality); + AnalyzeESD(fESD); + } else { // AOD + if(fAnalysisPbPb){ + AliCentrality *centObject = fAOD->GetCentrality(); + if(centObject) + centrality = centObject->GetCentralityPercentile("V0M"); + if((centrality>fMaxCent)||(centralityFill(centrality); + AnalyzeAOD(fAOD); + } + } + + if( fTreeOption) { + + fEvent->process = fMcProcessType; + fEvent->trig = fTriggeredEventMB; + fEvent->zvtxMC = fZvtxMC; + fEvent->cent = centrality; + + fTree->Fill(); + fTrackArrayGlobalPar->Clear(); + if(fTPCBranch) + fTrackArrayTPCPar->Clear(); + fVZEROArray->Clear(); + + if (fAnalysisMC) + fTrackArrayMC->Clear(); + } + + // Post output data. + PostData(1, fListOfObjects); +} + +//________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent) +{ + fRun = esdEvent->GetRunNumber(); + fEventId = 0; + if(esdEvent->GetHeader()) + fEventId = GetEventIdAsLong(esdEvent->GetHeader()); + + Short_t isPileup = esdEvent->IsPileupFromSPD(); + + // Int_t event = esdEvent->GetEventNumberInFile(); + UInt_t time = esdEvent->GetTimeStamp(); + // ULong64_t trigger = esdEvent->GetTriggerMask(); + Float_t magf = esdEvent->GetMagneticField(); + + + + + + if(fTriggeredEventMB) {// Only MC case can we have not triggered events + + // accepted event + fEvents->Fill(0); + + + if(fVtxStatus!=1) return; // accepted vertex + Int_t nESDTracks = esdEvent->GetNumberOfTracks(); + + if(nESDTracks<1)return; + cout<<"Multiplicity="<Fill(1); + AliESDVZERO *esdV0 = esdEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad + + for (Int_t iCh=0; iCh<64; ++iCh) { + Float_t multv=esdV0->GetMultiplicity(iCh); + Int_t intexv=iCh; + VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell(); + cellv0->cellindex=intexv; + cellv0->cellmult= multv; + } + + + + } // end if triggered + + if(fTreeOption) { + + fEvent->run = fRun; + fEvent->eventid = fEventId; + fEvent->time = time; + //fEvent->cent = centrality; + fEvent->mag = magf; + fEvent->zvtx = fZvtx; + fEvent->vtxstatus = fVtxStatus; + fEvent->pileup = isPileup; + + } + + + + +} + +//________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent) +{ + fRun = aodEvent->GetRunNumber(); + fEventId = 0; + if(aodEvent->GetHeader()) + fEventId = GetEventIdAsLong(aodEvent->GetHeader()); + + UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp(); + Float_t magf = aodEvent->GetMagneticField(); + + //Int_t trackmult = 0; // no pt cuts + //Int_t nadded = 0; + + Short_t isPileup = aodEvent->IsPileupFromSPD(); + + + + + if(fTriggeredEventMB) {// Only MC case can we have not triggered events + + // accepted event + fEvents->Fill(0); + + if(fVtxStatus!=1) return; // accepted vertex + Int_t nAODTracks = aodEvent->GetNumberOfTracks(); + if(nAODTracks<1)return; + + ProduceArrayTrksAOD( aodEvent, kGlobalTrk ); + if(fTPCBranch) + ProduceArrayTrksAOD( aodEvent, kTPCTrk ); + fEvents->Fill(1); + + AliAODVZERO *esdV0 = aodEvent->GetVZEROData();// loop sobre canales del V0 para obtener las multiplicidad + + for (Int_t iCh=0; iCh<64; ++iCh) { + Float_t multv=esdV0->GetMultiplicity(iCh); + Int_t intexv=iCh; + VZEROCell* cellv0 = new((*fVZEROArray)[iCh]) VZEROCell(); + cellv0->cellindex=intexv; + cellv0->cellmult= multv; + } + + + + + } // end if triggered + + if(fTreeOption) { + + //Sort(fTrackArrayGlobalPar, kFALSE); + + fEvent->run = fRun; + fEvent->eventid = fEventId; + fEvent->time = time; + //fEvent->cent = centrality; + fEvent->mag = magf; + fEvent->zvtx = fZvtx; + fEvent->vtxstatus = fVtxStatus; + //fEvent->trackmult = trackmult; + //fEvent->n = nadded; + fEvent->pileup = isPileup; + } +} + +//_____________________________________________________________________________ +Float_t AliAnalysisTaskHighPtDeDx::GetVertex(const AliVEvent* event) const +{ + Float_t zvtx = -999; + + const AliVVertex* primaryVertex = event->GetPrimaryVertex(); + + if(primaryVertex->GetNContributors()>0) + zvtx = primaryVertex->GetZ(); + + return zvtx; +} + +//_____________________________________________________________________________ +Short_t AliAnalysisTaskHighPtDeDx::GetPidCode(Int_t pdgCode) const +{ + // return our internal code for pions, kaons, and protons + + Short_t pidCode = 6; + + switch (TMath::Abs(pdgCode)) { + case 211: + pidCode = 1; // pion + break; + case 321: + pidCode = 2; // kaon + break; + case 2212: + pidCode = 3; // proton + break; + case 11: + pidCode = 4; // electron + break; + case 13: + pidCode = 5; // muon + break; + default: + pidCode = 6; // something else? + }; + + return pidCode; +} + + +//_____________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::ProcessMCTruthESD() +{ + // Fill the special MC histogram with the MC truth info + + Short_t trackmult = 0; + Short_t nadded = 0; + const Int_t nTracksMC = fMCStack->GetNtrack(); + + for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) { + + //Cuts + if(!(fMCStack->IsPhysicalPrimary(iTracks))) + continue; + + TParticle* trackMC = fMCStack->Particle(iTracks); + + Double_t chargeMC = trackMC->GetPDG()->Charge(); + if (chargeMC == 0) + continue; + + if (TMath::Abs(trackMC->Eta()) > fEtaCut ) + continue; + + trackmult++; + + // "charge:pt:p:eta:phi:pidCode" + Float_t ptMC = trackMC->Pt(); + Float_t pMC = trackMC->P(); + Float_t etaMC = trackMC->Eta(); + Float_t phiMC = trackMC->Phi(); + + Int_t pdgCode = trackMC->GetPdgCode(); + Short_t pidCodeMC = 0; + pidCodeMC = GetPidCode(pdgCode); + + // Here we want to add some of the MC histograms! + + // And therefore we first cut here! + if (trackMC->Pt() < fMinPt) { + + // Keep small fraction of low pT tracks + if(fRandom->Rndm() > fLowPtFraction) + continue; + } // else { + // Here we want to add the high pt part of the MC histograms! + // } + + if(fTreeOption) { + + DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC(); + nadded++; + + track->pMC = pMC; + track->ptMC = ptMC; + track->etaMC = etaMC; + track->phiMC = phiMC; + track->qMC = Short_t(chargeMC); + track->pidMC = pidCodeMC; + track->pdgMC = pdgCode; + } + + }//MC track loop + + if(fTreeOption) { + + Sort(fTrackArrayMC, kTRUE); + + fEvent->trackmultMC = trackmult; + fEvent->nMC = nadded; + } + +} + +//_____________________________________________________________________________ +void AliAnalysisTaskHighPtDeDx::ProcessMCTruthAOD() +{ + // Fill the special MC histogram with the MC truth info + + Short_t trackmult = 0; + Short_t nadded = 0; + const Int_t nTracksMC = fMCArray->GetEntriesFast(); + + for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) { + + AliAODMCParticle* trackMC = dynamic_cast(fMCArray->At(iTracks)); + + //Cuts + if(!(trackMC->IsPhysicalPrimary())) + continue; + + + Double_t chargeMC = trackMC->Charge(); + if (chargeMC == 0) + continue; + + if (TMath::Abs(trackMC->Eta()) > fEtaCut ) + continue; + + trackmult++; + + // "charge:pt:p:eta:phi:pidCode" + Float_t ptMC = trackMC->Pt(); + Float_t pMC = trackMC->P(); + Float_t etaMC = trackMC->Eta(); + Float_t phiMC = trackMC->Phi(); + + Int_t pdgCode = trackMC->PdgCode(); + Short_t pidCodeMC = 0; + pidCodeMC = GetPidCode(pdgCode); + + // Here we want to add some of the MC histograms! + + // And therefore we first cut here! + if (trackMC->Pt() < fMinPt) { + + // Keep small fraction of low pT tracks + if(fRandom->Rndm() > fLowPtFraction) + continue; + } // else { + // Here we want to add the high pt part of the MC histograms! + // } + + if(fTreeOption) { + + DeDxTrackMC* track = new((*fTrackArrayMC)[nadded]) DeDxTrackMC(); + nadded++; + + track->pMC = pMC; + track->ptMC = ptMC; + track->etaMC = etaMC; + track->phiMC = phiMC; + track->qMC = Short_t(chargeMC); + track->pidMC = pidCodeMC; + track->pdgMC = pdgCode; + } + + }//MC track loop + + if(fTreeOption) { + + Sort(fTrackArrayMC, kTRUE); + + fEvent->trackmultMC = trackmult; + fEvent->nMC = nadded; + } + +} + +//_____________________________________________________________________________ +Short_t AliAnalysisTaskHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) { + // + // Get the process type of the event. PYTHIA + // + // source PWG0 dNdpt + + Short_t globalType = -1; //init + + if(pythiaType==92||pythiaType==93){ + globalType = 2; //single diffractive + } + else if(pythiaType==94){ + globalType = 3; //double diffractive + } + //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB?? + else { + globalType = 1; //non diffractive + } + return globalType; +} + +//_____________________________________________________________________________ +Short_t AliAnalysisTaskHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) { + // + // get the process type of the event. PHOJET + // + //source PWG0 dNdpt + // can only read pythia headers, either directly or from cocktalil header + Short_t globalType = -1; + + if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction + globalType = 1; + } + else if (dpmJetType==5 || dpmJetType==6) { + globalType = 2; + } + else if (dpmJetType==7) { + globalType = 3; + } + return globalType; +} + +//_____________________________________________________________________________ +ULong64_t AliAnalysisTaskHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const +{ + // To have a unique id for each event in a run! + // Modified from AliRawReader.h + return ((ULong64_t)header->GetBunchCrossNumber()+ + (ULong64_t)header->GetOrbitNumber()*3564+ + (ULong64_t)header->GetPeriodNumber()*16777215*3564); +} + + +//____________________________________________________________________ +TParticle* AliAnalysisTaskHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label) +{ + // + // Finds the first mother among the primary particles of the particle identified by