PWGLF/SPECTRA/PiKaPr/TestAOD/corrusingfitsigma.C

   1 //This macro calculates the correction factor for the difference of the n sigma TOF distributions in data and mc
   2
   3
   4 class AliSpectraBothHistoManager;
   5 AliSpectraBothHistoManager* managerdata=0x0;
   6 AliSpectraBothHistoManager* managermc=0x0;
   7
   8 Double_t Gaus_plus_tail(Double_t* x, Double_t* par )
   9 {
  10
  11         Double_t mean=par[0];
  12         Double_t rms=par[1];
  13         Double_t c=par[2];
  14         Double_t slope=par[3]/par[1];//df/dx continues
  15         Double_t cut=par[3];
  16
  17
  18
  19         Double_t one_over_sqrt_2pi=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
  20         Double_t returnvalue=0.0;
  21
  22         Double_t n=0.5*(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/(TMath::Abs(rms)*slope);
  23         if (x[0]<mean+cut*rms)
  24                 returnvalue=TMath::Exp(-1.0*(x[0]-mean)*(x[0]-mean)/(2.0*rms*rms))*one_over_sqrt_2pi/TMath::Abs(rms);
  25         else
  26                 returnvalue=TMath::Exp(slope*(mean+cut*rms-x[0]))*TMath::Exp(-cut*cut*0.5)*one_over_sqrt_2pi/TMath::Abs(rms);
  27
  28         return c*returnvalue/n;
  29 }
  30
  31 //Double_t background(Double_t x, Double_t scale, Double_t slope, Double_t startpoint)
  32 Double_t background(Double_t* x,Double_t* par)
  33 {
  34         //at start ponit it is eqaul to scale
  35         Double_t scale=par[0];
  36         Double_t slope=par[1];
  37         Double_t startpoint=par[2];
  38         return scale*TMath::Exp((startpoint-x[0])*slope);
  39 }
  40
  41
  42 Double_t Gaus_on_background(Double_t* x,Double_t* par)
  43 {
  44         //par[0] peak position
  45         //par[1] sigma gaus
  46         //par[2] Normalization of signal
  47         //par[3] cut
  48         //par[4] scale of bg
  49         //par[5] slope of bg
  50         //par[6] start point of bg const
  51         //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<par[6]<<par[7]<<" "<<x[0]<< endl;
  52         Double_t returnvalue=0.0;
  53          returnvalue+=Gaus_plus_tail(x,par);
  54         //returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);
  55         returnvalue+=background(x,&(par[4]));
  56         return returnvalue;
  57 }
  58
  59 Int_t   OpenFile(TString nameFile,TString outputname,Int_t mg1)
  60 {
  61
  62
  63         TFile *file = TFile::Open(nameFile.Data());
  64         if(!file)
  65         {
  66                 cout<<"no file"<<endl;
  67                 return -1;
  68         }
  69         TDirectoryFile *dir=(TDirectoryFile*)file->Get(outputname.Data());
  70         if(!dir)
  71         {
  72                 cout<<"no dir "<<outputname.Data()<<endl;
  73                 return -1;
  74         }
  75         if(!dir->Get("SpectraHistos"))
  76                 return -1;
  77         if(mg1==1)
  78         {
  79                 managerdata= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
  80                 return managerdata->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();
  81
  82         }
  83         else
  84         {
  85                 managermc= (AliSpectraBothHistoManager*) dir->Get("SpectraHistos");
  86                 return managermc->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->GetEntries();
  87         }
  88         return -1;
  89 }
  90
  91
  92 void corrusingfitsigma(TString nameFile1,TString outputname1,TString nameFile2,TString outputname2, Float_t nsigmacut=3.0, Float_t minpt=0.6,Float_t maxpt=1.4, Bool_t saveoption=kFALSE)
  93 {
  94         ofstream oufiletxt("TOFcorrPID.txt");
  95         Float_t nsigmacut=TMath::Abs(nsigmacut);
  96         TString pidmethods[3]={"TPC","TOF","TPCTOF"};
  97         TString Particle[]={"Pion","Kaon","Proton"};
  98         gStyle->SetOptStat(0);
  99         TH1::AddDirectory(kFALSE);
 100         gSystem->Load("libCore.so");
 101         gSystem->Load("libPhysics.so");
 102         gSystem->Load("libTree");
 103         gSystem->Load("libMatrix");
 104         gSystem->Load("libSTEERBase");
 105         gSystem->Load("libESD");
 106         gSystem->Load("libAOD");
 107         gSystem->Load("libANALYSIS");
 108         gSystem->Load("libOADB");
 109         gSystem->Load("libANALYSISalice");
 110         gSystem->Load("libTender");
 111         gSystem->Load("libCORRFW");
 112         gSystem->Load("libPWGTools");
 113         gSystem->Load("libPWGLFspectra");
 114         Int_t ne1= OpenFile(nameFile1,outputname1,1);
 115         Int_t ne2= OpenFile(nameFile2,outputname2,2);
 116         Int_t imethod=1;
 117                 oufiletxt<<"file 1 "<<nameFile1.Data()<<endl;
 118         oufiletxt<<"outputname 1 "<<outputname1.Data()<<endl;
 119         oufiletxt<<"file 2 "<<nameFile2.Data()<<endl;
 120         oufiletxt<<"outputname 2 "<<outputname2.Data()<<endl;
 121         oufiletxt<<"nsigma "<<nsigmacut<<endl;
 122         oufiletxt<<"minpt "<<minpt<<endl;
 123         oufiletxt<<"maxpt "<<maxpt<<endl;
 124
 125         TList* lout=new TList();
 126
 127         for(Int_t ipart=0;ipart<3;ipart++)
 128         {
 129                 if(!managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
 130                         continue;
 131                 TH2F *nsig_data = (TH2F*)((TH2F*)managerdata->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
 132                 if(!nsig_data)
 133                         continue;
 134
 135                 if(!managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))
 136                         continue;
 137                 TH2F *nsig_mc = (TH2F*)((TH2F*)managermc->GetNSigHistogram(Form("hHistNSig%sPt%s",Particle[ipart].Data(),pidmethods[imethod].Data())))->Clone();
 138                 if(!nsig_mc)
 139                         continue;
 140                 Int_t firstbin=nsig_data->GetXaxis()->FindBin(minpt);
 141                 Int_t lastbin=nsig_data->GetXaxis()->FindBin(maxpt);
 142                 Float_t ptlow=nsig_data->GetXaxis()->GetBinLowEdge(firstbin);
 143                 Float_t ptup=nsig_data->GetXaxis()->GetBinUpEdge(lastbin);
 144
 145                 TH1D* hist1data= new TH1D(Form("datahist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 146                 TH1D* hist2data= new TH1D(Form("datahist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 147                 TH1D* sigmadata= new TH1D(Form("sigmadata%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 148
 149                 hist1data->SetLineColor(kBlack);
 150                 hist2data->SetLineColor(kBlack);
 151                 sigmadata->SetMarkerColor(kBlack);
 152                 sigmadata->SetMarkerStyle(22);
 153                 hist1data->SetLineStyle(1);
 154                 hist2data->SetLineStyle(2);
 155                 hist1data->SetLineWidth(2);
 156                 hist2data->SetLineWidth(2);
 157
 158
 159                 TH1D* hist1mc= new TH1D(Form("mchist3sigmafun4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 160                 TH1D* hist2mc= new TH1D(Form("mchist3sigmahist4sigma%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 161                 TH1D* sigmamc= new TH1D(Form("sigmamc%s",Particle[ipart].Data()),"",lastbin-firstbin+1,ptlow,ptup);
 162
 163                 hist1mc->SetLineColor(kRed);
 164                 hist2mc->SetLineColor(kRed);
 165                 sigmamc->SetMarkerColor(kRed);
 166                 sigmamc->SetMarkerStyle(23);
 167                 sigmamc->SetLineColor(kRed);
 168                 hist1mc->SetLineStyle(1);
 169                 hist2mc->SetLineStyle(2);
 170                 hist1mc->SetLineWidth(2);
 171                 hist2mc->SetLineWidth(2);
 172
 173
 174
 175                 for(int ibin=firstbin;ibin<lastbin;ibin++)
 176                 {
 177
 178                         TH1F *nsig_data_Proj1=(TH1F*)nsig_data->ProjectionY(Form("%s%sdata[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
 179                         nsig_data_Proj1->Sumw2();
 180                         nsig_data_Proj1->GetXaxis()->SetRangeUser(-6,6);
 181                         nsig_data_Proj1->SetLineColor(kBlack);
 182                         TH1F *nsig_mc_Proj1=(TH1F*)nsig_mc->ProjectionY(Form("%s%smc[%.2f,%.2f]",Particle[ipart].Data(),pidmethods[imethod].Data(),nsig_mc->GetXaxis()->GetBinLowEdge(ibin),nsig_mc->GetXaxis()->GetBinUpEdge(ibin)),ibin,ibin);
 183                         nsig_mc_Proj1->Sumw2();
 184                         nsig_mc_Proj1->GetXaxis()->SetRangeUser(-6,6);
 185                         nsig_mc_Proj1->SetLineColor(kRed);
 186
 187                         TF1* gfundata=new TF1("fun_for_fitdata",Gaus_on_background,-5,5,7);
 188                         gfundata->SetParameter(0,0.0);
 189                         gfundata->SetParameter(1,1.0);
 190                         gfundata->SetParameter(2,nsig_data_Proj1->GetBinContent(nsig_data_Proj1->GetXaxis()->FindBin(0.0)));
 191                         gfundata->FixParameter(3,1.0);
 192                         gfundata->FixParameter(4,0.0);
 193                         gfundata->FixParameter(5,1.0);
 194                         gfundata->FixParameter(6,-5);
 195                         gfundata->SetLineColor(kBlack);
 196                         nsig_data_Proj1->Fit(gfundata,"SR+","E");
 197                         Float_t sigmadatav=gfundata->GetParameter(1);
 198                         Float_t meandatav=gfundata->GetParameter(0);
 199                         sigmadata->SetBinContent(ibin-firstbin+1,gfundata->GetParameter(1));
 200                         sigmadata->SetBinError(ibin-firstbin+1,gfundata->GetParError(1));
 201
 202
 203                         TF1* gfunmc=new TF1("fun_for_fitmc",Gaus_on_background,-5,5,7);
 204                         gfunmc->SetParameter(0,0.0);
 205                         gfunmc->SetParameter(1,1.0);
 206                         gfunmc->SetParameter(2,nsig_mc_Proj1->GetBinContent(nsig_mc_Proj1->GetXaxis()->FindBin(0.0)));
 207                         gfunmc->FixParameter(3,1.0);
 208                         gfunmc->FixParameter(4,0.0);
 209                         gfunmc->FixParameter(5,1.0);
 210                         gfunmc->FixParameter(6,-5);
 211                         gfunmc->SetLineColor(kRed);
 212                         nsig_mc_Proj1->Fit(gfunmc,"SR+","E");
 213                         sigmamc->SetBinContent(ibin-firstbin+1,gfunmc->GetParameter(1));
 214                         sigmamc->SetBinError(ibin-firstbin+1,gfunmc->GetParError(1));
 215                         Float_t sigmamcv=gfunmc->GetParameter(1);
 216                         Float_t meanmcv=gfundata->GetParameter(0);
 217
 218
 219                         hist1data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));
 220                         Float_t width=nsig_data_Proj1->GetXaxis()->GetBinWidth(-1.0*nsigmacut);
 221                         Int_t lowbind=nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav);
 222                         Int_t upbind=nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav);
 223                         Float_t lowbinpartd=TMath::Abs(nsig_data_Proj1->GetXaxis()->GetBinLowEdge(lowbind)+1.0*nsigmacut*sigmadatav-meandatav);
 224                         Float_t upbinpartd=TMath::Abs(nsigmacut*sigmadatav+meandatav-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbind));
 225
 226
 227
 228                         Float_t corrd=(lowbinpartd*nsig_data_Proj1->GetBinContent(lowbind)+upbinpartd*nsig_data_Proj1->GetBinContent(upbind))/width;
 229
 230                         hist2data->SetBinContent(ibin-firstbin+1,nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))-corrd);
 231                         cout<<corrd/nsig_data_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmadatav+meandatav),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut*sigmadatav+meandatav))<<endl;
 232
 233                         hist1mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut),nsig_data_Proj1->GetXaxis()->FindBin(nsigmacut-0.0001)));
 234                         Int_t lowbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv);
 235                         Int_t upbinmc=nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv);
 236                         Float_t lowbinpartmc=TMath::Abs(nsig_mc_Proj1->GetXaxis()->GetBinLowEdge(lowbinmc)+1.0*nsigmacut*sigmamcv-meanmcv);
 237                         Float_t upbinpartmc=TMath::Abs(nsigmacut*sigmamcv+meanmcv-nsig_data_Proj1->GetXaxis()->GetBinUpEdge(upbinmc));
 238                         Float_t corrmc=(lowbinpartmc*nsig_mc_Proj1->GetBinContent(lowbinmc)+upbinpartmc*nsig_mc_Proj1->GetBinContent(upbinmc))/width;
 239                         cout<<corrmc/nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))<<" aaa "<<endl;
 240
 241                         hist2mc->SetBinContent(ibin-firstbin+1,nsig_mc_Proj1->Integral(nsig_data_Proj1->GetXaxis()->FindBin(-1.0*nsigmacut*sigmamcv+meanmcv),nsig_mc_Proj1->GetXaxis()->FindBin(nsigmacut*sigmamcv+meanmcv))-corrmc);
 242                         TCanvas* ctmp=new TCanvas(Form("c%d%s",ibin,Particle[ipart].Data()),Form("c%d%s",ibin,Particle[ipart].Data()),800,800);
 243                         ctmp->cd()->SetLogy();
 244
 245                         nsig_data_Proj1->GetXaxis()->SetTitle("#sigma_{TOF}");
 246                         nsig_data_Proj1->GetYaxis()->SetTitle("N_{entries}");
 247                         nsig_data_Proj1->SetTitle(Form("PionTOF[%.2f,%.2f]",nsig_data->GetXaxis()->GetBinLowEdge(ibin),nsig_data->GetXaxis()->GetBinUpEdge(ibin)));
 248                         nsig_data_Proj1->Draw("E1");
 249                         nsig_mc_Proj1->Draw("E1same");
 250                         if(saveoption)
 251                                 ctmp->SaveAs(Form("c%d%s.png",ibin,Particle[ipart].Data()));
 252                         delete ctmp;
 253                         delete gfundata;
 254                         delete gfunmc;
 255                         delete nsig_data_Proj1;
 256                         delete nsig_mc_Proj1;
 257
 258                 }
 259
 260                         hist1data->Sumw2();
 261                         hist1data->Print("all");
 262                         hist1data->Sumw2();
 263                         hist1mc->Sumw2();
 264                         hist2mc->Sumw2();
 265
 266                 hist2data->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 267                 hist2data->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
 268                 hist2mc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 269                 hist2mc->GetYaxis()->SetTitle("N(-3fit,3fit)/N(-3,3)");
 270                 hist2data->Divide(hist2data,hist1data,1,1,"B");
 271                 hist2mc->Divide(hist2mc,hist1mc,1,1,"B");
 272                 lout->Add(hist2data);
 273                 lout->Add(hist2mc);
 274
 275                 sigmadata->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 276                 sigmadata->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
 277                 sigmamc->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 278                 sigmamc->GetYaxis()->SetTitle("#sigma_{TOF} fitted");
 279
 280
 281
 282                 lout->Add(sigmadata);
 283                 lout->Add(sigmamc);
 284
 285                 TH1D* corr=(TH1D*)hist2data->Clone(Form("corr%s",Particle[ipart].Data()));
 286                 corr->Divide(hist2mc);
 287                 corr->SetLineColor(kBlue);
 288                 corr->GetYaxis()->SetTitle("correction");
 289                 corr->Fit("pol0","R0","",minpt,maxpt);
 290
 291                 oufiletxt<<Particle[ipart].Data()<<" "<<((TF1*)corr->GetListOfFunctions()->At(0))->GetParameter(0)<<endl;
 292                 lout->Add(corr);
 293                 TCanvas* corrcanvas=new TCanvas(Form("corrcanvas%s",Particle[ipart].Data()), Form("corrcanvas%s",Particle[ipart].Data()),600,600);
 294                 corrcanvas->SetMargin(0.12,0.05,0.15,0.1);
 295                 corr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 296                 corr->GetYaxis()->SetTitle("TOF PID mismatch corr.");
 297                 corr->GetYaxis()->SetTitleOffset(1.5*corr->GetYaxis()->GetTitleOffset());
 298                 corr->GetYaxis()->SetRangeUser(0.98,1.04);
 299
 300                 corrcanvas->cd();
 301                 corr->Draw("E1");
 302                 corrcanvas->SaveAs(Form("Corr%s.eps",Particle[ipart].Data()));
 303                 corrcanvas->SaveAs(Form("Corr%s.png",Particle[ipart].Data()));
 304
 305         }
 306         TFile* fileout=TFile::Open("TOFcorrPID.root","Recreate");
 307
 308         lout->Write("output",TObject::kSingleKey);
 309         fileout->Close();
 310 }