| 624a6bb0 |
1 | #include <TString.h> |
| 2 | #include <TStyle.h> |
| 3 | #include <TSystem.h> |
| 4 | #include <TFile.h> |
| 5 | #include <iostream> |
| 6 | #include "TFile.h" |
| 7 | #include "TH1I.h" |
| 8 | #include "TH2.h" |
| 9 | #include "TCanvas.h" |
| 10 | #include "TLegend.h" |
| 11 | #include "TROOT.h" |
| 12 | |
| 13 | TString Charge[]={"Pos","Neg"}; |
| 14 | TString Sign[]={"Plus","Minus"}; |
| 15 | TString Particle[]={"Pion","Kaon","Proton"}; |
| 16 | Int_t Color[3]={1,2,4}; |
| 17 | Int_t Marker[6]={20,21,22,24,25,26}; |
| 18 | Double_t projl[3]={0,20,80}; |
| 19 | Double_t proju[3]={10,40,90}; |
| 20 | |
| 21 | |
| 22 | void MainAnalysis() { |
| acef4f19 |
23 | //TString fold="5SigmaPIDFilterBit6"; |
| 24 | TString fold="5SigmaPID"; |
| 624a6bb0 |
25 | Int_t ibinToCompare=-1; |
| 26 | |
| 27 | //TString sname="Cent0to100_QVec0.0to100.0"; |
| 28 | |
| 29 | TString sname="Cent0to5_QVec0.0to100.0";ibinToCompare=0; |
| 30 | //TString sname="Cent5to10_QVec0.0to100.0";ibinToCompare=1; |
| 31 | //TString sname="Cent10to20_QVec0.0to100.0";ibinToCompare=2; |
| 32 | //TString sname="Cent20to30_QVec0.0to100.0";ibinToCompare=3; |
| 33 | //TString sname="Cent30to40_QVec0.0to100.0";ibinToCompare=4; |
| 34 | |
| 35 | //TString sname="Cent20to40_QVec0.0to2.0"; |
| 36 | //TString sname="Cent20to40_QVec3.0to100.0"; |
| 37 | |
| 38 | TString dataFile = Form("outputAOD%s/Pt.AOD.1._data_ptcut_%s.root",fold.Data(),sname.Data()); |
| 39 | TString mcFile =Form("outputAOD%s/Pt.AOD.1._MC_%s.root",fold.Data(),sname.Data()); |
| 40 | |
| acef4f19 |
41 | gSystem->Load("libCore.so"); |
| 624a6bb0 |
42 | gSystem->Load("libGeom.so"); |
| 624a6bb0 |
43 | gSystem->Load("libPhysics.so"); |
| acef4f19 |
44 | gSystem->Load("libVMC"); |
| 45 | gSystem->Load("libTree"); |
| 46 | gSystem->Load("libProof"); |
| 47 | gSystem->Load("libMatrix"); |
| 48 | gSystem->Load("libSTEERBase"); |
| 49 | gSystem->Load("libESD"); |
| 50 | gSystem->Load("libAOD"); |
| 624a6bb0 |
51 | gSystem->Load("libANALYSIS"); |
| acef4f19 |
52 | gSystem->Load("libOADB"); |
| 624a6bb0 |
53 | gSystem->Load("libANALYSISalice"); |
| acef4f19 |
54 | gSystem->Load("libTENDER"); |
| 55 | gSystem->Load("libCORRFW"); |
| 56 | //gSystem->Load("libPWG0base"); |
| 57 | gSystem->Load("libMinuit"); |
| 58 | gSystem->Load("libPWGTools"); |
| 59 | gSystem->Load("libPWGLFSPECTRA"); |
| 624a6bb0 |
60 | gSystem->AddIncludePath("-I$ALICE_ROOT/include"); |
| 61 | gStyle->SetPalette(1); |
| 62 | |
| 63 | // Open root MC file and get classes |
| 64 | cout << "Analysis Macro" << endl; |
| 65 | cout << " > Reading MC data" << endl; |
| 66 | TFile *_mc = TFile::Open(mcFile.Data()); |
| 67 | AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos"); |
| 68 | AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts"); |
| 69 | AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts"); |
| 70 | // print info about mc track and Event cuts |
| 71 | cout << " -- Info about MC -- "<< endl; |
| 72 | ecuts_mc->PrintCuts(); |
| 73 | // tcuts_mc->PrintCuts(); |
| 74 | // proceed likewise for data |
| 75 | TFile *_data = TFile::Open(dataFile.Data()); |
| 76 | AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos"); |
| 77 | AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts"); |
| 78 | AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts"); |
| 79 | // print info about track and Event cuts |
| 80 | cout << " -- Info about data -- " << endl; |
| 81 | ecuts_data->PrintCuts(); |
| 82 | tcuts_data->PrintCuts(); |
| 83 | |
| 84 | //dedx in data and MC |
| 85 | TCanvas *cPIDSig=new TCanvas("cPIDSig","cPIDSig",700,500); |
| 86 | cPIDSig->Divide(2,2); |
| 87 | cPIDSig->cd(1); |
| 88 | TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTPC"))->Clone(); |
| 89 | gPad->SetLogz(); |
| 90 | gPad->SetGridy(); |
| 91 | gPad->SetGridx(); |
| 92 | PIDSig_data->DrawClone("colz"); |
| 93 | cPIDSig->cd(2); |
| 94 | TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTPC"))->Clone(); |
| 95 | gPad->SetLogz(); |
| 96 | gPad->SetGridy(); |
| 97 | gPad->SetGridx(); |
| 98 | PIDSig_mc->DrawClone("colz"); |
| 99 | cPIDSig->cd(3); |
| 100 | TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTOF"))->Clone(); |
| 101 | gPad->SetLogz(); |
| 102 | gPad->SetGridy(); |
| 103 | gPad->SetGridx(); |
| 104 | PIDSig_data->DrawClone("colz"); |
| 105 | cPIDSig->cd(4); |
| 106 | TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTOF"))->Clone(); |
| 107 | gPad->SetLogz(); |
| 108 | gPad->SetGridy(); |
| 109 | gPad->SetGridx(); |
| 110 | PIDSig_mc->DrawClone("colz"); |
| 111 | |
| 112 | //nsig in data and MC |
| 113 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 114 | TCanvas *cnsig=new TCanvas(Form("cnsig%s",Particle[ipart].Data()),Form("cnsig%s",Particle[ipart].Data()),700,500); |
| 115 | cnsig->Divide(2,3); |
| 116 | cnsig->cd(1); |
| 117 | TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone(); |
| 118 | nsig_data->SetXTitle("Pt (GeV/c)"); |
| 119 | gPad->SetLogz(); |
| 120 | gPad->SetGridy(); |
| 121 | gPad->SetGridx(); |
| 122 | nsig_data->DrawClone("colz"); |
| 123 | cnsig->cd(2); |
| 124 | TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone(); |
| 125 | nsig_mc->SetXTitle("Pt (GeV/c)"); |
| 126 | gPad->SetLogz(); |
| 127 | gPad->SetGridy(); |
| 128 | gPad->SetGridx(); |
| 129 | nsig_mc->DrawClone("colz"); |
| 130 | cnsig->cd(3); |
| 131 | TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone(); |
| 132 | nsig_data->SetXTitle("Pt (GeV/c)"); |
| 133 | gPad->SetLogz(); |
| 134 | gPad->SetGridy(); |
| 135 | gPad->SetGridx(); |
| 136 | nsig_data->DrawClone("colz"); |
| 137 | cnsig->cd(4); |
| 138 | TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone(); |
| 139 | nsig_mc->SetXTitle("Pt (GeV/c)"); |
| 140 | gPad->SetLogz(); |
| 141 | gPad->SetGridy(); |
| 142 | gPad->SetGridx(); |
| 143 | nsig_mc->DrawClone("colz"); |
| 144 | cnsig->cd(5); |
| 145 | TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone(); |
| 146 | nsig_data->SetXTitle("Pt (GeV/c)"); |
| 147 | gPad->SetLogz(); |
| 148 | gPad->SetGridy(); |
| 149 | gPad->SetGridx(); |
| 150 | nsig_data->DrawClone("colz"); |
| 151 | cnsig->cd(6); |
| 152 | TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone(); |
| 153 | nsig_mc->SetXTitle("Pt (GeV/c)"); |
| 154 | gPad->SetLogz(); |
| 155 | gPad->SetGridy(); |
| 156 | gPad->SetGridx(); |
| 157 | nsig_mc->DrawClone("colz"); |
| 158 | } |
| 159 | //qvec in data and MC |
| 160 | for (Int_t icharge=0;icharge<2;icharge++){ |
| 161 | TCanvas *cqvec=new TCanvas(Form("cqvec%s",Charge[icharge].Data()),Form("cqvec%s",Charge[icharge].Data()),700,500); |
| 162 | cqvec->Divide(2,1); |
| 163 | cqvec->cd(1); |
| 164 | TString hname=Form("histq%s",Charge[icharge].Data()); |
| 165 | TH2F *qvec_data = (TH2F*)((TH2F*)hman_data->GetqVecHistogram(hname.Data()))->Clone(); |
| 166 | gPad->SetLogz(); |
| 167 | gPad->SetGridy(); |
| 168 | gPad->SetGridx(); |
| 169 | qvec_data->DrawClone("colz"); |
| 170 | cqvec->cd(2); |
| 171 | TH2F *qvec_mc = (TH2F*)((TH2F*)hman_mc->GetqVecHistogram(hname.Data()))->Clone(); |
| 172 | gPad->SetLogz(); |
| 173 | gPad->SetGridy(); |
| 174 | gPad->SetGridx(); |
| 175 | qvec_mc->DrawClone("colz"); |
| 176 | |
| 177 | TCanvas *cProjqvec=new TCanvas(Form("cProjqvec%s",Charge[icharge].Data()),Form("cProjqvec%s",Charge[icharge].Data()),700,500); |
| 178 | cProjqvec->Divide(3,2); |
| 179 | for(Int_t iproj=0;iproj<3;iproj++){ |
| 180 | TH1F *proj_data=(TH1F*)qvec_data->ProjectionX("data",qvec_data->GetYaxis()->FindBin(projl[iproj]),qvec_data->GetYaxis()->FindBin(proju[iproj])); |
| 181 | if(proj_data->GetEntries()==0)continue; |
| 182 | proj_data->Scale(1/proj_data->GetEntries()); |
| 183 | proj_data->SetTitle(Form("data q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj])); |
| 184 | TH1F *proj_mc=(TH1F*)qvec_mc->ProjectionX("mc",qvec_mc->GetYaxis()->FindBin(projl[iproj]),qvec_mc->GetYaxis()->FindBin(proju[iproj])); |
| 185 | proj_mc->Scale(1/proj_mc->GetEntries()); |
| 186 | proj_mc->SetTitle(Form("mc q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj])); |
| 187 | proj_mc->SetLineColor(2); |
| 188 | cProjqvec->cd(iproj+1); |
| 189 | gPad->SetGridy(); |
| 190 | gPad->SetGridx(); |
| 191 | proj_data->DrawClone(); |
| 192 | proj_mc->DrawClone("same"); |
| 193 | gPad->BuildLegend(); |
| 194 | proj_data->Divide(proj_mc); |
| 195 | cProjqvec->cd(iproj+4); |
| 196 | proj_data->DrawClone(); |
| 197 | } |
| 198 | } |
| 199 | |
| 200 | //efficiencies |
| 201 | Printf("\n\n-> Calculating MC Correction Factors"); |
| 202 | TH1F *CorrFact[6]; |
| 203 | TCanvas *ceff=new TCanvas("ceff","ceff",700,500); |
| 204 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 205 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 206 | Int_t index=ipart+3*icharge; |
| acef4f19 |
207 | //TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 624a6bb0 |
208 | //TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 209 | //TString hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| acef4f19 |
210 | TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 624a6bb0 |
211 | Printf("Getting %s",hname.Data()); |
| 212 | CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 213 | CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 214 | CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 215 | CorrFact[index]->SetMarkerStyle(Marker[index]); |
| 216 | CorrFact[index]->SetMarkerColor(Color[ipart]); |
| 217 | CorrFact[index]->SetLineColor(Color[ipart]); |
| 218 | hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 219 | Printf("... and divide it by %s",hname.Data()); |
| 220 | CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error |
| 221 | if(index==0)CorrFact[index]->DrawClone(); |
| 222 | else CorrFact[index]->DrawClone("same"); |
| 223 | } |
| 224 | } |
| acef4f19 |
225 | TFile *fESD=new TFile("EffAlex/pionEffPbPb.root"); |
| 226 | TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0"); |
| 227 | hEffESD->DrawClone("same"); |
| 624a6bb0 |
228 | gPad->BuildLegend(); |
| 229 | |
| 230 | //Normalization |
| 231 | printf("\n\n-> Spectra Normalization"); |
| 232 | AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts"); |
| 233 | Double_t events_data = ecuts_data->NumberOfEvents(); |
| 234 | Printf(": accepted events: %.0f",events_data); |
| 235 | |
| 624a6bb0 |
236 | //divide RAW for Correction Factor |
| 237 | Printf("\n\n-> Using MC correction factor to correct RAW spectra"); |
| 238 | TH1F *Spectra[6]; |
| 239 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 240 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 241 | Int_t index=ipart+3*icharge; |
| 242 | TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 243 | printf("Getting %s",hname.Data()); |
| 244 | Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 245 | Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 246 | Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 247 | Spectra[index]->SetMarkerStyle(Marker[index]); |
| 248 | Spectra[index]->SetMarkerColor(Color[ipart]); |
| 249 | Spectra[index]->SetLineColor(Color[ipart]); |
| 250 | Printf("... and divide it by %s",hname.Data()); |
| acef4f19 |
251 | Spectra[index]->Divide(CorrFact[index]);//////////////////////////////////////////////////////////////////////////////////////////FIXME |
| 252 | // if(index!=3)Spectra[index]->Divide(CorrFact[index]); |
| 253 | // else{ |
| 254 | // Spectra[index]=AliPWGHistoTools::MyDivideHistosDifferentBins(Spectra[index],hEffESD); |
| 255 | // } |
| 624a6bb0 |
256 | Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION |
| 257 | } |
| 258 | } |
| 259 | |
| acef4f19 |
260 | |
| 624a6bb0 |
261 | //Geant/Fluka Correction |
| 262 | Printf("\nGF correction for Kaons"); |
| 263 | TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root"; |
| 264 | TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data()); |
| 265 | TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles"); |
| 266 | TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles"); |
| 267 | for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){ |
| 268 | Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK))); |
| 269 | Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK))); |
| 270 | Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos); |
| 271 | Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg); |
| 272 | Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos); |
| 273 | Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg); |
| 274 | } |
| 275 | |
| 276 | //Geant Fluka for P |
| 277 | Printf("\nGF correction for Protons"); |
| 278 | const Int_t kNCharge=2; |
| 279 | Int_t kPos=0; |
| 280 | Int_t kNeg=1; |
| 281 | TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root"); |
| 282 | TH2D * hCorrFluka[kNCharge]; |
| 283 | TH2D * hCorrFluka[2]; |
| 284 | hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons"); |
| 285 | hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons"); |
| 286 | for(Int_t icharge = 0; icharge < kNCharge; icharge++){ |
| 287 | Int_t nbins = Spectra[2]->GetNbinsX(); |
| 288 | Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY(); |
| 289 | for(Int_t ibin = 0; ibin < nbins; ibin++){ |
| 290 | Float_t pt = Spectra[2]->GetBinCenter(ibin); |
| 291 | Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1); |
| 292 | Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1); |
| 293 | if (pt < minPtCorrection) pt = minPtCorrection+0.0001; |
| 294 | if (pt > maxPtCorrection) pt = maxPtCorrection; |
| 295 | Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt)); |
| 296 | //cout<<correction<<" charge "<<icharge<<endl; |
| 297 | if(icharge==0){ |
| 298 | if (correction != 0) {// If the bin is empty this is a 0 |
| 299 | Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction); |
| 300 | Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError (ibin)*correction); |
| 301 | }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user |
| 302 | cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries," << endl; |
| 303 | cout << " Bin content: " << Spectra[2]->GetBinContent(ibin) << endl; |
| 304 | } |
| 305 | } |
| 306 | if(icharge==1){ |
| 307 | if (correction != 0) {// If the bin is empty this is a 0 |
| 308 | Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction); |
| 309 | Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError (ibin)*correction); |
| 310 | }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user |
| 311 | cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries," << endl; |
| 312 | cout << " Bin content: " << Spectra[5]->GetBinContent(ibin) << endl; |
| 313 | } |
| 314 | } |
| 315 | } |
| 316 | } |
| 317 | |
| 318 | |
| 319 | DCACorrection(Spectra,hman_data,hman_mc); |
| 320 | |
| 321 | TH1F *MCTruth[6]; |
| 322 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 323 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 324 | Int_t index=ipart+3*icharge; |
| 325 | hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 326 | MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),0.9,1.1))->Clone(); |
| 327 | MCTruth[index]->Scale(1./events_data,"width");//NORMALIZATION |
| 328 | } |
| 329 | } |
| 330 | |
| 331 | //Contamination |
| 332 | Printf("\n\n-> Contamination from MC"); |
| 333 | TH1F *Cont[6]; |
| 334 | TCanvas *ccont=new TCanvas("ccont","ccont",700,500); |
| 335 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 336 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 337 | Int_t index=ipart+3*icharge; |
| 338 | TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 339 | Printf("Getting %s",hname.Data()); |
| 340 | Cont[index] =(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 341 | Cont[index]->SetName(Form("Cont_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 342 | Cont[index]->SetTitle(Form("RecTrue/RecSigma %s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 343 | Cont[index]->SetMarkerStyle(Marker[index]); |
| 344 | Cont[index]->SetMarkerColor(Color[ipart]); |
| 345 | Cont[index]->SetLineColor(Color[ipart]); |
| 346 | hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 347 | Printf("... and divide it by %s",hname.Data()); |
| 348 | Cont[index]->Divide((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1)); |
| 349 | if(index==0)Cont[index]->DrawClone(); |
| 350 | else Cont[index]->DrawClone("same"); |
| 351 | //Spectra[index]->Multiply(Cont[index]); |
| 352 | } |
| 353 | } |
| 354 | gPad->BuildLegend(); |
| 355 | |
| 356 | |
| 357 | //Drawing Final Spectra |
| 358 | TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500); |
| 359 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 360 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 361 | Int_t index=ipart+3*icharge; |
| 362 | if(index==0)Spectra[index]->DrawClone(); |
| 363 | else Spectra[index]->DrawClone("same"); |
| 364 | //MCTruth[index]->DrawClone("same"); |
| 365 | } |
| 366 | } |
| 367 | gPad->BuildLegend(); |
| 368 | |
| 369 | //saving spectra |
| 370 | Printf("\n\n-> Saving spectra in Out%s_%s.root",sname.Data(),fold.Data()); |
| 371 | TFile * fout=new TFile(Form("results/Out_%s_%s.root",sname.Data(),fold.Data()),"RECREATE"); |
| 372 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 373 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 374 | Int_t index=ipart+3*icharge; |
| 375 | Spectra[index]->Write(); |
| 376 | CorrFact[index]->Write(); |
| 377 | MCTruth[index]->Write(); |
| 378 | } |
| 379 | } |
| 380 | |
| 381 | |
| 382 | //if Bin 0-5% with no cut ratio with combined analysis |
| 624a6bb0 |
383 | if(ibinToCompare!=-1){ |
| 384 | TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500); |
| 385 | CratioComb->Divide(3,2); |
| 386 | TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare), |
| 387 | Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)}; |
| 388 | TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120323.root"); |
| 389 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 390 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 391 | Int_t index=ipart+3*icharge; |
| 392 | Spectra[index]->GetXaxis()->SetRangeUser(0,2.5); |
| 393 | TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone(); |
| 394 | CratioComb->cd(ipart+1); |
| 395 | gPad->SetGridy(); |
| 396 | gPad->SetGridx(); |
| 397 | if(icharge==0)Spectra[index]->DrawClone(); |
| 398 | else Spectra[index]->DrawClone("same"); |
| 399 | hcomb->DrawClone("same"); |
| 400 | Spectra[index]->Divide(hcomb); |
| 401 | Spectra[index]->SetMaximum(1.3); |
| 402 | Spectra[index]->SetMinimum(0.7); |
| 403 | CratioComb->cd(ipart+4); |
| 404 | gPad->SetGridy(); |
| 405 | gPad->SetGridx(); |
| 406 | if(icharge==0)Spectra[index]->DrawClone(); |
| 407 | else Spectra[index]->DrawClone("same"); |
| 408 | } |
| 409 | } |
| 410 | } |
| acef4f19 |
411 | |
| 412 | //comparison with charged hadron |
| 413 | Printf("\n\n-> ChargedHadron comparison"); |
| 414 | TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500); |
| 415 | cAllCh->Divide(1,4); |
| 416 | TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("histPtRec",-1,-1))->Clone(); |
| 417 | hChHad_data->Scale(1./events_data,"width");//NORMALIZATION |
| 418 | //fraction of sec in MC |
| 419 | TH1F *hSecAllMC=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0,0))->Clone(); |
| 420 | TH1F *hAllMC=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0,1))->Clone(); |
| 421 | hSecAllMC->Divide(hAllMC); |
| 422 | cAllCh->cd(1); |
| 423 | hSecAllMC->DrawClone(); |
| 424 | cAllCh->cd(2); |
| 425 | hChHad_data->DrawClone(); |
| 426 | for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){ |
| 427 | Double_t en=hChHad_data->GetBinContent(ibin); |
| 428 | Double_t sec=hSecAllMC->GetBinContent(ibin); |
| 429 | hChHad_data->SetBinContent(ibin,en-(en*sec*0.2)); |
| 430 | //Printf("%f %f %d",en,sec,ibin); |
| 431 | //Printf("%f",hChHad_data->GetBinContent(ibin)); |
| 432 | } |
| 433 | hChHad_data->DrawClone("lhistsame"); |
| 434 | //efficiency for primaries |
| 435 | TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtRec",-1,-1))->Clone(); |
| 436 | hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0.5,1.5))->Clone()); |
| 437 | cAllCh->cd(3); |
| 438 | hEff_mc->DrawClone(); |
| 439 | hChHad_data->Divide(hEff_mc); |
| 440 | cAllCh->cd(4); |
| 441 | hChHad_data->Draw(); |
| 442 | |
| 443 | //Comparison of efficiency with TPCTOF ESD analysis |
| 444 | Printf("\n\n-> Calculating Efficiency to be compared with ESD analysis"); |
| 445 | TH1F *EffTRUEPions; |
| 446 | TH1F *EffSIGMAPions; |
| 447 | TCanvas *cEffESD=new TCanvas("cEffESD","cEffESD",700,500); |
| 448 | cEffESD->Divide(1,2); |
| 449 | Int_t icharge=1; |
| 450 | Int_t ipart=0; |
| 451 | //using MC truth |
| 452 | //TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 453 | TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 454 | Printf("Getting %s",hname.Data()); |
| 455 | EffTRUEPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 456 | EffTRUEPions->SetName(Form("Eff TRUE_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 457 | EffTRUEPions->SetTitle(Form("Eff TRUE %s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 458 | EffTRUEPions->SetMarkerStyle(Marker[icharge]); |
| 459 | EffTRUEPions->SetMarkerColor(Color[ipart]); |
| 460 | EffTRUEPions->SetLineColor(Color[ipart]); |
| 461 | hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 462 | Printf("... and divide it by %s",hname.Data()); |
| 463 | EffTRUEPions->Divide(EffTRUEPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error |
| 464 | //using NSigma |
| 465 | //hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 466 | hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 467 | Printf("Getting %s",hname.Data()); |
| 468 | EffSIGMAPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 469 | EffSIGMAPions->SetName(Form("Eff SIGMA_%s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 470 | EffSIGMAPions->SetTitle(Form("Eff SIGMA %s%s",Particle[ipart].Data(),Sign[icharge].Data())); |
| 471 | EffSIGMAPions->SetMarkerStyle(Marker[icharge]); |
| 472 | EffSIGMAPions->SetMarkerColor(Color[ipart+1]); |
| 473 | EffSIGMAPions->SetLineColor(Color[ipart+1]); |
| 474 | hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()); |
| 475 | Printf("... and divide it by %s",hname.Data()); |
| 476 | EffSIGMAPions->Divide(EffSIGMAPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error |
| 477 | cEffESD->cd(1); |
| 478 | //if(icharge==0)EffTRUEPions->DrawClone(); |
| 479 | //else EffTRUEPions->DrawClone("same"); |
| 480 | EffTRUEPions->DrawClone("lhist"); |
| 481 | EffSIGMAPions->DrawClone("lhistsame"); |
| 482 | hEffESD->Draw("lhistsame"); |
| 483 | gPad->BuildLegend(); |
| 484 | cEffESD->cd(2); |
| 485 | TH1F *hRatioTRUE=AliPWGHistoTools::MyDivideHistosDifferentBins(EffTRUEPions,hEffESD); |
| 486 | hRatioTRUE->Draw("lhist"); |
| 487 | TH1F *hRatioSIGMA=AliPWGHistoTools::MyDivideHistosDifferentBins(EffSIGMAPions,hEffESD); |
| 488 | hRatioSIGMA->Draw("lhistsame"); |
| 489 | |
| 490 | |
| 491 | //Muon over Pion Ratio |
| 492 | TCanvas *cMu=new TCanvas("cMu","cMu"); |
| 493 | TH1F *hMuOverPi[2]; |
| 494 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 495 | TString hname=Form("histPtRecTruePrimaryMuon%s",Sign[icharge].Data()); |
| 496 | hMuOverPi[icharge]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(); |
| 497 | hname=Form("histPtRecTruePrimaryPion%s",Sign[icharge].Data()); |
| 498 | hMuOverPi[icharge]->Divide((TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone()); |
| 499 | if(icharge==0)hMuOverPi[icharge]->DrawClone(); |
| 500 | else hMuOverPi[icharge]->DrawClone("same"); |
| 501 | } |
| 502 | |
| 624a6bb0 |
503 | } |
| 504 | |
| 505 | |
| 506 | |
| 507 | void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc){ |
| 508 | |
| 509 | Double_t FitRange[2]={-3,3}; |
| 510 | Int_t nrebin=20; |
| 511 | Printf("\DCACorr"); |
| 512 | TH1F *hcorrection[2]; |
| 513 | TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500); |
| 514 | ccorrection->Divide(2,1); |
| 515 | TString sample[2]={"data","mc"}; |
| 516 | for(Int_t icharge=0;icharge<2;icharge++){ |
| 517 | for(Int_t ipart=0;ipart<3;ipart++){ |
| 518 | Int_t index=ipart+3*icharge; |
| 519 | for(Int_t isample=0;isample<2;isample++){ |
| 520 | TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500); |
| 521 | hcorrection[isample]=(TH1F*)Spectra[index]->Clone(); |
| 522 | hcorrection[isample]->Reset("all"); |
| 523 | cDCA->Divide(8,4); |
| 524 | for(Int_t ibin_data=6;ibin_data<38;ibin_data++){ |
| 525 | Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data); |
| 526 | Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data); |
| 527 | if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone(); |
| 528 | if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone(); |
| 529 | TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone(); |
| 530 | TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone(); |
| 531 | TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone(); |
| 532 | Double_t minentries=5; |
| 533 | if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue; |
| 534 | hmc1->Rebin(nrebin); |
| 535 | hmc2->Rebin(nrebin); |
| 536 | hmc3->Rebin(nrebin); |
| 537 | hToFit->Rebin(nrebin); |
| 538 | //Data and MC can have different stat |
| 539 | hToFit->Sumw2(); |
| 540 | hmc1->Sumw2(); |
| 541 | hmc2->Sumw2(); |
| 542 | hmc3->Sumw2(); |
| 543 | hToFit->Scale(1./hToFit->GetEntries()); |
| 544 | Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries(); |
| 545 | hmc1->Scale(1./normMC); |
| 546 | hmc2->Scale(1./normMC); |
| 547 | hmc3->Scale(1./normMC); |
| 548 | cDCA->cd(ibin_data-5); |
| 549 | gPad->SetGridy(); |
| 550 | gPad->SetGridx(); |
| 551 | gPad->SetLogy(); |
| 552 | hToFit->DrawClone("lhist"); |
| 553 | hmc3->DrawClone("lhist"); |
| 554 | TObjArray *mc = new TObjArray(3); // MC histograms are put in this array |
| 555 | mc->Add(hmc1); |
| 556 | mc->Add(hmc2); |
| 557 | mc->Add(hmc3); |
| 558 | TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise |
| 559 | fit->Constrain(0,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 560 | fit->Constrain(1,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 561 | fit->Constrain(2,0.0,1.0); // constrain fraction 1 to be between 0 and 1 |
| 562 | fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])); |
| 563 | hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1])); |
| 564 | hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data())); |
| 565 | Int_t status = fit->Fit(); // perform the fit |
| 566 | cout << "fit status: " << status << endl; |
| 567 | if (status == 0) { // check on fit status |
| 568 | TH1F* result = (TH1F*) fit->GetPlot(); |
| 569 | TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0); |
| 570 | TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1); |
| 571 | TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2); |
| 572 | //Drawing section |
| 573 | PrimMCPred->SetLineColor(2); |
| 574 | secStMCPred->SetLineColor(6); |
| 575 | secMCPred->SetLineColor(4); |
| 576 | hToFit->DrawClone("lhist"); |
| 577 | result->SetTitle("Fit result"); |
| 578 | result->DrawClone("lhistsame"); |
| 579 | PrimMCPred->DrawClone("lhistsame"); |
| 580 | secStMCPred->DrawClone("lhistsame"); |
| 581 | secMCPred->DrawClone("lhistsame"); |
| 582 | Double_t v1=0,v2=0,v3=0; |
| 583 | Double_t ev1=0,ev2=0,ev3=0; |
| 584 | fit->GetResult(0,v1,ev1); |
| 585 | fit->GetResult(1,v2,ev2); |
| 586 | fit->GetResult(2,v3,ev3); |
| 587 | Printf("\n\n\n\n\nv1:%f v2:%f v3:%f ev1:%f ev2:%f ev3:%f ",v1,v2,v3,ev1,ev2,ev3); |
| 588 | hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3)); |
| 589 | //hcorrection[isample]->SetBinError(ibin_data,TMath::Sqrt(ev1*ev1+ev2*ev2+ev3*ev3)); |
| 590 | hcorrection[isample]->SetBinError(ibin_data,0); |
| 591 | } |
| 592 | else{ |
| 593 | hcorrection[isample]->SetBinContent(ibin_data,1); |
| 594 | hcorrection[isample]->SetBinError(ibin_data,0); |
| 595 | } |
| 596 | Printf("deleting"); |
| 597 | delete hToFit; |
| 598 | } |
| 599 | |
| 600 | ccorrection->cd(icharge+1); |
| 601 | gPad->SetGridy(); |
| 602 | gPad->SetGridx(); |
| 603 | hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data())); |
| 604 | hcorrection[isample]->SetLineWidth(2); |
| 605 | hcorrection[isample]->SetLineColor(Color[ipart]); |
| 606 | hcorrection[isample]->SetLineStyle(isample+1); |
| 607 | hcorrection[isample]->SetMarkerColor(Color[ipart]); |
| 608 | hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,2.5); |
| 609 | if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist"); |
| 610 | else hcorrection[isample]->DrawClone("lhistsame"); |
| 611 | // smooth the DCA correction |
| 612 | // TF1 *fitfun = new TF1("fitfun","[0]+[1]/x^[2]",0.2,0.8); |
| 613 | // hcorrection[isample]->Fit(fitfun,"WRVMN","",0.35,0.8); |
| 614 | // fitfun->SetLineWidth(1.5); |
| 615 | // fitfun->SetLineColor(ipart+1); |
| 616 | // fitfun->SetLineStyle(ipart+1); |
| 617 | // fitfun->DrawClone("same"); |
| 618 | // for(Int_t ibin=1;ibin<30;ibin++){ |
| 619 | // hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin))); |
| 620 | // } |
| 621 | |
| 622 | } |
| 623 | Spectra[index]->Multiply(hcorrection[0]);//multiply for data |
| acef4f19 |
624 | //Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo |
| 624a6bb0 |
625 | } |
| 626 | } |
| 627 | ccorrection->cd(1); |
| 628 | gPad->BuildLegend(); |
| 629 | ccorrection->cd(2); |
| 630 | gPad->BuildLegend(); |
| 631 | |
| 632 | } |
| 633 | |
| 634 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame