PWGLF/SPECTRA/Nuclei/PtSpectra/macros/MakeFinalResults.C

   1 //////////////////////////////////////////////////////////////
   2 //                                                          //
   3 // This macro produces all final results and QA plots for   //
   4 // the DEUTERON nuclei analysis.                            //
   5 //                                                          //
   6 //                                                          //
   7 //                                                          //
   8 //////////////////////////////////////////////////////////////
   9
  10 #include "TFile.h"
  11 #include "TList.h"
  12 #include "TH1D.h"
  13 #include "TH2D.h"
  14 #include "TH3D.h"
  15 #include "THnSparse.h"
  16 #include "TF1.h"
  17 #include "TCanvas.h"
  18 #include "TFractionFitter.h"
  19 #include "TColor.h"
  20
  21
  22
  23 using namespace std;
  24
  25 //
  26 // global variables
  27 //
  28
  29 Char_t * fileNameMcEnhancedNew = "MC/McEnhancedNew.root";
  30 Char_t * fileNameMcEnhancedOld = "MC/McEnhancedOld2.root";
  31 Char_t * fileNameMcUnEnhanced  = "MC/McUnEnhanced.root";
  32 Char_t * fileNameData          = "data/dataFinal.root";
  33
  34 //
  35 // definition of functions
  36 //
  37 void       MakeEfficiency();
  38 TH1D    *  ExtractEfficiency(TList * listMC, Int_t sign = -1, Int_t centrality = 0, Bool_t isTOF = kTRUE);
  39 TF1     *  MakeEtaRapidityGenCorrection(Bool_t drawQA = kFALSE);
  40 //
  41 void       MakeFinalSpectra();
  42 TH1D    *  MakeRawSpectraTPC(TList * list, Int_t sign = -1, Int_t centralityBin = 0);
  43 TH1D    *  MakeRawSpectra(TList * list, Int_t sign = -1, Int_t centralityBin = 0);
  44 void       NormalizeSpectrum(TH1D * spectrum, Float_t dy, Float_t numberOfEvents);
  45 void       PlotSpectraComparison();
  46 void       CombinationAndSystematics();
  47 TH1D    *  GetSpectrumWithSyst(TH1D * histStatError);
  48 //
  49 void       MakeQAplots();
  50 TCanvas *  PlotTpcQA(Char_t * fileName = "", Bool_t isMC = kTRUE,  Bool_t isMCold = kFALSE);
  51 void       PlotSpectraQA(TList * list, Int_t particleType = 0, Int_t sign = -1);
  52 //
  53 void       MakeMaterialCorrection();
  54 Float_t    GetMaterialCorrection(Float_t ptBin, Int_t centralityBin, Bool_t pureMC = kFALSE);
  55
  56
  57
  58 //
  59 // implementation
  60 //
  61 //_______________________________________________________________________
  62 void CombinationAndSystematics() {
  63   //
  64   // produce the final spectra for the paper
  65   //
  66   Int_t kMaxCentrality = 5;
  67   //
  68   TFile fileOut("output/spectraDeuteronCombined.root", "RECREATE");
  69   fileOut.Close();
  70   //
  71   TFile * fileIn = TFile::Open("output/spectraDeuteron.root");
  72   //
  73   TCanvas * canvDeuteron = new TCanvas("canvDeuteron","deuteron spectra");
  74   //
  75   for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
  76     //
  77     TH1D * tpc = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Pos"));
  78     TH1D * tof = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Pos"));
  79     //
  80     TH1D * deuteron = (TH1D *) tof->Clone();
  81     deuteron->SetNameTitle(Form("deuteron_centrality%i",iCentr), Form("deuteron_centrality%i",iCentr));
  82     deuteron->GetYaxis()->SetTitle("d#it{N}/(d#it{y}d#it{p_{T}})");
  83     deuteron->GetXaxis()->SetTitle("p_{T} (GeV/#it{c})");
  84     //
  85     // merge the two
  86     //
  87     deuteron->SetBinContent(4, tpc->GetBinContent(4));  deuteron->SetBinError(4, tpc->GetBinError(4));
  88     deuteron->SetBinContent(5, tpc->GetBinContent(5));  deuteron->SetBinError(5, tpc->GetBinError(4));
  89     //
  90     // get systematics
  91     //
  92     TH1D * deuteronSyst = GetSpectrumWithSyst(deuteron);
  93     //
  94     // drawing and saving
  95     //
  96     if (iCentr == 0) {
  97       deuteronSyst->DrawCopy("E2");
  98       deuteron->DrawCopy("SAME");
  99     } else {
 100       deuteronSyst->DrawCopy("E2,SAME");
 101       deuteron->DrawCopy("SAME");
 102     }
 103     //
 104
 105     //
 106     TFile fileIn2("output/fitsBW.root");
 107     TF1 * fit = (TF1 * ) fileIn2.Get(Form("fit%i_signNeg", iCentr));
 108     if (fit) {
 109       fit->Draw("SAME");
 110     }
 111     fileIn2.Close();
 112     //
 113     TFile fileOut("output/spectraDeuteronCombined.root", "UPDATE");
 114     deuteron->Write();
 115     deuteronSyst->Write();
 116     if (fit) fit->Write();
 117     fileOut.Close();
 118
 119   }
 120
 121
 122
 123 }
 124
 125 //_______________________________________________________________________
 126 TH1D * GetSpectrumWithSyst(TH1D * histStatError) {
 127   //
 128   // add systematic error to the points
 129   //
 130   TH1D * histSyst = (TH1D *) histStatError->Clone();
 131   histSyst->SetNameTitle(Form("%s_SYST", histStatError->GetName()),
 132                          Form("%s_SYST", histStatError->GetName()));
 133   histSyst->SetFillColor(histStatError->GetMarkerColor());
 134   histSyst->SetFillStyle(0);
 135   histSyst->SetDrawOption("E2");
 136   //
 137   for(Int_t i=0; i <  histSyst->GetXaxis()->GetNbins(); i++) { // begin loop over pt-bins
 138       Float_t syst = 0.;
 139       Float_t pt = histSyst->GetXaxis()->GetBinCenter(i);
 140       //
 141       // (1.) tracking and matching normally 4%, but we add 2% for the pt-correction
 142       //
 143       syst += 0.06*0.06;
 144       //
 145       // (2.) PID -- difference between bin counting and fit in TOF
 146       //
 147       if (pt > 1.4) syst += 0.05*0.05;
 148       //
 149       // (3.) material knock-out -- 20% of correction
 150       //
 151       Float_t corr = 1.5*TMath::Exp(1.27259 - 2.527*pt);
 152       syst += (0.2*corr)*(0.2*corr);
 153       //
 154       // (4.) absorption in material (shifted proton exponential, difference between Eulogio's and geant production)
 155       // ==> for anti-deuterons correspondingly more. ==> see plot by Natasha (backup slide APW)
 156       //
 157       syst += 0.03*0.03;
 158       //
 159       // (5.) matching efficiency
 160       //
 161       //if (pt > 1.) syst += 0.05*0.05;
 162       //
 163       // put it to the spectrum
 164       //
 165       Float_t  sum = TMath::Sqrt(syst)*histSyst->GetBinContent(i);
 166       histSyst->SetBinError(i,sum);
 167     } // end loop over pt-bins
 168
 169   //
 170   return histSyst;
 171
 172
 173 }
 174
 175
 176 //_______________________________________________________________________
 177 void PlotSpectraComparison() {
 178   //
 179   // plot TPC & TOF spectra for comparison
 180   //
 181   Int_t kMaxCentrality = 5;
 182   //
 183   //
 184   TFile * fileIn = TFile::Open("output/spectraDeuteron.root");
 185   //
 186   TCanvas * canvRatio = new TCanvas("canvRatio","canvRatio");
 187   //
 188   for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
 189     //
 190     TCanvas * canvComparison = new TCanvas(Form("canvComparison_%i",iCentr),Form("spectra in centrality %i", iCentr));
 191     //
 192     TH1D * tpc = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Pos"));
 193     TH1D * tof = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Pos"));
 194     //
 195     TH1D * tpcNeg = (TH1D *) fileIn->Get(Form("cent%i_sign%s_TPC",iCentr,"Neg"));
 196     TH1D * tofNeg = (TH1D *) fileIn->Get(Form("cent%i_sign%s",iCentr,"Neg"));
 197     tpcNeg->SetMarkerStyle(21);
 198     tofNeg->SetMarkerStyle(25);
 199     //
 200     tof->DrawCopy();
 201     tpc->DrawCopy("SAME");
 202     tpcNeg->DrawCopy("SAME");
 203     tofNeg->DrawCopy("SAME");
 204     //
 205     TFile fileIn2("output/fitsBW.root");
 206     TF1 * fit = (TF1 * ) fileIn2.Get(Form("fit%i_signNeg", iCentr));
 207     fit->Draw("SAME");
 208     fileIn2.Close();
 209     //
 210     //TCanvas * canvRatio = new TCanvas(Form("canvRatio_%i",iCentr),Form("Ratio in centrality %i", iCentr));
 211     canvRatio->cd();
 212     tofNeg->Divide(tof);
 213     if (iCentr == 0) {
 214       tofNeg->DrawCopy();
 215     } else {
 216       tofNeg->DrawCopy("SAME");
 217     }
 218   }
 219
 220
 221 }
 222
 223
 224 //_______________________________________________________________________
 225 void MakeFinalSpectra() {
 226   //
 227   // make final spectra for all centrality bins
 228   //
 229   Int_t kMaxCentrality = 5;
 230   TFile fileOut("output/spectraDeuteron.root", "RECREATE");
 231   fileOut.Close();
 232   //
 233   // get efficiencies and data
 234   //
 235   TFile * inFileEff = TFile::Open("output/efficiencies.root");
 236   TH1D  * efficiencyTrackingNegNew = (TH1D*) inFileEff->Get("efficiencyTrackingNegNew");
 237   TH1D  * efficiencyTrackingPosNew = (TH1D*) inFileEff->Get("efficiencyTrackingPosNew");
 238   TH1D  * efficiencyTofPosNew      = (TH1D*) inFileEff->Get("efficiencyTofPosNew");
 239   TH1D  * efficiencyTofNegNew      = (TH1D*) inFileEff->Get("efficiencyTofNegNew");
 240   //
 241   TFile inFileData(fileNameData);
 242   TList * list = (TList * ) inFileData.Get("akalweit_Nuclei");
 243   //
 244   // positive particles
 245   //
 246   TCanvas * canvAllPos = new TCanvas("canvAllPos","deuteron spectra");
 247   for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
 248     //
 249     TH1D * rawSpecTpc = MakeRawSpectraTPC(list, +1, iCentr);
 250     TH1D * rawSpecTof = MakeRawSpectra(list, +1, iCentr);
 251     rawSpecTpc->Divide(efficiencyTrackingPosNew);
 252     rawSpecTof->Divide(efficiencyTofPosNew);
 253     //
 254     canvAllPos->cd();
 255     if (iCentr == 0) {
 256       rawSpecTpc->DrawCopy("EP");
 257       rawSpecTof->DrawCopy("EPSAME");
 258     } else {
 259       rawSpecTpc->DrawCopy("EPSAME");
 260       rawSpecTof->DrawCopy("EPSAME");
 261     }
 262     //
 263     // apply material correction
 264     //
 265     TFile matFile("output/materialCorrection.root");
 266     TH1D * materialCorrection = (TH1D *) matFile.Get(Form("matCorr_%i",iCentr));
 267     for(Int_t iBin = 0; iBin < rawSpecTpc->GetXaxis()->GetNbins(); iBin++) {
 268       Double_t pT = rawSpecTpc->GetXaxis()->GetBinCenter(iBin);
 269       if (pT > 0.4 && pT < 2.) {
 270         Float_t primFactor = 1. - materialCorrection->GetBinContent(iBin);
 271         rawSpecTof->SetBinContent(iBin, rawSpecTof->GetBinContent(iBin)*primFactor);
 272         rawSpecTof->SetBinError(iBin, rawSpecTof->GetBinError(iBin)*primFactor);
 273         rawSpecTpc->SetBinContent(iBin, rawSpecTpc->GetBinContent(iBin)*primFactor);
 274         rawSpecTpc->SetBinError(iBin, rawSpecTpc->GetBinError(iBin)*primFactor);
 275       }
 276     }
 277     //
 278     matFile.Close();
 279     //
 280     TFile fileOutput("output/spectraDeuteron.root", "UPDATE");
 281     rawSpecTpc->Write();
 282     rawSpecTof->Write();
 283     fileOutput.Close();
 284   }
 285   canvAllPos->Print("QAplots/posSpectra.pdf");
 286   //
 287   // negative particles
 288   //
 289   TCanvas * canvAllNeg = new TCanvas("canvAllNeg","anti-deuteron spectra");
 290   for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
 291     //
 292     TH1D * rawSpecTpc = MakeRawSpectraTPC(list, -1, iCentr);
 293     TH1D * rawSpecTof = MakeRawSpectra(list, -1, iCentr);
 294     rawSpecTpc->Divide(efficiencyTrackingNegNew);
 295     rawSpecTof->Divide(efficiencyTofNegNew);
 296     //
 297     canvAllPos->cd();
 298     if (iCentr == 0) {
 299       rawSpecTpc->DrawCopy("EP");
 300       rawSpecTof->DrawCopy("EPSAME");
 301     } else {
 302       rawSpecTpc->DrawCopy("EPSAME");
 303       rawSpecTof->DrawCopy("EPSAME");
 304     }
 305     TFile fileOutput("output/spectraDeuteron.root", "UPDATE");
 306     rawSpecTpc->Write();
 307     rawSpecTof->Write();
 308     fileOutput.Close();
 309   }
 310   canvAllNeg->Print("QAplots/negSpectra.pdf");
 311
 312
 313
 314 }
 315
 316
 317 //_______________________________________________________________________
 318 TH1D * MakeRawSpectraTPC(TList * list, Int_t sign, Int_t centralityBin) {
 319   //
 320   // Make raw spectra for given sign and centrality bin
 321   //
 322   // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
 323   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 324   // (2.) pT
 325   // (3.) sign
 326   // (4.) rapidity --> filled 4xW
 327   // (5.)  pull TPC dEx --> filled 4x
 328   // (6.) has valid TOF pid signal
 329   // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
 330   // (8..) dca_xy
 331   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
 332   //
 333   TH1D * histCentr = (TH1D *)  list->FindObject("fHistCentrality");
 334   //
 335   THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
 336   hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
 337   //
 338   // centrality 0: 0-5% / 1: 5-10% / 2: 10-20% / 3: 20-30% / 4: 30-40% / 5: 40-50% / 6: 50-60% / 7: 60-70% / 8: 70-80%
 339   //
 340   Int_t centrality = 0;
 341   if (centralityBin == 0) centrality = 0; //  0-10%
 342   if (centralityBin == 1) centrality = 2; // 10-20%
 343   if (centralityBin == 2) centrality = 3; // 20-40%
 344   if (centralityBin == 3) centrality = 5; // 40-60%
 345   if (centralityBin == 4) centrality = 7; // 60-80%
 346   //
 347   hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
 348   Float_t norm = histCentr->GetBinContent(centrality+2);
 349   if (centralityBin != 1) {
 350     hist->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
 351     norm += histCentr->GetBinContent(centrality+3);
 352   }
 353   //
 354   // apply cuts -- Has to be consistent with efficiency
 355   //
 356   hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
 357   hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
 358   hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
 359   hist->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut
 360   //
 361   TH1D * spec    = hist->Projection(2);
 362   spec->SetNameTitle("spec","SPEC");
 363   TH1D * rawSpec = hist->Projection(2);
 364   rawSpec->Reset();
 365   for(Int_t iBin = 0; iBin < spec->GetXaxis()->GetNbins(); iBin++) {
 366     Double_t pT = spec->GetXaxis()->GetBinCenter(iBin);
 367     if (pT > 0.6 && pT < 1.2) {
 368       rawSpec->SetBinContent(iBin, spec->GetBinContent(iBin));
 369       rawSpec->SetBinError(iBin, spec->GetBinError(iBin));
 370     }
 371   }
 372   //
 373   //
 374   NormalizeSpectrum(rawSpec,1,norm);
 375   //
 376   if (sign < 0) rawSpec->SetNameTitle(Form("cent%i_sign%s_TPC",centralityBin,"Neg"),
 377                                       Form("cent%i_sign%s_TPC",centralityBin,"Neg"));
 378   if (sign > 0) rawSpec->SetNameTitle(Form("cent%i_sign%s_TPC",centralityBin,"Pos"),
 379                                       Form("cent%i_sign%s_TPC",centralityBin,"Pos"));
 380   //
 381   Int_t colorList[9] = {600+4, 600+3, 600+2, 600+1, 600, 600-4, 600-7, 600-9, 600-10};
 382   rawSpec->SetMarkerStyle(20);
 383   rawSpec->SetMarkerSize(1.4);
 384   rawSpec->SetMarkerColor(colorList[centrality]);
 385   rawSpec->SetLineColor(colorList[centrality]);
 386   //
 387   hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut
 388   hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 389   hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 390   hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 391
 392   //
 393   return rawSpec;
 394
 395 }
 396
 397
 398
 399 //_______________________________________________________________________
 400 TH1D * MakeRawSpectra(TList * list, Int_t sign, Int_t centralityBin) {
 401   //
 402   // Make raw spectra for given sign and centrality bin
 403   //
 404   //
 405   // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
 406   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 407   // (2.) pT
 408   // (3.) sign
 409   // (4.) rapidity --> filled 4xW
 410   // (5.)  pull TPC dEx --> filled 4x
 411   // (6.) has valid TOF pid signal
 412   // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
 413   // (8..) dca_xy
 414   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
 415   //
 416   TH1D * histCentr = (TH1D *)  list->FindObject("fHistCentrality");
 417   //
 418   THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
 419   hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
 420   //
 421   // centrality 0: 0-5% / 1: 5-10% / 2: 10-20% / 3: 20-30% / 4: 30-40% / 5: 40-50% / 6: 50-60% / 7: 60-70% / 8: 70-80%
 422   //
 423   Int_t centrality = 0;
 424   if (centralityBin == 0) centrality = 0; //  0-10%
 425   if (centralityBin == 1) centrality = 2; // 10-20%
 426   if (centralityBin == 2) centrality = 3; // 20-40%
 427   if (centralityBin == 3) centrality = 5; // 40-60%
 428   if (centralityBin == 4) centrality = 7; // 60-80%
 429   //
 430   hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
 431   Float_t norm = histCentr->GetBinContent(centrality+2);
 432   if (centralityBin != 1) {
 433     hist->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
 434     norm += histCentr->GetBinContent(centrality+3);
 435   }
 436   //
 437   //
 438   //
 439   hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
 440   hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
 441   //
 442   hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
 443   //
 444   hist->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut
 445   //hist->GetAxis(7)->SetRangeUser(-2.,2.5);
 446   //
 447   // PID in the TOF
 448   //
 449   hist->GetAxis(6)->SetRangeUser(1,1);        // TOF-PID cut: require hasTOF
 450   //
 451   TH1D * rawSpec = hist->Projection(2);
 452   rawSpec->Reset();
 453   if (sign < 0) rawSpec->SetNameTitle(Form("cent%i_sign%s",centralityBin,"Neg"),
 454                                       Form("cent%i_sign%s",centralityBin,"Neg"));
 455   if (sign > 0) rawSpec->SetNameTitle(Form("cent%i_sign%s",centralityBin,"Pos"),
 456                                       Form("cent%i_sign%s",centralityBin,"Pos"));
 457   //
 458   // --> simple PID: hist->GetAxis(7)->SetRangeUser(-0.8.,0.8);  // TOF-PID cut
 459   //TF1 fitFunc("fitFunc","gaus(0) + [3] + [4]*x",-1.5,1.5);
 460   TF1 fTOFsignal("fTOFsignal", "(x <= ([3] + [1])) * [0] * TMath::Gaus(x, [1], [2]) + (x > ([3] + [1])) * [0] * TMath::Gaus([3] + [1], [1], [2]) * TMath::Exp(-([3]) * (x - [3] - [1]) / ([2] * [2]))    +    [4] + [5]*x   +   [6]*TMath::Exp(-[7]*x) ", -2.2,2.2); //0:norm, 1:mean, 2:sigma, 3:tail
 461   //
 462   //
 463   TH2D * m2fit = (TH2D *) hist->Projection(7,2);
 464   m2fit->SetNameTitle(Form("m2fits_%i_%i",sign,centralityBin),
 465                       Form("m**2 fits for sign %i and centrality %i",sign,centralityBin));
 466   hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut
 467   hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 468   hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 469   hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 470   //
 471   // rebin because of statistics
 472   //
 473   m2fit->RebinY(2);
 474   //
 475   TCanvas * canvM2fits = new TCanvas(Form("canvM2fits_%i_%i",sign,centralityBin),
 476                                      Form("m**2 fits for sign %i and centrality %i",sign,centralityBin));
 477   canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf(",sign,centralityBin));
 478   //
 479   Double_t oldParStartValue[8] = {664.78 , 0.0178744 , 0.212441 , 0.229438 , 171.889, 0. , 0.0203292};
 480   //
 481   for(Int_t iBin = 1; iBin < m2fit->GetXaxis()->GetNbins(); iBin++) {
 482     //
 483     Float_t ptBin = m2fit->GetXaxis()->GetBinCenter(iBin);
 484     if (ptBin < 0.6) continue;
 485     //
 486     TH1D * m2distr = m2fit->ProjectionY(Form("m2_%f",ptBin),iBin,iBin);
 487     //    m2distr->RebinX(2);
 488     //
 489     m2distr->SetMarkerColor(kBlue);
 490     m2distr->SetMarkerStyle(20);
 491     m2distr->SetMarkerSize(1.4);
 492     //
 493     // initialize signal with background parameters
 494     //
 495     fTOFsignal.SetParameters(oldParStartValue);
 496     //
 497     // improve fit convergence
 498     //
 499     fTOFsignal.SetParLimits(1,-0.2,0.2);
 500     fTOFsignal.SetParLimits(2,0.05,0.8);
 501     fTOFsignal.SetParLimits(3,0.,1.5);
 502     //
 503     fTOFsignal.SetParameter(4, m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(-1.)));
 504     fTOFsignal.SetParLimits(4, 0.,  m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(0.)));
 505     //
 506     fTOFsignal.SetParLimits(7,0.,10.);
 507     fTOFsignal.SetParLimits(6,0., 5.);
 508     //
 509     if (centralityBin < 3 && ptBin < 2.4) {
 510       fTOFsignal.SetRange(-1.3,2.2);
 511     } else{
 512       fTOFsignal.SetRange(-2.2,2.2);
 513     }
 514     //
 515     Int_t fitStatus = m2distr->Fit(&fTOFsignal, "QNR");
 516     cout << fTOFsignal.GetParameter(3) << endl;
 517     //
 518     for(Int_t ival =0; ival < 7; ival++) oldParStartValue[ival] = fTOFsignal.GetParameter(ival);
 519     //
 520     canvM2fits->cd();
 521     m2distr->GetYaxis()->SetRangeUser(0.1, m2distr->GetBinContent(m2distr->GetXaxis()->FindBin(0.))*1.5);
 522     m2distr->DrawCopy("EP");
 523     fTOFsignal.DrawCopy("same");
 524     //
 525     //
 526     // bin counting at low pt and subtract the background at high pt
 527     //
 528     Float_t yield =1;//  m2distr->Integral(m2distr->GetXaxis()->FindBin(-0.5), m2distr->GetXaxis()->FindBin(+0.5)) -
 529       (m2distr->Integral(m2distr->GetXaxis()->FindBin(-1.), m2distr->GetXaxis()->FindBin(-0.5)) + m2distr->Integral(m2distr->GetXaxis()->FindBin(0.5), m2distr->GetXaxis()->FindBin(1.)));
 530     //
 531     //
 532     //
 533     cout << "pT: " << ptBin
 534          << " , " << fTOFsignal.GetParameter(0)
 535          << " , " << fTOFsignal.GetParameter(1)
 536          << " , " << fTOFsignal.GetParameter(2)
 537          << " , " << fTOFsignal.GetParameter(3)
 538          << " , " << fTOFsignal.GetParameter(4)
 539          << " , " << fTOFsignal.GetParameter(5)
 540          << " , " << fTOFsignal.GetParameter(6)
 541          << " , " << fTOFsignal.GetParameter(7) << endl;
 542     //
 543     //
 544     fTOFsignal.SetParameter(4,0.);
 545     fTOFsignal.SetParameter(5,0.);
 546     fTOFsignal.SetParameter(6,0.);
 547     //
 548     //if (ptBin > 2.3)
 549     yield = fTOFsignal.Integral(-1.,1.5)/m2distr->GetBinWidth(10);
 550     //     yield = m2distr->Integral(m2distr->GetXaxis()->FindBin(-1.), m2distr->GetXaxis()->FindBin(+1.));
 551     //
 552     Float_t yieldErr = TMath::Sqrt(yield);
 553     //    if (fTOFsignal.GetParameter(0) != 0 ) yieldErr = yield*fTOFsignal.GetParError(0)/fTOFsignal.GetParameter(0);
 554     //
 555     Float_t maxPt = 4.2;
 556     if (centralityBin > 2) maxPt = 3.2;
 557     if (ptBin > 0.5 && ptBin < maxPt && yield > 0) {
 558       rawSpec->SetBinContent(iBin, yield);
 559       rawSpec->SetBinError(iBin, yieldErr);
 560     }
 561     //
 562     //
 563     canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf",sign,centralityBin));
 564   }
 565   canvM2fits->Print(Form("QA/M2fits_%i_%i.pdf)",sign,centralityBin));
 566   //
 567   rawSpec->Sumw2();
 568   //
 569   NormalizeSpectrum(rawSpec,1,norm);
 570   //
 571   Int_t colorList[9] = {600+4, 600+3, 600+2, 600+1, 600, 600-4, 600-7, 600-9, 600-10};
 572   rawSpec->SetMarkerStyle(24);
 573   rawSpec->SetMarkerSize(1.4);
 574   rawSpec->SetMarkerColor(colorList[centrality]);
 575   rawSpec->SetLineColor(colorList[centrality]);
 576   //
 577   return rawSpec;
 578
 579 }
 580
 581 //_______________________________________________________________________
 582 void NormalizeSpectrum(TH1D * spectrum, Float_t dy, Float_t numberOfEvents) {
 583   //
 584   // make (1 / Nev)*  dN/(dy*dpt)
 585   //
 586   spectrum->Scale(1./dy);
 587   spectrum->Scale(1./numberOfEvents);
 588   //
 589   Int_t nBins = spectrum->GetNbinsX();
 590   for(Int_t i=0; i < nBins+1; i++) {
 591     Double_t Content =spectrum->GetBinContent(i);
 592     Double_t error =spectrum->GetBinError(i);
 593     if (spectrum->GetBinWidth(i) == 0) continue;
 594     spectrum->SetBinContent(i, Content/spectrum->GetBinWidth(i));
 595     spectrum->SetBinError(i, error/spectrum->GetBinWidth(i));
 596   }
 597
 598 }
 599
 600
 601
 602 //_______________________________________________________________________
 603 void MakeEfficiency() {
 604   //
 605   // make efficiencies for TPC and TOF part.
 606   //
 607   TFile inFileMCold(fileNameMcEnhancedOld);
 608   TList * listMCold = (TList * ) inFileMCold.Get("akalweit_Nuclei");
 609   //
 610   //
 611   // TPC and tracking as well as matching efficiencies
 612   //
 613   //
 614   TH1D * efficiencyTrackingNegOld = ExtractEfficiency(listMCold, -1, 0, kFALSE);
 615   TH1D * efficiencyTrackingPosOld = ExtractEfficiency(listMCold, +1, 0, kFALSE);
 616   efficiencyTrackingNegOld->SetNameTitle("efficiencyTrackingNegOld","efficiencyTrackingNegOld");
 617   efficiencyTrackingPosOld->SetNameTitle("efficiencyTrackingPosOld","efficiencyTrackingPosOld");
 618   //
 619   TH1D * efficiencyTofNegOld = ExtractEfficiency(listMCold, -1, 0, kTRUE);
 620   TH1D * efficiencyTofPosOld = ExtractEfficiency(listMCold, +1, 0, kTRUE);
 621   efficiencyTofNegOld->SetNameTitle("efficiencyTofNegOld","efficiencyTofNegOld");
 622   efficiencyTofPosOld->SetNameTitle("efficiencyTofPosOld","efficiencyTofPosOld");
 623   //
 624   TF1 * etaGenCorrection = MakeEtaRapidityGenCorrection(kFALSE);
 625   efficiencyTrackingPosOld->Divide(etaGenCorrection);
 626   efficiencyTrackingNegOld->Divide(etaGenCorrection);
 627   efficiencyTofPosOld->Divide(etaGenCorrection);
 628   efficiencyTofNegOld->Divide(etaGenCorrection);
 629   //
 630   // new simulation
 631   //
 632   TFile inFileMCnew(fileNameMcEnhancedNew);
 633   TList * listMCnew = (TList * ) inFileMCnew.Get("akalweit_Nuclei");
 634   //
 635   TH1D * efficiencyTrackingNegNew = ExtractEfficiency(listMCnew, -1, 0, kFALSE);
 636   TH1D * efficiencyTrackingPosNew = ExtractEfficiency(listMCnew, +1, 0, kFALSE);
 637   efficiencyTrackingNegNew->SetNameTitle("efficiencyTrackingNegNew","efficiencyTrackingNegNew");
 638   efficiencyTrackingPosNew->SetNameTitle("efficiencyTrackingPosNew","efficiencyTrackingPosNew");
 639   efficiencyTrackingNegNew->SetMarkerStyle(21);
 640   efficiencyTrackingPosNew->SetMarkerStyle(21);
 641   //
 642   TH1D * efficiencyTofNegNew = ExtractEfficiency(listMCnew, -1, 0, kTRUE);
 643   TH1D * efficiencyTofPosNew = ExtractEfficiency(listMCnew, +1, 0, kTRUE);
 644   efficiencyTofNegNew->SetNameTitle("efficiencyTofNegNew","efficiencyTofNegNew");
 645   efficiencyTofPosNew->SetNameTitle("efficiencyTofPosNew","efficiencyTofPosNew");
 646   efficiencyTofNegNew->SetMarkerStyle(24);
 647   efficiencyTofPosNew->SetMarkerStyle(24);
 648   //
 649   TCanvas * canvEfficiencyTracking = new TCanvas("canvEfficiencyTracking","Efficiencies Tracking");
 650   efficiencyTrackingNegOld->GetYaxis()->SetRangeUser(0.,1.3);
 651   efficiencyTrackingNegOld->DrawCopy("E");
 652   efficiencyTrackingPosOld->DrawCopy("ESame");
 653   //
 654   efficiencyTrackingNegNew->DrawCopy("ESame");
 655   efficiencyTrackingPosNew->DrawCopy("ESame");
 656   //
 657   TCanvas * canvEfficiencyTof = new TCanvas("canvEfficiencyTof","Efficiencies Tracking");
 658   efficiencyTofNegNew->GetYaxis()->SetRangeUser(0.,1.3);
 659   efficiencyTofNegNew->DrawCopy("E");
 660   efficiencyTofPosNew->DrawCopy("ESame");
 661   efficiencyTofNegOld->DrawCopy("ESAME");
 662   efficiencyTofPosOld->DrawCopy("ESame");
 663   //
 664   canvEfficiencyTracking->Print("QAplots/efficiencyTracking.pdf");
 665   canvEfficiencyTof->Print("QAplots/efficiencyTof.pdf");
 666   //
 667   // store output
 668   //
 669   TFile * outFile = new TFile("output/efficiencies.root","RECREATE");
 670   efficiencyTrackingNegNew->Write();
 671   efficiencyTrackingPosNew->Write();
 672   efficiencyTofNegNew->Write();
 673   efficiencyTofPosNew->Write();
 674   outFile->Close();
 675   //
 676   //
 677   TCanvas * canvEfficiencyMatching = new TCanvas("canvEfficiencyMatching","matching efficiency");
 678   efficiencyTofNegNew->Divide(efficiencyTrackingNegNew);
 679   efficiencyTofPosNew->Divide(efficiencyTrackingPosNew);
 680   efficiencyTofPosNew->DrawCopy();
 681   efficiencyTofNegNew->DrawCopy("SAME");
 682   //
 683   TCanvas * canvEfficiencyRatio = new TCanvas("canvEfficiencyRatio","Efficiencies ratio");
 684   efficiencyTrackingPosNew->Divide(efficiencyTrackingPosOld);
 685   efficiencyTrackingPosNew->DrawCopy();
 686   efficiencyTrackingNegNew->Divide(efficiencyTrackingNegOld);
 687   efficiencyTrackingNegNew->DrawCopy("SAME");
 688
 689
 690
 691 }
 692
 693
 694
 695 //_______________________________________________________________________
 696 TH1D * ExtractEfficiency(TList * listMC, Int_t sign, Int_t centrality, Bool_t isTOF) {
 697   //
 698   // Extract the efficiency
 699   //
 700   //
 701   // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
 702   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 703   // (2.) pT
 704   // (3.) sign
 705   // (4.) rapidity --> filled 4x
 706   // (5.)  pull TPC dEx --> filled 4x
 707   // (6.) has valid TOF pid signal
 708   // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
 709   // (8..) dca_xy
 710   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
 711   //
 712   THnSparseF * hist = (THnSparseF *) listMC->FindObject("fHistMCparticles");
 713   //
 714   // common selections for generated and reconstructed particles
 715   //
 716   hist->GetAxis(0)->SetRangeUser(0,0);                     // select deuterons
 717   hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
 718   hist->GetAxis(3)->SetRangeUser(sign,sign);               // select anti-deuterons
 719   //assume eff. indepedent of centrality  ---> hist->GetAxis(1)->SetRangeUser(centrality,centrality);   // select centrality bin
 720   //
 721   hist->GetAxis(9)->SetRangeUser(0,0);
 722   TH1D * generated = hist->Projection(2);
 723   generated->SetNameTitle(Form("generated_cent%i_sign%i",centrality,sign),
 724                           Form("generated_cent%i_sign%i",centrality,sign));
 725   //
 726   // further cuts on reconstructed particles
 727   //
 728   hist->GetAxis(8)->SetRangeUser(-0.5,0.5);  // DCA-cut / TODO: do proper unfolding
 729   //
 730   hist->GetAxis(5)->SetRangeUser(-4.5,4.5);  // TPC-PID cut
 731   //
 732   if (isTOF) {
 733     //    hist->GetAxis(5)->SetRangeUser(-2.,2.);  // TPC-PID cut
 734     hist->GetAxis(6)->SetRangeUser(1.,1.);  // TOF-PID cut
 735     hist->GetAxis(7)->SetRangeUser(-1.2,1.2);  // TOF-PID cut
 736   }
 737   //
 738   hist->GetAxis(9)->SetRangeUser(1,1);
 739   TH1D * rawSpec = hist->Projection(2);
 740   rawSpec->Sumw2();
 741   rawSpec->SetNameTitle(Form("eff_cent%i_sign%i",centrality,sign),
 742                         Form("eff_cent%i_sign%i",centrality,sign));
 743   //
 744   rawSpec->Divide(generated);
 745   //
 746   if (sign > 0) {
 747     rawSpec->SetMarkerColor(kRed);
 748     rawSpec->SetLineColor(kRed);
 749   } else{
 750     rawSpec->SetMarkerColor(kBlue);
 751     rawSpec->SetLineColor(kBlue);
 752   }
 753   //
 754   // reset ranges
 755   //
 756   hist->GetAxis(5)->SetRangeUser(-1000.,1000);  // TPC-PID cut
 757   hist->GetAxis(6)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 758   hist->GetAxis(7)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 759   hist->GetAxis(8)->SetRangeUser(-1000.,1000);  // TOF-PID cut
 760   //
 761   //
 762   //
 763   return rawSpec;
 764
 765 }
 766
 767 //__________________________________________________________________
 768 TF1 * MakeEtaRapidityGenCorrection(Bool_t drawQA) {
 769   //
 770   // the enhanced sample LHC11b9_1 has a stupid cut on the generated level
 771   //
 772   TF1 * funcDeut = new TF1("funcDeut","TMath::ASinH(TMath::Sqrt(1 + (1.876*1.876)/(x*x))*TMath::SinH(0.5))",0.5,5);
 773   TF1 * etaCut = new TF1("etaCut","0.9",0.5,5);
 774   //
 775   funcDeut->SetLineColor(kBlue);
 776   //
 777   if (drawQA) {
 778     TCanvas * canv1 = new TCanvas("canv1","original functions");
 779     funcDeut->DrawCopy();
 780     funcDeut->GetHistogram()->GetYaxis()->SetRange(0,1.8);
 781     etaCut->DrawCopy("same");
 782   }
 783   //
 784   TF1 * funcCorr = new TF1("funcCorr","funcDeut/etaCut*(x - 1.1 < 0) + 1*(x - 1.1 >= 0)",0.5,2);
 785   if (drawQA) {
 786     TCanvas * canv2 = new TCanvas("canv2","correction function");
 787     funcCorr->Draw();
 788   }
 789   //
 790   return funcCorr;
 791
 792 }
 793
 794
 795 //__________________________________________________________________
 796 void MakeQAplots() {
 797   //
 798   // make the QA plots
 799   //
 800   TCanvas * canvData          =  PlotTpcQA(fileNameData,          kFALSE);
 801   TCanvas * canvMcEnhancedNew =  PlotTpcQA(fileNameMcEnhancedNew, kTRUE);
 802   TCanvas * canvMcUnEnhanced  =  PlotTpcQA(fileNameMcUnEnhanced,  kTRUE, kTRUE);
 803   TCanvas * canvMcEnhancedOld =  PlotTpcQA(fileNameMcEnhancedOld, kTRUE, kTRUE);
 804   //
 805   TCanvas * canvQA = new TCanvas("canvQA","canvQA");
 806   canvQA->Print("QAplots/QAplots.pdf(");
 807   //
 808   // add all QA plots to file
 809   //
 810   canvData->Print("QAplots/QAplots.pdf");
 811   canvMcEnhancedNew->Print("QAplots/QAplots.pdf");
 812   canvMcUnEnhanced->Print("QAplots/QAplots.pdf");
 813   canvMcEnhancedOld->Print("QAplots/QAplots.pdf");
 814   //
 815   canvQA->Print("QAplots/QAplots.pdf)");
 816   //
 817   // spectra QA
 818   //
 819   TFile inFile(fileNameData);
 820   TList * list = (TList * ) inFile.Get("akalweit_Nuclei");
 821   PlotSpectraQA(list);
 822
 823
 824 }
 825
 826
 827 //__________________________________________________________________
 828 TCanvas * PlotTpcQA(Char_t * fileName, Bool_t isMC, Bool_t isMCold) {
 829   //
 830   // plot QA of TPC-PID
 831   //
 832   TFile inFileMC(fileName);
 833   TList * list = (TList * ) inFileMC.Get("akalweit_Nuclei");
 834   //
 835   TF1 * fBBlines[6];
 836   const Double_t masses[6] = {0.000511, 0.138, 0.4936, 0.938, 2*0.938, 3*0.938};
 837   //
 838   // set BB parameters
 839   //
 840   Double_t parData[5] = {1.45802, 27.4992, 4.00313e-15, 2.48485, 8.31768};
 841   Double_t parMC[5] = {1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729};
 842   Double_t parMCold[5] = {1.17329, 27.4992, 4.00313e-15, 2.35563, 9.47569}; // OLD FOR LHC11b9_1 !!
 843   //
 844   for(Int_t ii = 0; ii < 6; ii++) {
 845     fBBlines[ii] = new TF1(Form("fBB_%i",ii),"AliExternalTrackParam::BetheBlochAleph(x/([0]),[1],[2],[3], [4], [5])",0.3,10);
 846
 847     fBBlines[ii]->SetParameters(masses[ii], parData[0], parData[1], parData[2], parData[3], parData[4]);
 848     if (isMC) {
 849       fBBlines[ii]->SetParameters(masses[ii], parMC[0], parMC[1], parMC[2], parMC[3], parMC[4]);
 850       if (isMCold) fBBlines[ii]->SetParameters(masses[ii], parMCold[0], parMCold[1], parMCold[2], parMCold[3], parMCold[4]);
 851     }
 852     fBBlines[ii]->SetLineColor(kRed);
 853     fBBlines[ii]->SetLineWidth(2);
 854   }
 855   //
 856   TCanvas * canvDeDx = new TCanvas(Form("canvDeDxTpcQA_%s",fileName),"control histogram for dE/dx");
 857   canvDeDx->Divide(1,2);
 858   canvDeDx->cd(1);
 859   TH3D * histTPC = (TH3D *) list->FindObject("fHistPidQA");
 860   histTPC->GetZaxis()->SetRangeUser(-1,-1);
 861   TH2D * histNegTPC = (TH2D *) histTPC->Project3D("yx");
 862   histNegTPC->SetNameTitle("histNegTPC", Form("histNegTPC_%s", fileName));
 863   gPad->SetLogx();
 864   histNegTPC->DrawCopy("colZ");
 865   for(Int_t ii = 0; ii < 6; ii++) fBBlines[ii]->DrawCopy("same");
 866   //
 867   canvDeDx->cd(2);
 868   histTPC->GetZaxis()->SetRangeUser(+1,+1);
 869   TH2D * histPosTPC = (TH2D*) histTPC->Project3D("yx");
 870   histPosTPC->SetNameTitle("histPosTPC", Form("histPosTPC_%s", fileName));
 871   gPad->SetLogx();
 872   histPosTPC->DrawCopy("colZ");
 873   //  for(Int_t ii = 0; ii < 6; ii++) fBBlines[ii]->DrawCopy("same");
 874   fBBlines[4]->DrawCopy("same"); // Draw deuteron line
 875   //
 876   // He-3 goes extra
 877   //
 878   TF1 * funcHel3 = new TF1("funcHel3","4*AliExternalTrackParam::BetheBlochAleph(2*x/(0.938*3),1.74962,27.4992,4.00313e-15,2.42485,8.31768)",0.2,6.);
 879   if (isMC) funcHel3 = new TF1("funcHel3MC","4*AliExternalTrackParam::BetheBlochAleph(2*x/(0.938*3), 1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729)",0.2,6.);
 880   //
 881   fBBlines[5]->DrawCopy("same"); // Draw deuteron line
 882   funcHel3->DrawCopy("same");
 883
 884   return canvDeDx;
 885
 886 }
 887
 888
 889
 890 //_______________________________________________________________________
 891 void PlotSpectraQA(TList * list, Int_t particleType, Int_t sign) {
 892   //
 893   // Make some basic QA plots for TOF related PID
 894   //
 895   THnSparseF * hist = (THnSparseF *) list->FindObject("fHistRealTracks");
 896   hist->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select Rapidity range
 897   hist->GetAxis(0)->SetRangeUser(particleType,particleType); // select deuterons
 898   hist->GetAxis(3)->SetRangeUser(sign,sign); // select deuterons
 899   //
 900   TCanvas * canvSpectraQA1 = new TCanvas("canvSpectraQA1","QA for nsigma TPC");
 901   TH2F * histTpcNsigma = (TH2F *) hist->Projection(5,2);
 902   histTpcNsigma->SetNameTitle("histTpcNsigma","histTpcNsigma");
 903   histTpcNsigma->DrawCopy("colZ");
 904   //
 905   TCanvas * canvSpectraQA2 = new TCanvas("canvSpectraQA2","QA for nsigma TOF");
 906   hist->GetAxis(6)->SetRangeUser(1,1);   // require proper match to TOF
 907   hist->GetAxis(5)->SetRangeUser(-4.,4); // pre-select tracks with TPC-PID
 908   TH2F * histTofNsigma = (TH2F *) hist->Projection(7,2);
 909   histTofNsigma->RebinY(4);
 910   histTofNsigma->SetNameTitle("histTofNsigma","histTofNsigma");
 911   histTofNsigma->DrawCopy("colZ");
 912
 913 }
 914
 915
 916
 917 //_______________________________________________________________________
 918 void MakeMaterialCorrection(){
 919   //
 920   // compare fitted values with MC pt-dependence
 921   //
 922   TFile outFile("output/materialCorrection.root","RECREATE");
 923   outFile.Close();
 924   //
 925   Int_t kMaxCentrality = 5;
 926   TFile * fileIn = TFile::Open("output/efficiencies.root");
 927   //
 928   //
 929   for(Int_t iCentr = 0; iCentr < kMaxCentrality; iCentr++) {
 930     //
 931     TH1D * matCorrMC = (TH1D *) fileIn->Get("efficiencyTrackingNegNew");
 932     matCorrMC->Reset();
 933     matCorrMC->SetNameTitle(Form("matCorrMC_%i",iCentr), Form("matCorrMC_i",iCentr));
 934     //
 935     TH1D * matCorr = (TH1D *) fileIn->Get("efficiencyTrackingNegNew");
 936     matCorr->SetLineColor(kRed);
 937     matCorr->Reset();
 938     matCorr->SetNameTitle(Form("matCorr_%i",iCentr), Form("matCorr_i",iCentr));
 939     //
 940     for(Int_t iBin = 0; iBin < matCorrMC->GetXaxis()->GetNbins(); iBin++) {
 941       Double_t pT = matCorrMC->GetXaxis()->GetBinCenter(iBin);
 942       if (pT > 0.6 && pT < 2.8) {
 943         Float_t corr = GetMaterialCorrection(pT,iCentr,kTRUE);
 944         matCorrMC->SetBinContent(iBin, corr);
 945         if (pT < 2.) {
 946           Float_t corrReal = GetMaterialCorrection(pT,iCentr,kFALSE);
 947           matCorr->SetBinContent(iBin, corrReal);
 948         }
 949       }
 950     }
 951     //
 952     //
 953     TCanvas * canvMaterial = new TCanvas("canvMaterial","canvMaterial");
 954     canvMaterial->cd();
 955     if (iCentr == 0) {
 956       matCorr->DrawCopy();
 957       matCorrMC->DrawCopy("SAME");
 958     } else {
 959       matCorrMC->DrawCopy("SAME");
 960     }
 961     TFile * outputFile = new TFile("output/materialCorrection.root","UPDATE");
 962     matCorr->Write();
 963     outputFile->Close();
 964     delete outputFile;
 965   }
 966
 967
 968 }
 969
 970
 971
 972 //_______________________________________________________________________
 973 Float_t GetMaterialCorrection(Float_t ptBin, Int_t centralityBin, Bool_t pureMC) {
 974   //
 975   // subtract the material contamination from positive particles
 976   //
 977   //
 978   Float_t dcaCutInAnalysis = 0.5;
 979   //
 980   TFile * inFileData = TFile::Open("data/dataFinal.root");
 981   TList * listData = (TList *) inFileData->Get("akalweit_Nuclei");
 982   THnSparse * histData = (THnSparse *) listData->FindObject("fHistRealTracks");
 983   //
 984   TFile * inFileMC = TFile::Open("MC/McCombined.root");
 985   TList * listMC = (TList *) inFileMC->Get("akalweit_Nuclei");
 986   THnSparse * histMC = (THnSparse *) listMC->FindObject("fHistMCparticles");
 987   //
 988   //
 989   // (1.) select deuterons
 990   //
 991   histData->GetAxis(0)->SetRangeUser(0,0);
 992   histMC->GetAxis(0)->SetRangeUser(0,0);
 993   histData->GetAxis(3)->SetRangeUser(+1,+1); // sign
 994   histMC->GetAxis(3)->SetRangeUser(+1,+1);   // sign
 995   //
 996   histData->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
 997   histMC->GetAxis(4)->SetRangeUser(-0.49, 0.49);             // select rapidity range
 998   //
 999   //  hist->GetAxis(8)->SetRangeUser(-1.,1.);  // DCA-cut / TODO: do proper unfolding
1000   histData->GetAxis(5)->SetRangeUser(-2.5,2.5);  // TPC-PID cut
1001   histData->GetAxis(6)->SetRangeUser(1.,1.);  // TOF-PID cut
1002   histData->GetAxis(7)->SetRangeUser(-0.5,0.5);  // TOF-PID cut
1003   //
1004   // (2.) select pt-range
1005   //
1006   histData->GetAxis(2)->SetRangeUser(ptBin,ptBin);
1007   histMC->GetAxis(2)->SetRangeUser(ptBin,ptBin);
1008   //
1009   // (4.a) get different MC templates
1010   //
1011   histMC->GetAxis(9)->SetRangeUser(1,1);
1012   TH1D * prim = (TH1D*) histMC->Projection(8);
1013   prim->SetNameTitle("prim","prim");
1014   prim->SetLineColor(kRed);
1015   //
1016   histMC->GetAxis(9)->SetRangeUser(4,4);
1017   //
1018   // centrality can oonly be selected after the primary yield
1019   //
1020   Int_t centrality = 0;
1021   if (centralityBin == 0) centrality = 0; //  0-10%
1022   if (centralityBin == 1) centrality = 2; // 10-20%
1023   if (centralityBin == 2) centrality = 3; // 20-40%
1024   if (centralityBin == 3) centrality = 5; // 40-60%
1025   if (centralityBin == 4) centrality = 7; // 60-80%
1026   //
1027   histData->GetAxis(1)->SetRangeUser(centrality,centrality);       // select centrality bin
1028   histMC->GetAxis(1)->SetRangeUser(centrality,centrality);       // select centrality bin
1029   if (centralityBin != 1) {
1030     histData->GetAxis(1)->SetRangeUser(centrality,centrality+1);   // select centrality bin
1031     histMC->GetAxis(1)->SetRangeUser(centrality,centrality+1);     // select centrality bin
1032   }
1033   //
1034   TH1D * material = (TH1D*) histMC->Projection(8);
1035   histMC->GetAxis(2)->SetRangeUser(ptBin,ptBin);
1036   material->SetNameTitle("material","material");
1037   material->SetLineColor(kGreen);
1038   //
1039   TH1D * dcaDistr = (TH1D*) histData->Projection(8);
1040   //
1041   // if only MC, we just return the values
1042   //
1043   if (pureMC) {
1044     Int_t fractionBinLow = dcaDistr->FindBin(-dcaCutInAnalysis);
1045     Int_t fractionBinUp  = dcaDistr->FindBin(+dcaCutInAnalysis);
1046     Float_t materialFraction = 0.;
1047     prim->Add(material); // has to be scaled with number of events somehow
1048     materialFraction = material->Integral(fractionBinLow,fractionBinUp)/prim->Integral(fractionBinLow,fractionBinUp);
1049     //
1050     delete listData;
1051     delete listMC;
1052     inFileData->Close();
1053     inFileMC->Close();
1054     delete inFileData;
1055     delete inFileMC;
1056     //
1057     return materialFraction;
1058   }
1059   //
1060   // (6.) do the fit
1061   //
1062   Float_t dcaUp = 0.5;//0.0182 + 0.035/TMath::Power(ptBin,1.01);
1063   Float_t dcaLow  = -1.*dcaUp;
1064   Int_t binLow = dcaDistr->GetXaxis()->FindBin(dcaLow);
1065   Int_t binUp  = dcaDistr->GetXaxis()->FindBin(dcaUp);
1066   //
1067   TObjArray *mcShapes = new TObjArray(3);        // MC histograms are put in this array
1068   mcShapes->Add(prim);
1069   mcShapes->Add(material);
1070   //
1071   TFractionFitter* fit = new TFractionFitter(dcaDistr, mcShapes); // initialise
1072   fit->SetRangeX(binLow,binUp);
1073   fit->Constrain(1,0.,1.);
1074   fit->Constrain(2,0.,1.);
1075   //
1076   //
1077   //
1078   Int_t status = fit->Fit();
1079   if (status==0) fit->GetPlot()->Draw("same");
1080   //
1081   // (7.) Rescaling of histograms and final plots
1082   //
1083   TH1D * result = 0x0;
1084   Double_t yieldPrim, yieldSec, yieldMat, error;
1085   if (status == 0) {
1086     fit->GetResult(0,yieldPrim,error);
1087     fit->GetResult(1,yieldMat,error);
1088     result = (TH1D*) fit->GetPlot();
1089     result->SetLineWidth(3);
1090   }
1091   //
1092   prim->Scale(1./prim->Integral(binLow,binUp));
1093   if (status == 0) prim->Scale(yieldPrim*result->Integral(binLow,binUp));
1094   //
1095   material->Scale(1./material->Integral(binLow,binUp));
1096   if (status == 0) material->Scale(yieldMat*result->Integral(binLow,binUp));
1097   //
1098   NormalizeSpectrum(dcaDistr,1,1);
1099   NormalizeSpectrum(prim,1,1);
1100   NormalizeSpectrum(material,1,1);
1101   if (status == 0) NormalizeSpectrum(result,1,1);
1102   //
1103   TH1D * sumHist = (TH1D*) prim->Clone(); // result is not equal to weighted sum ---> see FractionFitter doku and paper (?)
1104   sumHist->SetNameTitle("sumHist","sumHist");
1105   sumHist->Add(material);
1106   sumHist->SetLineColor(kOrange);
1107   //
1108   TCanvas * canvFeed = new TCanvas("canvFeed", "feed down and material correction");
1109   dcaDistr->SetMarkerStyle(24);  dcaDistr->SetMarkerSize(1.4);
1110   dcaDistr->DrawCopy("Ep");
1111   prim->DrawCopy("same");
1112   material->DrawCopy("same");
1113   sumHist->DrawCopy("same");
1114   //
1115   // Final-Results
1116   //
1117   Int_t fractionBinLow = dcaDistr->FindBin(-dcaCutInAnalysis);
1118   Int_t fractionBinUp  = dcaDistr->FindBin(+dcaCutInAnalysis);
1119   //
1120   Float_t materialFraction = 0.;
1121   if (status == 0) {
1122     result->DrawCopy("same");
1123     materialFraction = material->Integral(fractionBinLow,fractionBinUp)/dcaDistr->Integral(fractionBinLow,fractionBinUp);
1124   }
1125   Printf("RESULTING MATERIAL FRACTION: %f", materialFraction);
1126   Printf("TEMPLATE  MATERIAL FRACTION: %f", yieldMat);
1127   //
1128   delete sumHist;
1129   delete result;
1130   delete mcShapes;
1131   delete fit;
1132   //  delete histData;
1133   //delete histMC;
1134   delete listData;
1135   delete listMC;
1136   inFileData->Close();
1137   inFileMC->Close();
1138   delete inFileData;
1139   delete inFileMC;
1140   //
1141   return materialFraction;
1142
1143 }