PWGCF/FEMTOSCOPY/macros/MakeKFileTherm.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <complex>
   7
   8 #include "TVector2.h"
   9 #include "TFile.h"
  10 #include "TString.h"
  11 #include "TF1.h"
  12 #include "TH1.h"
  13 #include "TH2.h"
  14 #include "TH3.h"
  15 #include "TProfile.h"
  16 #include "TProfile2D.h"
  17 #include "TMath.h"
  18 #include "TText.h"
  19 #include "TRandom3.h"
  20 #include "TArray.h"
  21 #include "TLegend.h"
  22 #include "TStyle.h"
  23 #include "TMinuit.h"
  24
  25 using namespace std;
  26
  27 void MakeKFileTherm(){
  28   // 2-particles
  29
  30   const int kbVALUES=6;// 6 b values (2,3,5,7,8,9)
  31
  32   TFile *probeFile=new TFile("Therm_FSI_b2.root","READ");// all files have same binning.
  33   TH2D *probeHisto = (TH2D*)probeFile->Get("K2_ss");// kt x qinv
  34   int binsY= probeHisto->GetNbinsY();
  35   double lowY = probeHisto->GetYaxis()->GetBinLowEdge(1);
  36   double highY = probeHisto->GetYaxis()->GetBinUpEdge(binsY);
  37   probeFile->Close();
  38
  39   TFile *OutFile=new TFile("KFile_temp.root","RECREATE");
  40   TH2D *K2ssT = new TH2D("K2ssT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
  41   TH2D *K2osT = new TH2D("K2osT","",kbVALUES,0.5,kbVALUES+0.5, binsY,lowY,highY);// kt integrated
  42   TH3D *K2ssT_kt = new TH3D("K2ssT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
  43   TH3D *K2osT_kt = new TH3D("K2osT_kt","",kbVALUES,0.5,kbVALUES+0.5, 6,0.5,6+0.5, binsY,lowY,highY);// kt differential
  44   K2ssT_kt->GetXaxis()->SetTitle("b bin"); K2osT_kt->GetXaxis()->SetTitle("b bin");
  45   K2ssT_kt->GetYaxis()->SetTitle("kt bin"); K2osT_kt->GetYaxis()->SetTitle("kt bin");
  46   K2ssT_kt->GetZaxis()->SetTitle("qinv"); K2osT_kt->GetZaxis()->SetTitle("qinv");
  47
  48   for(int r=0; r<kbVALUES; r++){
  49
  50     TString *nameT=new TString("Therm_FSI_b");
  51
  52
  53     if(r==0) {*nameT += 2;}
  54     if(r==1) {*nameT += 3;}
  55     if(r==2) {*nameT += 5;}
  56     if(r==3) {*nameT += 7;}
  57     if(r==4) {*nameT += 8;}
  58     if(r==5) {*nameT += 9;}
  59
  60     nameT->Append(".root");
  61
  62
  63
  64     //
  65     TFile *file_T=new TFile(nameT->Data(),"READ");
  66     TH2D *Num2_ss = (TH2D*)file_T->Get("K2_ss");
  67     TH2D *Den2_ss = (TH2D*)file_T->Get("PlaneWF_ss");
  68     TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY();// kt integrated
  69     TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY();// kt integrated
  70     TH2D *Num2_os = (TH2D*)file_T->Get("K2_os");
  71     TH2D *Den2_os = (TH2D*)file_T->Get("PlaneWF_os");
  72     TH1D *Num_os = (TH1D*)Num2_os->ProjectionY();// kt integrated
  73     TH1D *Den_os = (TH1D*)Den2_os->ProjectionY();// kt integrated
  74     Num_ss->Divide(Den_ss);
  75     Num_os->Divide(Den_os);
  76     for(int i=1; i<=binsY; i++){
  77       K2ssT->SetBinContent(r+1, i, Num_ss->GetBinContent(i));
  78       K2osT->SetBinContent(r+1, i, Num_os->GetBinContent(i));
  79     }
  80
  81     // kt differential
  82     for(int ktbin=1; ktbin<=6; ktbin++){
  83       TString *name=new TString("pro_");
  84       *name += ktbin;
  85       *name += 1;
  86       TH1D *Num_ss = (TH1D*)Num2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
  87       *name += 2;
  88       TH1D *Den_ss = (TH1D*)Den2_ss->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
  89       *name += 3;
  90       TH1D *Num_os = (TH1D*)Num2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
  91       *name += 4;
  92       TH1D *Den_os = (TH1D*)Den2_os->ProjectionY(name->Data(), ktbin*2+3, ktbin*2+4);
  93       Num_ss->Divide(Den_ss);
  94       Num_os->Divide(Den_os);
  95        for(int i=1; i<=binsY; i++){
  96          K2ssT_kt->SetBinContent(r+1, ktbin, i, Num_ss->GetBinContent(i));
  97          K2osT_kt->SetBinContent(r+1, ktbin, i, Num_os->GetBinContent(i));
  98        }
  99     }
 100
 101     // 3-particles
 102
 103     TH3D *Num3_ss=(TH3D*)file_T->Get("K3ss_3D");
 104     TH3D *Den3_ss=(TH3D*)file_T->Get("PlaneWF3ss_3D");
 105     Num3_ss->Divide(Den3_ss);
 106     TString *OutNameSS = new TString("K3ss_");
 107     *OutNameSS += r;
 108     OutFile->cd();
 109     TH3D *K3ss = (TH3D*)Num3_ss->Clone();
 110     TH3D *temp3ss = (TH3D*)Num3_ss->Clone();
 111     for(int i=1; i<=K3ss->GetNbinsX(); i++){
 112       for(int j=1; j<=K3ss->GetNbinsY(); j++){
 113         for(int k=1; k<=K3ss->GetNbinsZ(); k++){
 114
 115           // GRS
 116           //double GRS = K2ssT->GetBinContent(r+1, i) * K2ssT->GetBinContent(r+1, j) * K2ssT->GetBinContent(r+1, k);
 117           //K3ss->SetBinContent(i,j,k, GRS);
 118
 119           // Omega0
 120           if(temp3ss->GetBinContent(i,j,k) > 1.0) K3ss->SetBinContent(i,j,k, 1.0);
 121           else{
 122             double mean=0;
 123             double terms=0;
 124             if(temp3ss->GetBinContent(i,j,k) < 1.0) {mean += temp3ss->GetBinContent(i,j,k); terms++;}
 125             if(temp3ss->GetBinContent(i,k,j) < 1.0) {mean += temp3ss->GetBinContent(i,k,j); terms++;}
 126             if(temp3ss->GetBinContent(j,i,k) < 1.0) {mean += temp3ss->GetBinContent(j,i,k); terms++;}
 127             if(temp3ss->GetBinContent(j,k,i) < 1.0) {mean += temp3ss->GetBinContent(j,k,i); terms++;}
 128             if(temp3ss->GetBinContent(k,i,j) < 1.0) {mean += temp3ss->GetBinContent(k,i,j); terms++;}
 129             if(temp3ss->GetBinContent(k,j,i) < 1.0) {mean += temp3ss->GetBinContent(k,j,i); terms++;}
 130
 131             if(terms > 0) {mean /= terms; K3ss->SetBinContent(i,j,k, mean);}
 132             else K3ss->SetBinContent(i,j,k, 0);
 133             }
 134
 135         }
 136       }
 137     }
 138     K3ss->Write(OutNameSS->Data());
 139     //
 140     TH3D *Num3_os=(TH3D*)file_T->Get("K3os_3D");
 141     TH3D *Den3_os=(TH3D*)file_T->Get("PlaneWF3os_3D");
 142     Num3_os->Divide(Den3_os);
 143     //
 144     TString *OutNameOS = new TString("K3os_");
 145     *OutNameOS += r;
 146     OutFile->cd();
 147     TH3D *K3os = (TH3D*)Num3_os->Clone();
 148     TH3D *temp3os = (TH3D*)Num3_os->Clone();
 149     for(int i=1; i<=K3os->GetNbinsX(); i++){
 150       for(int j=1; j<=K3os->GetNbinsY(); j++){
 151         for(int k=1; k<=K3os->GetNbinsZ(); k++){
 152
 153           // GRS
 154           //double GRS = K2ssT->GetBinContent(r+1, i) * K2osT->GetBinContent(r+1, j) * K2osT->GetBinContent(r+1, k);
 155           //K3os->SetBinContent(i,j,k, GRS);
 156
 157           // Omega0
 158           if(temp3os->GetBinContent(i,j,k) > 3.0) K3os->SetBinContent(i,j,k, 1.0);
 159           else {
 160             double mean=0;
 161             double terms=0;
 162             if(temp3os->GetBinContent(i,j,k) < 3.0) {mean += temp3os->GetBinContent(i,j,k); terms++;}
 163             if(temp3os->GetBinContent(i,k,j) < 3.0) {mean += temp3os->GetBinContent(i,k,j); terms++;}
 164
 165             if(terms > 0) {mean /= terms; K3os->SetBinContent(i,j,k, mean);}
 166             else K3os->SetBinContent(i,j,k, 0);
 167             }
 168
 169         }
 170       }
 171     }
 172
 173     K3os->Write(OutNameOS->Data());
 174
 175
 176     file_T->Close();
 177
 178   }// r loop
 179
 180   //
 181   OutFile->cd();
 182   //
 183   K2ssT->Write("K2ssT");
 184   K2osT->Write("K2osT");
 185   K2ssT_kt->Write("K2ssT_kt");
 186   K2osT_kt->Write("K2osT_kt");
 187   //
 188
 189   OutFile->Close();
 190
 191 }