PWGLF/SPECTRA/PiKaPr/ITSsa/MakeCorrectedITSsaSpectraMultiBin.C

   1 /////////////////////////////////////////////////////////
   2 //Macro to read the results of the 3gausfit for spectra
   3 //E. Biolcati, 13-feb-2010
   4 /////////////////////////////////////////////////////////
   5 #if !defined(__CINT__) || defined(__MAKECINT__)
   6 #include <Riostream.h>
   7 #include <TLatex.h>
   8 #include <TH1F.h>
   9 #include <TH2D.h>
  10 #include <TF1.h>
  11 #include <TMath.h>
  12 #include <TLegend.h>
  13 #include <TLegendEntry.h>
  14 #include <TGraphErrors.h>
  15 #include <TGraphAsymmErrors.h>
  16 #include <TTree.h>
  17 #include <TCanvas.h>
  18 #include <TSystem.h>
  19 #include <TFile.h>
  20 #include <TStyle.h>
  21 #include <TPad.h>
  22 #endif
  23
  24 void Labella(Int_t i);
  25 void Legenda(TH1 *h2, TH1 *h3, TH1 *h4);
  26
  27 //_______________________________________________________
  28 void MakeCorrectedITSsaSpectraMultiBin(Int_t multibin=0){
  29
  30         TString dirNameSIM, dirNameDATA;
  31         TFile *fiSIM  = new TFile(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d/SpectraReco.root",multibin));
  32         TFile *fiDATA = new TFile(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d/SpectraReco.root",multibin));
  33         TFile *fout   = new TFile(Form("ITSsaSpectraCorr_MultiBin%d.root",multibin),"recreate");
  34
  35         //dca correction
  36         TFile *fPanosPosP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_PosP_expostrange_MultBin%d.root",multibin));
  37         TFile *fPanosNegP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_NegP_expostrange_MultBin%d.root",multibin));
  38         TH1F *hPanosPosP= (TH1F*)fPanosPosP->Get("fHistSecTOTCorrDATAMC");
  39         TH1F *hPanosNegP= (TH1F*)fPanosNegP->Get("fHistSecTOTCorrDATAMC");
  40         for(Int_t pbin=0;pbin<=hPanosPosP->GetNbinsX();pbin++)hPanosPosP->SetBinError(pbin,0);
  41         for(Int_t pbin=0;pbin<=hPanosNegP->GetNbinsX();pbin++)hPanosNegP->SetBinError(pbin,0);
  42
  43         //nevts
  44         TH1F *hstat = (TH1F*)fiDATA->Get("fHistNEvents");
  45         Double_t NEvts = hstat->GetBinContent(hstat->FindBin(1.));
  46         cout<<"Event number used for the normalization "<<NEvts<<endl;
  47
  48         //canvas
  49         gStyle->SetOptStat(0);
  50         TCanvas *cs=new TCanvas("cs","cs",1200,700);
  51         cs->Divide(2,1);
  52         TCanvas *cmix=new TCanvas("cmix","cmix",800,900);
  53         cmix->Divide(1,3,0,0);
  54
  55         TH1F *hPieff[2];
  56         TH1F *hKeff[2];
  57         TH1F *hPeff[2];
  58         TH1F *hPiDATA[2];
  59         TH1F *hKDATA[2];
  60         TH1F *hPDATA[2];
  61         TH1F *hPiDATANorm[2];
  62         TH1F *hKDATANorm[2];
  63         TH1F *hPDATANorm[2];
  64         TH1F *hKsuPi[2];
  65         TH1F *hPsuPi[2];
  66         TH1F *hPsuK[2];
  67
  68         hPieff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg0");
  69         hKeff[0]  = (TH1F*)fiSIM->Get("hCorrFacNeg1");
  70         hPeff[0]  = (TH1F*)fiSIM->Get("hCorrFacNeg2");
  71         hPieff[1] = (TH1F*)fiSIM->Get("hCorrFacPos0");
  72         hKeff[1]  = (TH1F*)fiSIM->Get("hCorrFacPos1");
  73         hPeff[1]  = (TH1F*)fiSIM->Get("hCorrFacPos2");
  74         TCanvas *ceff=new TCanvas("ceff","ceff",800,500);
  75         ceff->Divide(2,1);
  76         ceff->cd(1);
  77         gPad->SetGridy();
  78         hPieff[0]->Draw();
  79         hKeff[0]->Draw("same");
  80         hPeff[0]->Draw("same");
  81         ceff->cd(2);
  82         gPad->SetGridy();
  83         hPieff[1]->Draw();
  84         hKeff[1]->Draw("same");
  85         hPeff[1]->Draw("same");
  86
  87         hPiDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg0");
  88         hKDATA[0]  = (TH1F*)fiDATA->Get("hSpectraNeg1");
  89         hPDATA[0]  = (TH1F*)fiDATA->Get("hSpectraNeg2");
  90         hPiDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos0");
  91         hKDATA[1]  = (TH1F*)fiDATA->Get("hSpectraPos1");
  92         hPDATA[1]  = (TH1F*)fiDATA->Get("hSpectraPos2");
  93
  94         hPiDATANorm[0] = new TH1F(*hPiDATA[0]);
  95         hPiDATANorm[0]->SetName("hSpectraPiPlusN");
  96         hPiDATANorm[1] = new TH1F(*hPiDATA[1]);
  97         hPiDATANorm[1]->SetName("hSpectraPiMinusN");
  98         hKDATANorm[0]  = new TH1F(*hKDATA[0]);
  99         hKDATANorm[0]->SetName("hSpectraKPlusN");
 100         hKDATANorm[1]  = new TH1F(*hKDATA[1]);
 101         hKDATANorm[1]->SetName("hSpectraKMinusN");
 102         hPDATANorm[0]  = new TH1F(*hPDATA[0]);
 103         hPDATANorm[0]->SetName("hSpectraPPlusN");
 104         hPDATANorm[1]  = new TH1F(*hPDATA[1]);
 105         hPDATANorm[1]->SetName("hSpectraPMinusN");
 106
 107         for(Int_t i=0;i<2;i++){ //0==> negative, 1==>positive
 108                 //line colors
 109                 hPiDATA[i] -> SetLineColor(2);
 110                 hKDATA[i]  -> SetLineColor(3);
 111                 hPDATA[i]  -> SetLineColor(4);
 112                 hPiDATA[i] -> SetMarkerStyle(27);
 113                 hKDATA[i]  -> SetMarkerStyle(23);
 114                 hPDATA[i]  -> SetMarkerStyle(24);
 115                 hPiDATA[i] -> SetMarkerColor(2);
 116                 hKDATA[i]  -> SetMarkerColor(3);
 117                 hPDATA[i]  -> SetMarkerColor(4);
 118                 hPiDATA[i] -> SetMarkerStyle(23);
 119                 hKDATA[i]  -> SetMarkerStyle(23);
 120                 hPDATA[i]  -> SetMarkerStyle(23);
 121
 122                 hPiDATANorm[i] -> SetLineColor(2);
 123                 hKDATANorm[i]  -> SetLineColor(3);
 124                 hPDATANorm[i]  -> SetLineColor(4);
 125                 hPiDATANorm[i] -> SetMarkerStyle(27);
 126                 hKDATANorm[i]  -> SetMarkerStyle(23);
 127                 hPDATANorm[i]  -> SetMarkerStyle(24);
 128                 hPiDATANorm[i] -> SetMarkerColor(2);
 129                 hKDATANorm[i]  -> SetMarkerColor(3);
 130                 hPDATANorm[i]  -> SetMarkerColor(4);
 131                 hPiDATANorm[i] -> SetMarkerStyle(23);
 132                 hKDATANorm[i]  -> SetMarkerStyle(23);
 133                 hPDATANorm[i]  -> SetMarkerStyle(23);
 134
 135                 //division for efficiency
 136                 hPiDATA[i] -> Divide(hPieff[i]);
 137                 hKDATA[i]  -> Divide(hKeff[i]);
 138                 hPDATA[i]  -> Divide(hPeff[i]);
 139                 hPiDATANorm[i] -> Divide(hPieff[i]);
 140                 hKDATANorm[i]  -> Divide(hKeff[i]);
 141                 hPDATANorm[i]  -> Divide(hPeff[i]);
 142
 143                 //normalization number of events
 144                 hPiDATANorm[i] -> Scale(1./NEvts);
 145                 hKDATANorm[i]  -> Scale(1./NEvts);
 146                 hPDATANorm[i]  -> Scale(1./NEvts);
 147
 148                 //correction factor based on fit to DCA distr for P and Pibar
 149                 hPDATA[0]->Multiply(hPanosPosP);
 150                 hPDATA[1]->Multiply(hPanosNegP);
 151
 152                 //drawing
 153                 cs->cd(i+1);
 154                 gPad->SetLogy();
 155                 gPad->SetGridy();
 156                 gPad->SetGridx();
 157                 Float_t minim=0.01;
 158                 Float_t maxim=10.;
 159                 hPiDATANorm[i]->Draw("e");
 160                 hPiDATANorm[i]->GetYaxis()->SetRangeUser(minim,maxim);
 161                 hPiDATANorm[i]->GetYaxis()->SetTitle("#frac{1}{N}#frac{d^{2}N}{dp_{t}dy}");
 162                 hPiDATANorm[i]->SetTitle("DATA from ITSsa - corrected");
 163                 hKDATANorm[i]->Draw("esames");
 164                 hPDATANorm[i]->Draw("esames");
 165                 Labella(i);
 166                 Legenda(hPiDATA[i],hKDATA[i],hPDATA[i]);
 167                 cs->Update();
 168         }
 169
 170         //fluka correction for pi
 171         TFile *fGeanFlukaPi= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root"));
 172         TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
 173         TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
 174         for(Int_t binPi=0;binPi<=hPiDATA[0]->GetNbinsX();binPi++){
 175                 Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hPiDATA[0]->GetBinCenter(binPi)));
 176                 Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hPiDATA[1]->GetBinCenter(binPi)));
 177                 //cout<<"PiPos  "<<FlukaCorrPiPos<<"  "<<hPiDATA[0]->GetBinCenter(binPi)<<"  " <<binPi <<endl;
 178                 //cout<<"PiNeg  "<<FlukaCorrPiNeg<<"  "<<hPiDATA[0]->GetBinCenter(binPi)<<"  " <<binPi <<endl;
 179                 hPiDATA[0]->SetBinContent(binPi,hPiDATA[0]->GetBinContent(binPi)*FlukaCorrPiPos);
 180                 hPiDATA[1]->SetBinContent(binPi,hPiDATA[1]->GetBinContent(binPi)*FlukaCorrPiNeg);
 181         }
 182         //fluka correction for pi
 183         TFile *fGeanFlukaK= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root"));
 184         TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
 185         TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
 186         for(Int_t binK=0;binK<=hKDATA[0]->GetNbinsX();binK++){
 187                 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hKDATA[0]->GetBinCenter(binK)));
 188                 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hKDATA[1]->GetBinCenter(binK)));
 189                 //cout<<"KPos :"<<FlukaCorrKPos<<"  "<<hKDATA[0]->GetBinCenter(binK)<<"  " <<binK <<endl;
 190                 //cout<<"KNeg :"<<FlukaCorrKNeg<<"  "<<hKDATA[0]->GetBinCenter(binK)<<"  " <<binK <<endl;
 191                 hKDATA[0]->SetBinContent(binK,hKDATA[0]->GetBinContent(binK)*FlukaCorrKPos);
 192                 hKDATA[1]->SetBinContent(binK,hKDATA[1]->GetBinContent(binK)*FlukaCorrKNeg);
 193         }
 194         //fluka correction for P
 195         //ITS specific file for protons/antiprotons
 196         Int_t kPos=0;
 197         Int_t kNeg=1;
 198         TFile* fITS = new TFile ("RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
 199         TH2D * hCorrFlukaITS[2];
 200         hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
 201         hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
 202         for(Int_t icharge = 0; icharge < 2; icharge++){
 203                 Int_t nbins = hPDATA[0]->GetNbinsX();
 204                 Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
 205                 for(Int_t ibin = 0; ibin < nbins; ibin++){
 206                         Float_t pt = hPDATA[0]->GetBinCenter(ibin);
 207                         Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
 208                         Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 209                         if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 210                         if (pt > maxPtCorrection) pt = maxPtCorrection;
 211                         Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
 212                         if(icharge==0){
 213                                 if (correction != 0) {// If the bin is empty this is a  0
 214                                         hPDATA[0]->SetBinContent(ibin,hPDATA[0]->GetBinContent(ibin)*correction);
 215                                         hPDATA[0]->SetBinError(ibin,hPDATA[0]->GetBinError  (ibin)*correction);
 216                                 }else if (hPDATA[0]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 217                                         cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
 218                                         cout << " Bin content: " << hPDATA[0]->GetBinContent(ibin) << endl;
 219                                 }
 220                         }
 221                         if(icharge==1){
 222                                 if (correction != 0) {// If the bin is empty this is a  0
 223                                         hPDATA[1]->SetBinContent(ibin,hPDATA[1]->GetBinContent(ibin)*correction);
 224                                         hPDATA[1]->SetBinError(ibin,hPDATA[1]->GetBinError  (ibin)*correction);
 225                                 }else if (hPDATA[1]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 226                                         cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
 227                                         cout << " Bin content: " << hPDATA[1]->GetBinContent(ibin) << endl;
 228                                 }
 229                         }
 230                 }
 231         }
 232
 233
 234         //mixed particle ratios
 235         for(Int_t i=0;i<2;i++){
 236                 hKsuPi[i] = (TH1F*)hKDATA[i]->Clone("KsuPi");
 237                 hPsuPi[i] = (TH1F*)hPDATA[i]->Clone("PsuPi");
 238                 hPsuK[i]  = (TH1F*)hPDATA[i]->Clone("PsuK");
 239                 hKsuPi[i] -> Divide(hPiDATA[i]);
 240                 hPsuPi[i] -> Divide(hPiDATA[i]);
 241                 hPsuK[i]  -> Divide(hKDATA[i]);
 242         }
 243
 244         //positive/negative ratios
 245         TH1F *hPiratio = (TH1F*)hPiDATA[1]->Clone("PionsRatio");
 246         TH1F *hKratio  = (TH1F*)hKDATA[1]->Clone("KaonsRatio");
 247         TH1F *  hPratio  = (TH1F*)hPDATA[1]->Clone("ProtonsRatio");
 248         hPiratio -> Divide(hPiDATA[0]);
 249         hKratio  -> Divide(hKDATA[0]);
 250         hPratio  -> Divide(hPDATA[0]);
 251
 252
 253         //drawing positive/negative ratios
 254         TCanvas *cratio=new TCanvas("cratio","",980,600);
 255         cratio->Divide(1,3,0,0);
 256         cratio->SetBottomMargin(0.08);
 257         cratio->cd(1);
 258         gPad->SetGridy();
 259         hPiratio->GetYaxis()->SetTitle("");
 260         hPiratio->GetYaxis()->SetRangeUser(0.7,1.3);
 261         hPiratio->SetTitle("");
 262         hPiratio->Draw("mp");
 263         TLatex *ll1=new TLatex(0.7,0.7,"#pi^{+}/#pi^{-}");
 264         ll1->SetNDC();
 265         ll1->SetTextSize(0.14);
 266         ll1->Draw();
 267         cratio->cd(2);
 268         gPad->SetGridy();
 269         hKratio->SetTitle("");
 270         hKratio->GetYaxis()->SetRangeUser(0.7,1.3);
 271         hKratio->Draw("mp");
 272         TLatex *ll2=new TLatex(0.7,0.7,"K^{+}/K^{-}");
 273         ll2->SetNDC();
 274         ll2->SetTextSize(0.14);
 275         ll2->Draw();
 276         cratio->cd(3);
 277         gPad->SetGridy();
 278         hPratio->SetTitle("");
 279         hPratio->GetYaxis()->SetRangeUser(0.7,1.3);
 280         hPratio->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 281         hPratio->Draw("mp");
 282         TLatex *ll3=new TLatex(0.7,0.7,"p/#bar{p}");
 283         ll3->SetNDC();
 284         ll3->SetTextSize(0.144);
 285         ll3->Draw();
 286
 287         //drawing mixed particle ratios
 288         gStyle->SetOptTitle(0);
 289         cmix->cd(1);
 290         hKsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 291         hKsuPi[0]->Draw();
 292         hKsuPi[1]->Draw("same");
 293         hKsuPi[0]->SetMinimum(0);
 294         hKsuPi[0]->SetMarkerStyle(23);
 295         hKsuPi[1]->SetMarkerStyle(24);
 296         TLegend *legm1=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
 297         legm1->AddEntry(hKsuPi[0],"K^{-}/#pi^{-}","p");
 298         legm1->AddEntry(hKsuPi[1],"K^{+}/#pi^{+}","p");
 299         legm1->SetFillColor(0);
 300         legm1->SetBorderSize(0);
 301         legm1->Draw();
 302         hKsuPi[0]->SetMarkerColor(2);
 303         hKsuPi[1]->SetMarkerColor(4);
 304         hKsuPi[0]->SetLineColor(2);
 305         hKsuPi[1]->SetLineColor(4);
 306
 307         cmix->cd(2);
 308         hPsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 309         hPsuPi[0]->Draw();
 310         hPsuPi[1]->Draw("same");
 311         hPsuPi[0]->SetMinimum(0);
 312         hPsuPi[0]->SetMarkerStyle(23);
 313         hPsuPi[1]->SetMarkerStyle(24);
 314         hPsuPi[0]->SetMarkerColor(2);
 315         hPsuPi[1]->SetMarkerColor(4);
 316         hPsuPi[0]->SetLineColor(2);
 317         hPsuPi[1]->SetLineColor(4);
 318         TLegend *legm2=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
 319         legm2->AddEntry(hPsuPi[0],"#bar{p}/#pi^{-}","p");
 320         legm2->AddEntry(hPsuPi[1],"p/#pi^{+}","p");
 321         legm2->SetFillColor(0);
 322         legm2->SetBorderSize(0);
 323         legm2->Draw();
 324
 325         cmix->cd(3);
 326         hPsuK[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 327         hPsuK[0]->Draw();
 328         hPsuK[1]->Draw("same");
 329         hPsuK[0]->SetMinimum(0);
 330         hPsuK[0]->SetMarkerStyle(23);
 331         hPsuK[1]->SetMarkerStyle(24);
 332         hPsuK[0]->SetMarkerColor(2);
 333         hPsuK[1]->SetMarkerColor(4);
 334         hPsuK[0]->SetLineColor(2);
 335         hPsuK[1]->SetLineColor(4);
 336         TLegend *legm3=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
 337         legm3->AddEntry(hPsuPi[0],"#bar{p}/K^{-}","p");
 338         legm3->AddEntry(hPsuPi[1],"p/K^{+}","p");
 339         legm3->SetFillColor(0);
 340         legm3->SetBorderSize(0);
 341         legm3->Draw();
 342
 343   //save histograms in the root files
 344         fout->cd();
 345         hPiDATA[0]->Write();
 346         hKDATA[0] ->Write();
 347         hPDATA[0] ->Write();
 348         hPiDATA[1]->Write();
 349         hKDATA[1] ->Write();
 350         hPDATA[1] ->Write();
 351         hPiDATANorm[0]->Write();
 352         hKDATANorm[0] ->Write();
 353         hPDATANorm[0] ->Write();
 354         hPiDATANorm[1]->Write();
 355         hKDATANorm[1] ->Write();
 356         hPDATANorm[1] ->Write();
 357         //fout->Close();
 358         return;
 359
 360 }//end of the main
 361
 362
 363 //_______________________________________________________
 364 void Labella(Int_t i){
 365         Char_t txt[50];
 366         if(i==0) sprintf(txt,"negative particles");
 367         else sprintf(txt,"positive particles");
 368         TLatex *ltx=new TLatex(0.4,0.3,txt);
 369         ltx->SetNDC();
 370         ltx->SetTextColor(6);
 371         ltx->SetTextFont(22);
 372         ltx->Draw();
 373         return;
 374 }
 375
 376 //_______________________________________________________
 377 void Legenda(TH1 *h2, TH1 *h3, TH1 *h4){
 378         TLegend *leg=new TLegend(0.51,0.11,0.84,0.25,NULL,"brNDC");
 379         leg->SetFillColor(0);
 380         leg->SetBorderSize(0);
 381         TLegendEntry *entry2=leg->AddEntry(h2,"pions","p");
 382         entry2->SetTextColor(2);
 383         TLegendEntry *entry3=leg->AddEntry(h3,"kaons","p");
 384         entry3->SetTextColor(3);
 385         TLegendEntry *entry4=leg->AddEntry(h4,"protons","p");
 386         entry4->SetTextColor(4);
 387         leg->Draw("same");
 388 }
 389
 390 //EOF
 391