PWGPP/pid/computePriors.C

   1 computePriors(){
   2   char filename[100];
   3   sprintf(filename,"yourfile.root");
   4   TFile *f = new TFile(filename);
   5   TTree *t = f->Get("yourtree");
   6
   7   Float_t mass[7] = {0.000511,0.139570,0.493677,0.938272,1.877837,2.817402,1.408054};
   8
   9   Int_t nev = t->GetEntries();
  10
  11   printf("ntracks = %i\n",nev);
  12
  13   const Int_t maxstep = 10;
  14
  15   TH1D *hpriors[7][maxstep];
  16   TH1D *hratio[7];
  17
  18   for(Int_t i=0;i<7;i++){
  19     for(Int_t j=0;j<maxstep;j++){
  20       char histoname[100];
  21       sprintf(histoname,"priors%istep%i",i,j);
  22       hpriors[i][j] = new TH1D(histoname,histoname,10,0,10);
  23       if(j==0){
  24         sprintf(histoname,"ratio%i",i);
  25         hratio[i] = new TH1D(histoname,histoname,10,0,10);
  26         for(Int_t k=1;k <= 100;k++){
  27           hpriors[i][j]->SetBinContent(k,1);
  28         }
  29       }
  30     }
  31   }
  32
  33   for(Int_t isp=0;isp < 7;isp++) hratio[isp]->Sumw2();
  34
  35   // loop variable definition
  36   Float_t pt, eta, weight=1.0;
  37   Float_t ppi, pka, ppr, pel, pmu, pde, ptr, phe, ptot;
  38
  39   for(Int_t j=1; j <maxstep;j++){
  40     for(Int_t isp=0;isp<7;isp++) hratio[isp]->Reset();
  41     printf("step %i\n",j);
  42     for(Int_t i=0;i<nev;i++){
  43       t->GetEvent(i);
  44       Int_t sw = t->GetLeaf("kTOF")->GetValue() * t->GetLeaf("kTPC")->GetValue();
  45       if(sw){
  46         pt = t->GetLeaf("pt")->GetValue();
  47         eta = t->GetLeaf("eta")->GetValue();
  48         //      weight = 1.0; // if you have fill the tree with weights
  49         if(pt > 0 && pt < 10){
  50           ppi = t->GetLeaf("tofW")->GetValue(2)*t->GetLeaf("tpcW")->GetValue(2)*hpriors[1][j-1]->Interpolate(pt);
  51           pka = t->GetLeaf("tofW")->GetValue(3)*t->GetLeaf("tpcW")->GetValue(3)*hpriors[2][j-1]->Interpolate(pt);
  52           ppr = t->GetLeaf("tofW")->GetValue(4)*t->GetLeaf("tpcW")->GetValue(4)*hpriors[3][j-1]->Interpolate(pt);
  53           pel = t->GetLeaf("tofW")->GetValue(0)*t->GetLeaf("tpcW")->GetValue(0)*hpriors[0][j-1]->Interpolate(pt);
  54           pmu = t->GetLeaf("tofW")->GetValue(1)*t->GetLeaf("tpcW")->GetValue(1)*hpriors[0][j-1]->Interpolate(pt);
  55           pde = t->GetLeaf("tofW")->GetValue(5)*t->GetLeaf("tpcW")->GetValue(5)*hpriors[4][j-1]->Interpolate(pt);
  56           ptr = t->GetLeaf("tofW")->GetValue(6)*t->GetLeaf("tpcW")->GetValue(6)*hpriors[5][j-1]->Interpolate(pt);
  57           phe = t->GetLeaf("tofW")->GetValue(7)*t->GetLeaf("tpcW")->GetValue(7)*hpriors[6][j-1]->Interpolate(pt);
  58           ptot = ppi+pka+ppr+pel+pde+ptr+phe+pmu;
  59
  60           if(ptot > 0){
  61             // fill ratio for |y|<0.5
  62             if(TMath::Abs(eta2y(pt,mass[1],eta)) < 0.5) hratio[1]->Fill(pt,ppi/ptot*weight);
  63             if(TMath::Abs(eta2y(pt,mass[2],eta)) < 0.5) hratio[2]->Fill(pt,pka/ptot*weight);
  64             if(TMath::Abs(eta2y(pt,mass[3],eta)) < 0.5) hratio[3]->Fill(pt,ppr/ptot*weight);
  65             if(TMath::Abs(eta2y(pt,mass[0],eta)) < 0.5) hratio[0]->Fill(pt,pel/ptot*weight);
  66             if(TMath::Abs(eta2y(pt,mass[4],eta)) < 0.5) hratio[4]->Fill(pt,pde/ptot*weight);
  67             if(TMath::Abs(eta2y(pt,mass[5],eta)) < 0.5) hratio[5]->Fill(pt,ptr/ptot*weight);
  68             if(TMath::Abs(eta2y(pt,mass[6],eta)) < 0.5) hratio[6]->Fill(pt,phe/ptot*weight);
  69
  70             // fill priors assuming #eta cut during the fill of the tree
  71             hpriors[1][j]->Fill(pt,ppi/ptot*weight);
  72             hpriors[2][j]->Fill(pt,pka/ptot*weight);
  73             hpriors[3][j]->Fill(pt,ppr/ptot*weight);
  74             hpriors[0][j]->Fill(pt,pel/ptot*weight);
  75             hpriors[4][j]->Fill(pt,pde/ptot*weight);
  76             hpriors[5][j]->Fill(pt,ptr/ptot*weight);
  77             hpriors[6][j]->Fill(pt,phe/ptot*weight);
  78           }
  79         }
  80       }
  81     }
  82
  83     for(Int_t isp=0;isp < 7;isp++) hpriors[isp][j]->Sumw2();
  84
  85     // Normalize abundances to the pion one
  86     hpriors[0][j]->Divide(hpriors[1][j]);
  87     hpriors[2][j]->Divide(hpriors[1][j]);
  88     hpriors[3][j]->Divide(hpriors[1][j]);
  89     hpriors[4][j]->Divide(hpriors[1][j]);
  90     hpriors[5][j]->Divide(hpriors[1][j]);
  91     hpriors[6][j]->Divide(hpriors[1][j]);
  92     hpriors[1][j]->Divide(hpriors[1][j]);
  93
  94     // make ratios w.r.t. pions (|y|<0.5)
  95     hratio[0]->Divide(hratio[1]);
  96     hratio[2]->Divide(hratio[1]);
  97     hratio[3]->Divide(hratio[1]);
  98     hratio[4]->Divide(hratio[1]);
  99     hratio[5]->Divide(hratio[1]);
 100     hratio[6]->Divide(hratio[1]);
 101     hratio[1]->Divide(hratio[1]);
 102   }
 103
 104   // Write output
 105   sprintf(filename,"priors.root");
 106   TFile *fout = new TFile(filename,"RECREATE");
 107   for(Int_t j=maxstep-2; j <maxstep;j++){ // last 2 steps
 108     hpriors[0][j]->Write();
 109     hpriors[1][j]->Write();
 110     hpriors[2][j]->Write();
 111     hpriors[3][j]->Write();
 112     hpriors[4][j]->Write();
 113     hpriors[5][j]->Write();
 114     hpriors[6][j]->Write();
 115   }
 116   hratio[0]->Write();
 117   hratio[1]->Write();
 118   hratio[2]->Write();
 119   hratio[3]->Write();
 120   hratio[4]->Write();
 121   hratio[5]->Write();
 122   hratio[6]->Write();
 123   fout->Close();
 124 }
 125 Double_t y2eta(Double_t pt, Double_t m, Double_t y){
 126   Double_t mt=TMath::Sqrt(m*m+pt*pt);
 127   return TMath::ASinH(mt/pt*TMath::SinH(y));
 128 }
 129 Double_t eta2y(Double_t pt, Double_t m, Double_t eta){
 130   Double_t mt=TMath::Sqrt(m*m+pt*pt);
 131   return TMath::ASinH(pt/mt*TMath::SinH(eta));
 132 }