PWGCF/FEMTOSCOPY/macros/Plot_muons.C

   1 #include <math.h>
   2 #include <time.h>
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5 #include <Riostream.h>
   6 #include <assert.h>
   7
   8 #include "TVector2.h"
   9 #include "TFile.h"
  10 #include "TString.h"
  11 #include "TF1.h"
  12 #include "TF3.h"
  13 #include "TH1.h"
  14 #include "TH2.h"
  15 #include "TH3.h"
  16 #include "TProfile.h"
  17 #include "TProfile2D.h"
  18 #include "TMath.h"
  19 #include "TText.h"
  20 #include "TRandom3.h"
  21 #include "TArray.h"
  22 #include "TLegend.h"
  23 #include "TStyle.h"
  24 #include "TMinuit.h"
  25 #include "TCanvas.h"
  26 #include "TLatex.h"
  27 #include "TPaveStats.h"
  28 #include "TASImage.h"
  29
  30 #define BohrR 1963.6885 // Bohr Radius for pions
  31 #define FmToGeV 0.19733 // conversion to fm
  32 #define PI 3.1415926
  33 #define masspiC 0.1395702 // pi+ mass (GeV/c^2)
  34 #define kappa3 0.16
  35 #define kappa4 0.40
  36
  37 using namespace std;
  38
  39
  40 void Plot_muons(){
  41
  42
  43   gStyle->SetOptStat(0);
  44   gStyle->SetOptDate(0);
  45   gStyle->SetOptFit(0111);
  46
  47   int ChComb=1;// 0=Same-charge, 1=Mixed-charge
  48   int kTbin_L=5, kTbin_H=6;
  49   int binKT3=1;//1-2
  50
  51   //TFile *_file0= new TFile("Results/PDC_12a17a_muons_R6.root","READ");
  52   //TFile *_file0= new TFile("Results/PDC_13b2_efix_p1_muons_R2.root","READ");
  53   TFile *_file0= new TFile("Results/PDC_12a17e_muons_R4.root","READ");
  54
  55   TList *MyList=(TList*)_file0->Get("MyList");
  56   _file0->Close();
  57
  58   TH1D *PionCandidates=(TH1D*)MyList->FindObject("fPionCandidates");
  59   PionCandidates->GetXaxis()->SetTitle("PDG code");
  60   //PionCandidates->Draw();
  61   //
  62   TH1D *MuonParentsPrimary=(TH1D*)MyList->FindObject("fMuonParents");
  63   MuonParentsPrimary->GetXaxis()->SetTitle("PDG code");
  64   MuonParentsPrimary->SetFillColor(1);
  65   //MuonParentsPrimary->Draw();
  66   //
  67   TH1D *MuonParentsSecondary=(TH1D*)MyList->FindObject("fSecondaryMuonParents");
  68   MuonParentsSecondary->GetXaxis()->SetTitle("PDG code");
  69   MuonParentsSecondary->SetFillColor(1);
  70   //MuonParentsSecondary->Draw();
  71   //
  72   // M0 R10-R6, M6 for R4, M17 for R2
  73   TH3D *PurityNum_3D = (TH3D*)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityNum");
  74   TH2D *PurityDen_2D = (TH2D*)MyList->FindObject("Explicit2_Charge1_1_Charge2_1_SC_0_M_6_ED_0_Term_1_PIDpurityDen");
  75   TH1D *PurityNum=PurityNum_3D->ProjectionX("PurityNum",kTbin_L,kTbin_H,1,20);
  76   double PurityNorm=PurityDen_2D->Integral(kTbin_L,kTbin_H,1,20);
  77   PurityNum->Scale(1/PurityNorm);
  78   char *namesAxis[15]={"e-e","e-mu","e-pi","e-k","e-p","mu-mu","mu-pi","mu-k","mu-p","pi-pi","pi-k","pi-p","k-k","k-p","p-p"};
  79   for(int i=1; i<=15; i++) PurityNum->GetXaxis()->SetBinLabel(i, namesAxis[i-1]);
  80   PurityNum->GetXaxis()->SetRange(1,15);
  81   PurityNum->GetYaxis()->SetTitle("Probability");
  82   PurityNum->Draw();
  83   //
  84   //
  85   TCanvas *can = new TCanvas("can", "can",800,0,800,800);// 11,53,700,500
  86   can->SetHighLightColor(2);
  87   gStyle->SetOptFit(0111);
  88   can->SetFillColor(10);//10
  89   can->SetBorderMode(0);
  90   can->SetBorderSize(2);
  91   can->SetFrameFillColor(0);
  92   can->SetFrameBorderMode(0);
  93   can->SetFrameBorderMode(0);
  94   can->cd();
  95   TPad *pad = new TPad("pad","pad",0.,0.,1.,1.);
  96   gPad->SetTickx();
  97   gPad->SetTicky();
  98   pad->SetGridx();
  99   pad->SetGridy();
 100   pad->SetTopMargin(0.02);//0.05
 101   pad->SetRightMargin(0.02);//3e-2
 102   pad->SetBottomMargin(0.1);//0.12
 103   pad->SetLeftMargin(0.1);
 104   pad->Draw();
 105   pad->cd();
 106   TLegend *legend = new TLegend(.5,.65, .9,.95,NULL,"brNDC");//.45 or .4 for x1
 107   legend->SetBorderSize(0);
 108   legend->SetFillColor(0);
 109   legend->SetTextFont(42);
 110   legend->SetTextSize(0.03);
 111   //
 112   TH3D *MuonSmearedNum2_3=(TH3D*)MyList->FindObject("fMuonContamSmearedNum2");
 113   TH3D *MuonSmearedDen2_3=(TH3D*)MyList->FindObject("fMuonContamSmearedDen2");
 114   TH3D *PionNum2_3=(TH3D*)MyList->FindObject("fMuonContamIdealNum2");
 115   TH3D *PionDen2_3=(TH3D*)MyList->FindObject("fMuonContamIdealDen2");
 116   TH3D *PionPionK2_3=(TH3D*)MyList->FindObject("fPionPionK2");
 117   //
 118   TH3D *MuonSmearedNum3_3=(TH3D*)MyList->FindObject("fMuonContamSmearedNum3");
 119   TH3D *MuonSmearedDen3_3=(TH3D*)MyList->FindObject("fMuonContamSmearedDen3");
 120   TH3D *PionNum3_3=(TH3D*)MyList->FindObject("fMuonContamIdealNum3");
 121   TH3D *PionDen3_3=(TH3D*)MyList->FindObject("fMuonContamIdealDen3");
 122   TH3D *PionPionK3_3=(TH3D*)MyList->FindObject("fPionPionK3");
 123   TH3D *MuonPionK3_3=(TH3D*)MyList->FindObject("fMuonPionK3");
 124
 125
 126   TH1D *MuonSmearedNum2[2];// SC/MC
 127   TH1D *MuonSmearedDen2[2];
 128   TH1D *PionNum2[2];
 129   TH1D *PionDen2[2];
 130   TH1D *PionPionK2[2];
 131   //
 132   TH1D *MuonSmearedNum3[2];// SC/MC
 133   TH1D *MuonSmearedDen3[2];
 134   TH1D *PionNum3[2];
 135   TH1D *PionDen3[2];
 136   TH1D *PionPionK3[2];
 137   //
 138   TH1D *C2muonSmeared[2];
 139   TH1D *C2pion[2];
 140   TH1D *C3muonSmeared[2];
 141   TH1D *C3pion[2];
 142   //
 143   for(int chtype=0; chtype<2; chtype++){
 144     TString *names[10];
 145     for(int i=0; i<10; i++) {names[i]=new TString("name_"); *names[i] += i; *names[i] += chtype;}
 146
 147     MuonSmearedNum2[chtype]=(TH1D*)MuonSmearedNum2_3->ProjectionZ(names[0]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
 148     MuonSmearedDen2[chtype]=(TH1D*)MuonSmearedDen2_3->ProjectionZ(names[1]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
 149     PionNum2[chtype]=(TH1D*)PionNum2_3->ProjectionZ(names[2]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
 150     PionDen2[chtype]=(TH1D*)PionDen2_3->ProjectionZ(names[3]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
 151     PionPionK2[chtype]=(TH1D*)PionPionK2_3->ProjectionZ(names[4]->Data(),chtype+1,chtype+1,kTbin_L,kTbin_H);
 152     PionPionK2[chtype]->Divide(PionDen2[chtype]);
 153     ////////////////
 154     MuonSmearedNum3[chtype]=(TH1D*)MuonSmearedNum3_3->ProjectionZ(names[5]->Data(),chtype+1,chtype+1,binKT3,binKT3);
 155     MuonSmearedDen3[chtype]=(TH1D*)MuonSmearedDen3_3->ProjectionZ(names[6]->Data(),chtype+1,chtype+1,binKT3,binKT3);
 156     PionNum3[chtype]=(TH1D*)PionNum3_3->ProjectionZ(names[7]->Data(),chtype+1,chtype+1,binKT3,binKT3);
 157     PionDen3[chtype]=(TH1D*)PionDen3_3->ProjectionZ(names[8]->Data(),chtype+1,chtype+1,binKT3,binKT3);
 158     PionPionK3[chtype]=(TH1D*)PionPionK3_3->ProjectionZ(names[9]->Data(),chtype+1,chtype+1,binKT3,binKT3);
 159     PionPionK3[chtype]->Divide(PionDen3[chtype]);
 160     //
 161     C2muonSmeared[chtype]=(TH1D*)MuonSmearedNum2[chtype]->Clone();
 162     C2pion[chtype]=(TH1D*)PionNum2[chtype]->Clone();
 163     C2muonSmeared[chtype]->Divide(MuonSmearedDen2[chtype]);
 164     C2pion[chtype]->Divide(PionDen2[chtype]);
 165     //
 166     C3muonSmeared[chtype]=(TH1D*)MuonSmearedNum3[chtype]->Clone();
 167     C3pion[chtype]=(TH1D*)PionNum3[chtype]->Clone();
 168     C3muonSmeared[chtype]->Divide(MuonSmearedDen3[chtype]);
 169     C3pion[chtype]->Divide(PionDen3[chtype]);
 170     //
 171     //
 172     C2pion[chtype]->SetLineColor(4);
 173     C2muonSmeared[chtype]->SetLineColor(2);
 174     C2pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
 175     C2pion[chtype]->SetMinimum(0.98); C2pion[chtype]->SetMaximum(1.35);
 176     C2pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
 177     C2pion[chtype]->GetXaxis()->SetTitle("q_{inv} (GeV/c)");
 178     C2pion[chtype]->GetYaxis()->SetTitle("C_{2}K_{2}");
 179     //
 180     C3pion[chtype]->SetLineColor(4);
 181     C3muonSmeared[chtype]->SetLineColor(2);
 182     C3pion[chtype]->GetXaxis()->SetRangeUser(0,0.15);
 183     C3pion[chtype]->SetMinimum(0.90); C3pion[chtype]->SetMaximum(2.0);
 184     C3pion[chtype]->GetYaxis()->SetTitleOffset(1.5);
 185     C3pion[chtype]->GetXaxis()->SetTitle("Q_{3} (GeV/c)");
 186     C3pion[chtype]->GetYaxis()->SetTitle("C_{3}K_{3}");
 187   }
 188   //
 189   //
 190
 191   C2pion[ChComb]->Draw();
 192   C2muonSmeared[ChComb]->Draw("same");
 193   legend->AddEntry(C2pion[ChComb],"Input Pion-Pion, C_{2}^{#pi-#pi,QS}K_{2}^{#pi-#pi,QS}","l");
 194   legend->AddEntry(C2muonSmeared[ChComb],"Pion-Muon residual, C_{2}^{#mu-#pi,QS}K_{2}^{#mu-#pi,QS}","l");
 195   legend->Draw("same");
 196   //PionPionK2->Draw("same");
 197   //
 198   //
 199
 200   //
 201   //C3pion[ChComb]->Draw();
 202   //C3muonSmeared[ChComb]->Draw("same");
 203   //PionPionK3[0]->Draw("same");
 204   //legend->AddEntry(C3pion[ChComb],"Input Pion-Pion-Pion, C_{3}^{#pi-#pi-#pi,QS}K_{3}^{#pi-#pi-#pi,QS}","l");
 205   //legend->AddEntry(C3muonSmeared[ChComb],"Muon-Pion-Pion residual, C_{3}^{#mu-#pi-#pi,QS}K_{3}^{#mu-#pi-#pi,QS}","l");
 206   //legend->Draw("same");
 207
 208
 209   // corrections
 210   TFile *fout=new TFile("MuonCorrection_temp.root","RECREATE");
 211   TH1D *C2muonCorrection[2];
 212   C2muonCorrection[0] = new TH1D("C2muonCorrection_SC","",100,0,0.5);
 213   C2muonCorrection[1] = new TH1D("C2muonCorrection_MC","",100,0,0.5);
 214   TH1D *WeightmuonCorrection = new TH1D("WeightmuonCorrection","",100,0,0.5);
 215   TH1D *C3muonCorrection[2];
 216   C3muonCorrection[0] = new TH1D("C3muonCorrection_SC","",50,0,0.5);
 217   C3muonCorrection[1] = new TH1D("C3muonCorrection_MC","",50,0,0.5);
 218   //
 219   C2muonCorrection[0]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[0]->GetYaxis()->SetTitle("x_{2}");
 220   C2muonCorrection[1]->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); C2muonCorrection[1]->GetYaxis()->SetTitle("x_{2}");
 221   C3muonCorrection[0]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[0]->GetYaxis()->SetTitle("x_{3}");
 222   C3muonCorrection[1]->GetXaxis()->SetTitle("Q_{3} (GeV/c)"); C3muonCorrection[1]->GetYaxis()->SetTitle("x_{3}");
 223   WeightmuonCorrection->GetXaxis()->SetTitle("q_{inv} (GeV/c)"); WeightmuonCorrection->GetYaxis()->SetTitle("x_{2}^{w}");
 224
 225   // 0.944 and 0.959
 226   float GoodPairFraction=1-0.93*(1-PurityNum->GetBinContent(10));
 227   cout<<"Pion Pair Purity = "<<PurityNum->GetBinContent(10)<<endl;
 228   cout<<"Effective Pion Pair Purity = "<<GoodPairFraction<<endl;
 229   float pionPurity=pow(GoodPairFraction,0.5);
 230   float muonPurity=1-pionPurity;
 231   for(int chtype=0; chtype<2; chtype++){
 232     for(int bin=1; bin<=100; bin++){
 233
 234       bool emptybin2=kFALSE, emptybin3=kFALSE;
 235       if(PionPionK2[chtype]->GetBinContent(bin)==0) {PionPionK2[chtype]->SetBinContent(bin, 1.00001);}
 236       if(PionPionK3[chtype]->GetBinContent(bin)==0) {PionPionK3[chtype]->SetBinContent(bin, 1.00001);}
 237       if(bin > C2pion[chtype]->GetNbinsX()) emptybin2=kTRUE;
 238       if(bin > C3pion[chtype]->GetNbinsX()) emptybin3=kTRUE;
 239
 240       double value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin);
 241       double den = (GoodPairFraction*C2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)*C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin);
 242       if(den > 0 && !emptybin2) C2muonCorrection[chtype]->SetBinContent(bin,value/den);
 243       else C2muonCorrection[chtype]->SetBinContent(bin, 1);
 244       //
 245       if(chtype==0){
 246         value = C2pion[chtype]->GetBinContent(bin)/PionPionK2[chtype]->GetBinContent(bin) - 1.0;
 247         den = ((GoodPairFraction*C2pion[chtype]->GetBinContent(bin) + (1-GoodPairFraction)*C2muonSmeared[chtype]->GetBinContent(bin))/PionPionK2[chtype]->GetBinContent(bin)) - 1.0;
 248         if(den > 0 && !emptybin2) WeightmuonCorrection->SetBinContent(bin,value/den);
 249       }
 250       //
 251       value = C3pion[chtype]->GetBinContent(bin)/PionPionK3[chtype]->GetBinContent(bin);
 252       den = (pow(pionPurity,3)*C3pion[chtype]->GetBinContent(bin) + (3*pow(pionPurity,2)*muonPurity)*C3muonSmeared[chtype]->GetBinContent(bin))/PionPionK3[chtype]->GetBinContent(bin);
 253       if(den > 0 && !emptybin3) C3muonCorrection[chtype]->SetBinContent(bin,value/den);
 254       else C3muonCorrection[chtype]->SetBinContent(bin, 1);
 255     }
 256     //C2muonCorrection[chtype]->SetBinContent(1, C2muonCorrection[chtype]->GetBinContent(2));
 257     //C3muonCorrection[chtype]->SetBinContent(1, C3muonCorrection[chtype]->GetBinContent(2));
 258     C2muonCorrection[chtype]->Write();
 259     C3muonCorrection[chtype]->Write();
 260     if(chtype==0) {
 261       //WeightmuonCorrection->SetBinContent(1, WeightmuonCorrection->GetBinContent(2));
 262       WeightmuonCorrection->Write();
 263     }
 264   }
 265
 266   //
 267   //C3muonCorrection[0]->SetMinimum(0.99);
 268   //C3muonCorrection[0]->SetMaximum(1.05);
 269   //C3muonCorrection[0]->GetYaxis()->SetTitleOffset(1.3);
 270   //C3muonCorrection[0]->Draw();
 271   //WeightmuonCorrection->GetYaxis()->SetTitleOffset(1.3);
 272   //WeightmuonCorrection->Draw();
 273   fout->Close();
 274
 275
 276 }