PWGCF/FEMTOSCOPY/macros/MakeWeightFile.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6
   7 #include "TVector2.h"
   8 #include "TFile.h"
   9 #include "TString.h"
  10 #include "TF1.h"
  11 #include "TH1.h"
  12 #include "TH2.h"
  13 #include "TH3.h"
  14 #include "TProfile.h"
  15 #include "TProfile2D.h"
  16 #include "TMath.h"
  17 #include "TText.h"
  18 #include "TRandom3.h"
  19 #include "TArray.h"
  20 #include "TLegend.h"
  21 #include "TStyle.h"
  22 #include "TMinuit.h"
  23
  24 using namespace std;
  25
  26
  27 void MakeWeightFile()
  28 {
  29
  30   gStyle->SetOptStat(0);
  31   gStyle->SetOptDate(0);
  32
  33
  34   //TFile *InputFile = new TFile("Results/RawWeightFile_11h.root","READ");
  35   //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11.root","READ");
  36   //TFile *InputFile = new TFile("Results/RawWeightFile_genSignal_Rinv11_Smeared.root","READ");
  37   //TFile *InputFile = new TFile("Results/RawWeightFile_11h_PID1p5.root","READ");
  38   //TFile *InputFile = new TFile("Results/RawWeightFile_FB5.root","READ");
  39   //TFile *InputFile = new TFile("Results/RawWeightFile_11h_Chi3p1.root","READ");
  40   //TFile *InputFile = new TFile("Results/PDC_11h_standard_and_Raw0p04TTC.root","READ");// standard weights
  41   //TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_new.root","READ");
  42   TFile *InputFile = new TFile("Results/RawWeightFile_11h_4ktbins_2sigma_3sigmaTTC.root","READ");// _1(2sigmaTTC), _2(3sigmaTTC)
  43
  44   //TFile *InputFile = new TFile("MyOutput.root","READ");
  45   TDirectoryFile *tdir = (TDirectoryFile*)InputFile->Get("PWGCF.outputChaoticityAnalysis.root");
  46   TList *MyList=(TList*)tdir->Get("ChaoticityOutput_2");
  47   //TList *MyList=(TList*)InputFile->Get("MyList");
  48   InputFile->Close();
  49
  50   const int KtBins=4;
  51   const int KyBins=1;
  52   const int EDBins=1;
  53   const int MBins=10;//10, make sure MB assignment below is correct!!!!!!!!!!!!!!!
  54
  55   TFile *OutFile = new TFile("WeightFile_temp.root","RECREATE");
  56   TH3F *WeightHistos[KtBins][MBins];
  57
  58   for(int ktB=1; ktB<=KtBins; ktB++){
  59     for(int MB=1; MB<=MBins; MB++){
  60       TString *InNameNum = new TString("TwoPart_num_Kt_");
  61       *InNameNum += ktB-1;
  62       InNameNum->Append("_Ky_0_M_");
  63       if(MB<=10) *InNameNum += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
  64       else *InNameNum += 1;
  65       InNameNum->Append("_ED_0");
  66       //
  67       TString *InNameDen = new TString("TwoPart_den_Kt_");
  68       *InNameDen += ktB-1;
  69       InNameDen->Append("_Ky_0_M_");
  70       if(MB<=10) *InNameDen += MB-1;// make sure MB assignment is correct!!!!!!!!!!!!!!!
  71       else *InNameDen += 1;
  72       InNameDen->Append("_ED_0");
  73       //
  74       TH3D *tempNum = (TH3D*)MyList->FindObject(InNameNum->Data());
  75       TH3D *tempDen = (TH3D*)MyList->FindObject(InNameDen->Data());
  76       //
  77       int NormBinStart = tempNum->GetXaxis()->FindBin(0.135);//0.135
  78       int NormBinEnd = tempNum->GetXaxis()->FindBin(0.2);//0.2
  79       double Norm = tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
  80       Norm /= tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd);
  81       cout<<"Normalization = "<<Norm<<endl;
  82       cout<<tempNum->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<"  "<<tempDen->Integral(NormBinStart,NormBinEnd, NormBinStart,NormBinEnd, NormBinStart,NormBinEnd)<<endl;
  83       //
  84       TString *OutNameWeight = new TString("Weight_Kt_");
  85       *OutNameWeight += ktB-1;
  86       OutNameWeight->Append("_Ky_0_M_");
  87       *OutNameWeight += MB-1;
  88       OutNameWeight->Append("_ED_0");
  89
  90       int Nbins=tempNum->GetNbinsX();
  91       double QLimit = tempNum->GetXaxis()->GetBinUpEdge(Nbins);
  92       WeightHistos[ktB-1][MB-1] = new TH3F(OutNameWeight->Data(),"r3 Weights", Nbins,0,QLimit, Nbins,0,QLimit, Nbins,0,QLimit);
  93       WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitle("out");
  94       WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitle("side");
  95       WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitle("long");
  96       WeightHistos[ktB-1][MB-1]->GetXaxis()->SetTitleOffset(1.8);
  97       WeightHistos[ktB-1][MB-1]->GetYaxis()->SetTitleOffset(1.8);
  98       WeightHistos[ktB-1][MB-1]->GetZaxis()->SetTitleOffset(1.8);
  99       for(int outB=1; outB<=Nbins; outB++){
 100         for(int sideB=1; sideB<=Nbins; sideB++){
 101           for(int longB=1; longB<=Nbins; longB++){
 102             double weight=1, weight_e=0;
 103             if(tempDen->GetBinContent(outB,sideB,longB) > 0 && tempNum->GetBinContent(outB,sideB,longB) > 0) {
 104               //if(outB==sideB && outB==longB) cout<<outB<<"  "<<tempNum->GetBinContent(outB,sideB,longB)<<"  "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
 105               weight = double(tempNum->GetBinContent(outB,sideB,longB))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm;
 106               weight_e = pow(sqrt(double(tempNum->GetBinContent(outB,sideB,longB)))/double(tempDen->GetBinContent(outB,sideB,longB)) / Norm,2);
 107               weight_e += pow(sqrt(double(tempDen->GetBinContent(outB,sideB,longB)))*double(tempNum->GetBinContent(outB,sideB,longB))/pow(double(tempDen->GetBinContent(outB,sideB,longB)),2) / Norm,2);
 108               weight_e = sqrt(weight_e);
 109               //if(weight < 0.8 && weight > 0.1) cout<<outB<<"  "<<sideB<<"  "<<longB<<"  "<<tempNum->GetBinContent(outB,sideB,longB)<<"  "<<tempDen->GetBinContent(outB,sideB,longB)<<endl;
 110             }
 111             if(weight==0){
 112               WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, 0);
 113               WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, 0);
 114             }else {
 115               WeightHistos[ktB-1][MB-1]->SetBinContent(outB,sideB,longB, weight-1.0);// difference from unity
 116               WeightHistos[ktB-1][MB-1]->SetBinError(outB,sideB,longB, weight_e);
 117             }
 118
 119           }
 120         }
 121       }
 122
 123       WeightHistos[ktB-1][MB-1]->Write();
 124
 125     }// MB
 126   }// ktB
 127
 128   int BOI1=2;
 129   int BOI2=2;
 130   TH3D *histoOI = WeightHistos[0][0]->Clone();
 131   histoOI->SetDirectory(0);
 132   //OutFile->Close();
 133
 134   TH1D *pro = WeightHistos[0][0]->ProjectionY("pro",BOI1,BOI1,BOI2,BOI2);
 135   pro->Draw();
 136
 137   /*TFile *projections=new TFile("projections.root","RECREATE");
 138   histoOI->GetXaxis()->SetRange(BOI1,BOI1);
 139   TH2D *proYZ = (TH2D*)histoOI->Project3D("yz");
 140   histoOI->GetXaxis()->SetRange(1,40);
 141   histoOI->GetYaxis()->SetRange(BOI1,BOI1);
 142   TH2D *proXZ = (TH2D*)histoOI->Project3D("xz");
 143   histoOI->GetYaxis()->SetRange(1,40);
 144   histoOI->GetZaxis()->SetRange(BOI1,BOI1);
 145   TH2D *proXY = (TH2D*)histoOI->Project3D("xy");
 146   proYZ->GetXaxis()->SetRangeUser(0,0.1);
 147   proYZ->GetYaxis()->SetRangeUser(0,0.1);
 148   proXZ->GetXaxis()->SetRangeUser(0,0.1);
 149   proXZ->GetYaxis()->SetRangeUser(0,0.1);
 150   proXY->GetXaxis()->SetRangeUser(0,0.1);
 151   proXY->GetYaxis()->SetRangeUser(0,0.1);
 152   proXY->Draw("lego");
 153
 154   proYZ->Write();
 155   proXZ->Write();
 156   proXY->Write();
 157   */
 158
 159 }