PWGLF/SPECTRA/PiKaPr/TestAOD/HighLevelQA/GenerateTPCTOFnsigma.C

   1 ////////////////////////////////////////////////////////////
   2 // GenerateTPCTOFnsigma.C                                 //
   3 // creates a histogram of the TPC+TOF nsigma distribution //
   4 // written by John Groh                                   //
   5 ////////////////////////////////////////////////////////////
   6
   7 void GenerateTPCTOFnsigma()
   8 {
   9   //-------------------------------------------
  10   // Get TPC+TOF nsigma projections from data -
  11   //-------------------------------------------
  12   const Int_t nPart = 3;
  13   TString Particle[nPart] = {"Pion","Kaon","Proton"};
  14   const Int_t Color[nPart] = {(Int_t)kBlack, (Int_t)kRed, (Int_t)kBlue};
  15   const Int_t nPt = 6;
  16   const Float_t Pt[nPt] = {0.6, 1.0, 1.4, 1.8, 2.2, 2.6};
  17
  18   // open file and get histo manager
  19   TFile * fin = TFile::Open("output/No_Outliers/AnResDATATrain12_NoOutliers.root");
  20   TDirectoryFile * dir = (TDirectoryFile*)fin->Get("OutputAODSpectraTask_Data_Cent30to40_QVec0.0to100.0_Eta-0.8to0.8_3.0SigmaPID_TrBit1024");
  21   AliSpectraAODHistoManager * hman = (AliSpectraAODHistoManager*)dir->Get("SpectraHistos");
  22
  23   // get nsigma projections
  24   TH2F * hTPCTOFnsigmaDATA[nPart];
  25   TH1F * hTPCTOFnsigmaDATAProj[nPart][nPt];
  26   for (Int_t ipart=0; ipart<nPart; ipart++)
  27     {
  28       // get 2d histo
  29       hTPCTOFnsigmaDATA[ipart] = (TH2F*)((TH2F*)hman->GetNSigHistogram(Form("hHistNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone();
  30       for (Int_t ipt=0; ipt<nPt; ipt++)
  31         {
  32           // project it for the given pt bin
  33           hTPCTOFnsigmaDATAProj[ipart][ipt] = (TH1F*)hTPCTOFnsigmaDATA[ipart]->ProjectionY(Form("%ss - p_{T} = %.1f GeV/c",Particle[ipart].Data(),Pt[ipt]),
  34                                                                                            hTPCTOFnsigmaDATA[ipart]->GetXaxis()->FindBin(Pt[ipt]),
  35                                                                                            hTPCTOFnsigmaDATA[ipart]->GetXaxis()->FindBin(Pt[ipt]));
  36         }
  37     }
  38
  39   //-----------------------------------------------
  40   // Create TPC+TOF nsigma distribution histogram -
  41   //-----------------------------------------------
  42   // file for saving histogram
  43   TFile * fout = new TFile("TPCTOFnsigma.root","RECREATE");
  44
  45   // TPC nsigma distributions - just  gaussians
  46   TF1 * fTPCnsigma[nPart];
  47   for (Int_t ipart=0; ipart<nPart; ipart++)
  48     {
  49       fTPCnsigma[ipart] = new TF1("fTPCnsigma","gaus",-25,25);
  50       fTPCnsigma[ipart]->SetParameter(0,1); // normalization
  51     }
  52   fTPCnsigma[0]->SetParameter(1,0.0); // pion mean
  53   fTPCnsigma[0]->SetParameter(2,0.98); // pion sigma
  54   fTPCnsigma[1]->SetParameter(1,-0.1); // kaon mean
  55   fTPCnsigma[1]->SetParameter(2,0.98); // kaon sigma
  56   fTPCnsigma[2]->SetParameter(1,0.2); // proton mean
  57   fTPCnsigma[2]->SetParameter(2,0.72); // proton sigma
  58
  59   // TOF nsigma distributions
  60   gROOT->LoadMacro("/opt/alice/aliroot/trunk/src/PWGLF/SPECTRA/PiKaPr/TOF/PbPb276/macros/TOFsignal.C");
  61   TF1 *fTOFnsigma[nPart];
  62   for (Int_t ipart=0; ipart<nPart; ipart++)
  63     {
  64       fTOFnsigma[ipart] = new TF1(Form("fTOFnsigma_%s",Particle[ipart].Data()), TOFsignal, -25, 25, 4);
  65       fTOFnsigma[ipart]->SetParameter(0,1.0); // normalization
  66       fTOFnsigma[ipart]->SetParameter(3, 1.0); // tail
  67     }
  68   fTOFnsigma[0]->SetParameter(1, 0.0); // pion mean
  69   fTOFnsigma[0]->SetParameter(2, 1.02); // pion sigma
  70   fTOFnsigma[1]->SetParameter(1, -0.04); // kaon mean
  71   fTOFnsigma[1]->SetParameter(2, 1.05); // kaon sigma
  72   fTOFnsigma[2]->SetParameter(1, -0.18); // proton mean
  73   fTOFnsigma[2]->SetParameter(2, 0.82); // proton sigma
  74
  75
  76   // combined TPC+TOF nsigma distributions
  77   //  // to ensure that the binning is the same I'm copying an actual nsigma projection and emptying it before doing anything
  78   TH1F * hTPCTOFnsigma[nPart];
  79   for (Int_t ipart=0; ipart<nPart; ipart++)
  80     {
  81       //hTPCTOFnsigma[ipart] = (TH1F*)hTPCTOFnsigmaDATAProj[0][0]->Clone();
  82       //hTPCTOFnsigma[ipart]->Reset();
  83       hTPCTOFnsigma[ipart] = new TH1F("","",1000,0,25);
  84       hTPCTOFnsigma[ipart]->SetName(Form("hTPCTOFnsigma%s",Particle[ipart].Data()));
  85       hTPCTOFnsigma[ipart]->SetTitle("Predicted NSigma distributions (TPC+TOF);n#sigma;");
  86
  87       for (Int_t i=0; i<1000000; i++)
  88         {
  89           // get random numbers from the nsig distributions
  90           Float_t TPCpoint = fTPCnsigma[ipart]->GetRandom(-25,25);
  91           Float_t TOFpoint = fTOFnsigma[ipart]->GetRandom(-25,25);
  92
  93           // add them in quadrature and divide by root 2
  94           Float_t TPCTOFcomb = TMath::Sqrt( (TMath::Power(TPCpoint,2) + TMath::Power(TOFpoint,2))/2 );
  95
  96           // fill hTPCTOFnsigma w/ this
  97           hTPCTOFnsigma[ipart]->Fill(TPCTOFcomb);
  98         }
  99
 100       // normalize it
 101       Float_t integral = hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0),hTPCTOFnsigma[ipart]->FindBin(20));
 102       hTPCTOFnsigma[ipart]->Scale(1./integral);
 103     }
 104
 105   // draw  them
 106   TCanvas * cTPCTOFnsigmaPredicted = new TCanvas("cTPCTOFnsigmaPredicted","cTPCTOFnsigmaPredicted");
 107   gPad->SetLogy();
 108   for (Int_t ipart=0; ipart<nPart; ipart++)
 109     {
 110       hTPCTOFnsigma[ipart]->GetXaxis()->SetTitle("n#sigma");
 111       hTPCTOFnsigma[ipart]->GetXaxis()->SetRangeUser(-5,20);
 112       hTPCTOFnsigma[ipart]->SetLineColor(Color[ipart]);
 113       if (ipart == 0)
 114         {
 115           hTPCTOFnsigma[ipart]->DrawCopy("hist");
 116           TLegend * lTPCTOFnsigma = new TLegend(.69,.69,.99,.99);
 117           lTPCTOFnsigma->SetFillColor(0);
 118         }
 119       else hTPCTOFnsigma[ipart]->DrawCopy("histsame");
 120       lTPCTOFnsigma->AddEntry(hTPCTOFnsigma[ipart],Form("%ss",Particle[ipart].Data()),"l");
 121     }
 122   lTPCTOFnsigma->DrawClone();
 123
 124   // do the same for MC
 125   // TPC nsigma distributions - just  gaussians
 126   TF1 * fTPCnsigmaMC[nPart];
 127   for (Int_t ipart=0; ipart<nPart; ipart++)
 128     {
 129       fTPCnsigmaMC[ipart] = new TF1("fTPCnsigmaMC","gaus",-25,25);
 130       fTPCnsigmaMC[ipart]->SetParameter(0,1); // normalization
 131     }
 132   fTPCnsigmaMC[0]->SetParameter(1,-0.03); // pion mean
 133   fTPCnsigmaMC[0]->SetParameter(2,1.38); // pion sigma
 134   fTPCnsigmaMC[1]->SetParameter(1,0.14); // kaon mean
 135   fTPCnsigmaMC[1]->SetParameter(2,1.27); // kaon sigma
 136   fTPCnsigmaMC[2]->SetParameter(1,-0.17); // proton mean
 137   fTPCnsigmaMC[2]->SetParameter(2,1.03); // proton sigma
 138
 139   // TOF nsigma distributions
 140   //gROOT->LoadMacro("/opt/alice/aliroot/trunk/src/PWGLF/SPECTRA/PiKaPr/TOF/PbPb276/macros/TOFsignal.C");
 141   TF1 *fTOFnsigmaMC[nPart];
 142   for (Int_t ipart=0; ipart<nPart; ipart++)
 143     {
 144       fTOFnsigmaMC[ipart] = new TF1(Form("fTOFnsigmaMC_%s",Particle[ipart].Data()), TOFsignal, -25, 25, 4);
 145       fTOFnsigmaMC[ipart]->SetParameter(0,1.0); // normalization
 146       fTOFnsigmaMC[ipart]->SetParameter(3, 1.0); // tail
 147     }
 148   fTOFnsigmaMC[0]->SetParameter(1, -0.01); // pion mean
 149   fTOFnsigmaMC[0]->SetParameter(2, 0.94); // pion sigma
 150   fTOFnsigmaMC[1]->SetParameter(1, 0.02); // kaon mean
 151   fTOFnsigmaMC[1]->SetParameter(2, 0.84); // kaon sigma
 152   fTOFnsigmaMC[2]->SetParameter(1, 0.04); // proton mean
 153   fTOFnsigmaMC[2]->SetParameter(2, 0.66); // proton sigma
 154
 155
 156   // combined TPC+TOF nsigma distributions
 157   //  // to ensure that the binning is the same I'm copying an actual nsigma projection and emptying it before doing anything
 158   TH1F * hTPCTOFnsigmaMC[nPart];
 159   for (Int_t ipart=0; ipart<nPart; ipart++)
 160     {
 161       //hTPCTOFnsigmaMC[ipart] = (TH1F*)hTPCTOFnsigmaMCDATAProj[0][0]->Clone();
 162       //hTPCTOFnsigmaMC[ipart]->Reset();
 163       hTPCTOFnsigmaMC[ipart] = new TH1F("","",1000,0,25);
 164       hTPCTOFnsigmaMC[ipart]->SetName(Form("hTPCTOFnsigmaMC%s",Particle[ipart].Data()));
 165       hTPCTOFnsigmaMC[ipart]->SetTitle("Predicted NSigma distributions (TPC+TOF);n#sigma;");
 166
 167       for (Int_t i=0; i<1000000; i++)
 168         {
 169           // get random numbers from the nsig distributions
 170           Float_t TPCpoint = fTPCnsigmaMC[ipart]->GetRandom(-25,25);
 171           Float_t TOFpoint = fTOFnsigmaMC[ipart]->GetRandom(-25,25);
 172
 173           // add them in quadrature and divide by root 2
 174           Float_t TPCTOFcomb = TMath::Sqrt( (TMath::Power(TPCpoint,2) + TMath::Power(TOFpoint,2))/2 );
 175
 176           // fill hTPCTOFnsigmaMC w/ this
 177           hTPCTOFnsigmaMC[ipart]->Fill(TPCTOFcomb);
 178         }
 179
 180       // normalize it
 181       Float_t integral = hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0),hTPCTOFnsigmaMC[ipart]->FindBin(20));
 182       hTPCTOFnsigmaMC[ipart]->Scale(1./integral);
 183     }
 184
 185   // draw  them
 186   TCanvas * cTPCTOFnsigmaPredictedMC = new TCanvas("cTPCTOFnsigmaPredictedMC","cTPCTOFnsigmaPredictedMC");
 187   gPad->SetLogy();
 188   for (Int_t ipart=0; ipart<nPart; ipart++)
 189     {
 190       hTPCTOFnsigmaMC[ipart]->GetXaxis()->SetTitle("n#sigma");
 191       hTPCTOFnsigmaMC[ipart]->GetXaxis()->SetRangeUser(-5,20);
 192       hTPCTOFnsigmaMC[ipart]->SetLineColor(Color[ipart]);
 193       if (ipart == 0)
 194         {
 195           hTPCTOFnsigmaMC[ipart]->DrawCopy("hist");
 196           TLegend * lTPCTOFnsigmaMC = new TLegend(.69,.69,.99,.99);
 197           lTPCTOFnsigmaMC->SetFillColor(0);
 198         }
 199       else hTPCTOFnsigmaMC[ipart]->DrawCopy("histsame");
 200       lTPCTOFnsigmaMC->AddEntry(hTPCTOFnsigmaMC[ipart],Form("%ss",Particle[ipart].Data()),"l");
 201     }
 202   lTPCTOFnsigmaMC->DrawClone();
 203
 204
 205
 206   //---------------------------------------------------------------------------------------------------
 207   // Compare to the actual TPC+TOF nsigma distributions for different pt bins and different particles -
 208   //---------------------------------------------------------------------------------------------------
 209   TCanvas * cCompareToData[nPt];
 210   Int_t xPos = 50; // for canvas positioning
 211   Int_t yPos = 25; // for canvas positioning
 212   TH1F * hRatio[nPart][nPt];
 213   for (Int_t ipt=0; ipt<nPt; ipt++)
 214     {
 215       cCompareToData[ipt] = new TCanvas(Form("cCompareToData_%.1fPt",Pt[ipt]),Form("cCompareToData_%.1fPt",Pt[ipt]),xPos += 50, yPos += 25, 700, 500);
 216       cCompareToData[ipt]->Divide(3,2);
 217
 218       // direct comparison of the histograms
 219       for (Int_t ipart=0; ipart<nPart; ipart++)
 220         {
 221           cCompareToData[ipt]->cd(ipart+1);
 222           gPad->SetLogy();
 223
 224           // draw histo from data
 225           hTPCTOFnsigmaDATAProj[ipart][ipt]->SetLineColor(Color[ipart]);
 226           hTPCTOFnsigmaDATAProj[ipart][ipt]->SetTitle(Form("TPC+TOF n#sigma for %ss at p_{T} = %.1f GeV/c;n#sigma;",Particle[ipart].Data(),Pt[ipt]));
 227           hTPCTOFnsigmaDATAProj[ipart][ipt]->GetXaxis()->SetRangeUser(-5,10);
 228           hTPCTOFnsigmaDATAProj[ipart][ipt]->DrawCopy("hist");
 229
 230           // scale generated histo and draw it as well
 231           Float_t maxDATA = -1;
 232
 233           for (Int_t ibin=hTPCTOFnsigmaDATAProj[ipart][ipt]->FindBin(0); ibin<hTPCTOFnsigmaDATAProj[ipart][ipt]->FindBin(1.5); ibin++)
 234             {
 235               if (hTPCTOFnsigmaDATAProj[ipart][ipt]->GetBinContent(ibin) > maxDATA) maxDATA = hTPCTOFnsigmaDATAProj[ipart][ipt]->GetBinContent(ibin);
 236             }
 237           hTPCTOFnsigma[ipart]->Scale(maxDATA / hTPCTOFnsigma[ipart]->GetMaximum());
 238           hTPCTOFnsigma[ipart]->SetLineColor(kGreen+2);
 239           hTPCTOFnsigma[ipart]->DrawCopy("histsame");
 240
 241           // legend
 242           TLegend * lCompareToData = new TLegend(.6,.8,.95,.92);
 243           lCompareToData->SetFillColor(0);
 244           lCompareToData->AddEntry(hTPCTOFnsigmaDATAProj[ipart][ipt],"Data","l");
 245           lCompareToData->AddEntry(hTPCTOFnsigma[ipart],"Predicted","l");
 246           lCompareToData->DrawClone();
 247         }
 248
 249       // ratios of the histograms
 250       for (Int_t ipart=0; ipart<nPart; ipart++)
 251         {
 252           cCompareToData[ipt]->cd(ipart+1+nPart);
 253
 254           // compute the ratio
 255           hRatio[ipart][ipt] = (TH1F*)hTPCTOFnsigmaDATAProj[ipart][ipt]->Clone();
 256           hRatio[ipart][ipt]->Divide(hTPCTOFnsigma[ipart]);
 257           hRatio[ipart][ipt]->SetName("Data/Predicted");
 258           hRatio[ipart][ipt]->SetTitle("Data/Predicted;n#sigma;");
 259           hRatio[ipart][ipt]->SetLineColor(Color[ipart]);
 260           hRatio[ipart][ipt]->GetXaxis()->SetRangeUser(0,3);
 261           //      hRatio[ipart][ipt]->GetYaxis()->SetRangeUser(-2,10);
 262           hRatio[ipart][ipt]->DrawCopy("hist");
 263         }
 264     }
 265
 266
 267   /*
 268   // crazy experiment: see if I can get a TF1 with the full convolution
 269   TF1 *fTOFnsig=new TF1("fTOFnsig","(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]))",0,25);
 270   fTOFnsig->SetParameter(0,1);
 271   fTOFnsig->SetParameter(1,1);
 272   fTOFnsig->SetParameter(2,1);
 273   fTOFnsig->SetParameter(3,3);
 274
 275   TF1 * fTPCnsig = new TF1("fTPCnsig","gaus",0,25);
 276   fTPCnsig->SetParameter(0,1);
 277   fTPCnsig->SetParameter(1,0);
 278   fTPCnsig->SetParameter(2,1);
 279
 280   TF1 * fTPCTOFnsig = new TF1("fTPCTOFnsig","TMath::Sqrt( (TMath::Power(fTPCnsigma,2) + TMath::Power(fTOFnsig,2))/2 )",0,25);
 281
 282   //  TCanvas * c1 = new TCanvas();
 283   hTPCTOFnsigma->Fit(fTPCTOFnsig,"NM","",0,3);
 284   fTPCTOFnsig->DrawCopy("same");
 285
 286   Float_t integralOfFit = fTPCTOFnsig->Integral(0,3);
 287   cout << "The integral is " << integralOfFit << endl;
 288
 289   Float_t integralOfHist = hTPCTOFnsigma->Integral(hTPCTOFnsigma->FindBin(0), hTPCTOFnsigma->FindBin(3));
 290   cout << "The integral of the histo is " << integralOfHist << endl;
 291   */
 292
 293
 294   // print to a pdf for easy viewing
 295   cTPCTOFnsigmaPredicted->Print("TPCTOFnsigma.pdf(","pdf");
 296   for (Int_t ipt=0; ipt<nPt; ipt++)
 297     {
 298       if (ipt == nPt-1) cCompareToData[ipt]->Print("TPCTOFnsigma.pdf)","pdf");
 299       else cCompareToData[ipt]->Print("TPCTOFnsigma.pdf","pdf");
 300     }
 301
 302   // save it all to file
 303   fout->cd();
 304   cTPCTOFnsigmaPredicted->Write();
 305   cTPCTOFnsigmaPredictedMC->Write();
 306   for (Int_t ipt=0; ipt<nPt; ipt++) cCompareToData[ipt]->Write();
 307   fout->Close();
 308 }