+++ /dev/null
-const Double_t kMean=0.135 ; //Approximate peak position to facilitate error estimate
-
-//-----------------------------------------------------------------------------
-void MakeMmixPi0(const Int_t centrality=0,
- const char* cModule="")
-{
- //---------------------------------------------------------------------------
- // This macro processes PWGPP QA output of the analysis task PHOSPbPbQA
- // (see analysis code in the class AliAnalysisTaskPHOSPbPbQA).
- // It fits Real/Mixed ratio, normalize Mixed and subtract it from Real
- // Arguments:
- // centrality: 0 or 1 for centralities 0-20% or 20-100%
- // cModule: string with the PHOS module to analyze: "", "SM1", "SM2" or "SM3"
- //-
- // Yuri Kharlov. 19.11.2011
- //---------------------------------------------------------------------------
-/* $Id$ */
-
- TFile * f = new TFile("PHOSPbPb_all.root") ;
- TList *histESD = (TList*) f->Get("PHOSPbPbQAResults");
- char key[125] ;
-
- const char* pid="All";
- const char* txtCent[] = {"0-20%","20-100%"};
- TH2F * hCentrality = (TH2F*)f->Get("hCenPHOS") ;
- TH1D * hCentrality1 = hCentrality->ProjectionX();
- TString inputKey;
- TString outputKey = Form("%s10_cent%d",pid,centrality);
-
- TH2F *h , *hAdd;
- TH2F *hm, *hmAdd;
- if (centrality == 0 || centrality == 1) {// centrality 0-20%, 20-100%
- inputKey = Form("hPi0%s%s_cen%d" ,pid,cModule,centrality);
- TH2F *h = (TH2F*)f->Get(inputKey) ;
- inputKey = Form("hMiPi0%s%s_cen%d",pid,cModule,centrality);
- TH2F *hm= (TH2F*)f->Get(inputKey) ;
- if (h==0) {
- printf("Missing histogram %s\n",inputKey);
- return;
- }
- }
- else {
- printf("Wrong centrality %d. Allowed values are 0,1,2,3,10.\n",centrality);
- return;
- }
-
- // Int_t nPadX=2,nPadY=2;
- // Int_t nbin=4 ;
- // Double_t xa[]={0.8, 2.0, 3.5, 6.0, 20.} ;
-
- Int_t nPadX=1,nPadY=1;
- Int_t nbin=1 ;
- Double_t xa[]={2.0, 20.} ;
-
- PPRstyle();
- gStyle->SetOptFit(111);
- gStyle->SetPadLeftMargin(0.14);
- gStyle->SetPadRightMargin(0.01);
- gStyle->SetPadTopMargin(0.05);
- gStyle->SetPadBottomMargin(0.08);
-
- TString txtModule;
- if (strcmp(cModule,"") ==0) txtModule="All PHOS modules";
- else if (strcmp(cModule,"SM1")==0) txtModule="PHOS module 4";
- else if (strcmp(cModule,"SM2")==0) txtModule="PHOS module 3";
- else if (strcmp(cModule,"SM3")==0) txtModule="PHOS module 2";
-
- TPaveText *label = new TPaveText(0.15,0.80,0.40,0.90,"NDC");
- label->SetFillColor(kWhite);
- label->SetBorderSize(1);
- label->AddText(Form("Centrality %s" ,txtCent[centrality]));
- label->AddText(Form("%s",txtModule.Data()));
-
- //Fit real only
- //Linear Bg
- char kkey[55];
- sprintf(kkey,outputKey.Data()) ;
- char key2[155];
- sprintf(key,"Mix%s",kkey) ;
- sprintf(key2,"%s_mr1",key) ;
- TH1D * mr1 = new TH1D(key2,"Mass",nbin,xa) ;
- sprintf(key2,"%s_sr1",key) ;
- TH1D * sr1 = new TH1D(key2,"Width",nbin,xa) ;
- sprintf(key2,"%s_ar1",key) ;
- TH1D * ar1 = new TH1D(key2,"a",nbin,xa) ;
- sprintf(key2,"%s_br1",key) ;
- TH1D * br1 = new TH1D(key2,"a",nbin,xa) ;
- sprintf(key2,"%s_yr1",key) ;
- TH1D * nr1 = new TH1D(key2,"Raw yield",nbin,xa) ;
- sprintf(key2,"%s_yr1int",key) ;
- TH1D * nr1int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;
-
- //Quadratic Bg
- sprintf(key2,"%s_mr2",key) ;
- TH1D * mr2 = new TH1D(key2,"Mass",nbin,xa) ;
- sprintf(key2,"%s_sr2",key) ;
- TH1D * sr2 = new TH1D(key2,"Width",nbin,xa) ;
- sprintf(key2,"%s_ar2",key) ;
- TH1D * ar2 = new TH1D(key2,"a",nbin,xa) ;
- sprintf(key2,"%s_br2",key) ;
- TH1D * br2 = new TH1D(key2,"a",nbin,xa) ;
- sprintf(key2,"%s_cr2",key) ;
- TH1D * cr2 = new TH1D(key2,"a",nbin,xa) ;
- sprintf(key2,"%s_yr2",key) ;
- TH1D * nr2 = new TH1D(key2,"Raw yield",nbin,xa) ;
- sprintf(key2,"%s_yr2int",key) ;
- TH1D * nr2int = new TH1D(key2,"Raw yield, integrated",nbin,xa) ;
-
- TF1 * fit1 = new TF1("fit",CB,0.,1.,6) ;
- fit1->SetParName(0,"A") ;
- fit1->SetParName(1,"m_{0}") ;
- fit1->SetParName(2,"#sigma") ;
- fit1->SetParName(3,"a_{0}") ;
- fit1->SetParName(4,"a_{1}") ;
- fit1->SetParName(5,"a_{2}") ;
- fit1->SetLineWidth(2) ;
- fit1->SetLineColor(2) ;
- TF1 * fgs = new TF1("gs",CBs,0.,1.,3) ;
- fgs->SetParName(0,"A") ;
- fgs->SetParName(1,"m_{0}") ;
- fgs->SetParName(2,"#sigma") ;
- fgs->SetLineColor(2) ;
- fgs->SetLineWidth(2) ;
-
- TF1 * fit2 = new TF1("fit2",CB2,0.,1.,7) ;
- fit2->SetParName(0,"A") ;
- fit2->SetParName(1,"m_{0}") ;
- fit2->SetParName(2,"#sigma") ;
- fit2->SetParName(3,"a_{0}") ;
- fit2->SetParName(4,"a_{1}") ;
- fit2->SetParName(5,"a_{2}") ;
- fit2->SetParName(6,"a_{3}") ;
- fit2->SetLineWidth(2) ;
- fit2->SetLineColor(4) ;
- fit2->SetLineStyle(2) ;
-
- TF1 * fbg1 = new TF1("bg1",BG1,0.,1.,3) ;
- TF1 * fbg2 = new TF1("bg2",BG2,0.,1.,4) ;
-
- TCanvas * c3 = new TCanvas("mggFit1_Signal","mggFitCB",10,10,1200,800) ;
- c3->Divide(nPadX,nPadY) ;
-
- TCanvas * c1 = new TCanvas("mggFit1","mggFit1",10,10,1200,800) ;
- c1->Divide(nPadX,nPadY) ;
- c1->cd(0) ;
-
- TCanvas *c4, *c2;
- TAxis * pta=h->GetYaxis() ;
- TAxis * ma=h->GetXaxis() ;
- for(Int_t i=1;i<=nbin;i++){
- if(i<=15)
- c1->cd(i) ;
- else{
- if(c2==0){
- c2 = new TCanvas("mggFit2","mggFit2",10,10,1200,800) ;
- c2->Divide(nPadX,nPadY) ;
- }
- c2->cd(i-16) ;
- }
- printf("\t%.1f<pt<%.1f GeV/c\n",xa[i-1],xa[i]);
- Int_t imin=pta->FindBin(xa[i-1]+0.0001);
- Int_t imax=pta->FindBin(xa[i]-0.0001) ;
- Double_t pt=(xa[i]+xa[i-1])/2. ;
- TH1D * hp = h->ProjectionX("re",imin,imax) ;
- hp->Sumw2() ;
- TH1D * hpm= hm->ProjectionX("mi",imin,imax) ;
- hpm->Sumw2() ;
- // if(i<=11){
- if(i<=10){
- hp ->Rebin(4) ;
- hpm->Rebin(4) ;
- }
- else{
- hp ->Rebin(5) ;
- hpm->Rebin(5) ;
- }
- for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hp ->GetBinContent(ib)==0)hp ->SetBinError(ib,1.);}
- for(Int_t ib=1; ib<=hp->GetNbinsX();ib++){if(hpm->GetBinContent(ib)==0)hpm->SetBinError(ib,1.);}
- TH1D * hpm2 = (TH1D*)hpm->Clone("Bg1") ;
- TH1D * hpcopy = (TH1D*)hp ->Clone("hpcopy") ;
- TH1D * hp2 = (TH1D*)hp ->Clone("hp2") ;
- hpcopy->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
- hp2 ->SetXTitle("M_{#gamma#gamma} (GeV/c^{2})");
- hpcopy->Divide(hpm) ;
- sprintf(key,"%3.1f<p_{T}<%3.1f GeV/c",xa[i-1],xa[i]) ;
- hpcopy->SetTitle(key) ;
- hpcopy->SetMarkerStyle(20) ;
- hpcopy->SetMarkerSize(0.7) ;
-
- fit1->SetParameters(0.001,0.136,0.005,0.0002,-0.002,0.0) ;
- fit1->SetParLimits(0,0.000,1.000) ;
- fit1->SetParLimits(1,0.125,0.145) ;
- fit1->SetParLimits(2,0.002,0.015) ;
-
- Double_t rangeMin=0.06;
- Double_t rangeMax=0.26;
-
- hpcopy->Fit(fit1,"Q" ,"",rangeMin,rangeMax) ;
- hpcopy->Fit(fit1,"MQ","",rangeMin,rangeMax) ;
-
- ar1->SetBinContent(i,fit1->GetParameter(3)) ;
- ar1->SetBinError (i,fit1->GetParError(3)) ;
- br1->SetBinContent(i,fit1->GetParameter(4)) ;
- br1->SetBinError (i,fit1->GetParError(4)) ;
-
- fit2->SetParameters(fit1->GetParameters()) ;
- fit2->SetParLimits(0,0.000,1.000) ;
- fit2->SetParLimits(1,0.125,0.145) ;
- fit2->SetParLimits(2,0.002,0.015) ;
-
- hpcopy->Fit(fit2,"+NQ","",rangeMin,rangeMax) ;
- hpcopy->Fit(fit2,"+MQ","",rangeMin,rangeMax) ;
-
- ar2->SetBinContent(i,fit2->GetParameter(3)) ;
- ar2->SetBinError (i,fit2->GetParError(3)) ;
- br2->SetBinContent(i,fit2->GetParameter(4)) ;
- br2->SetBinError (i,fit2->GetParError(4)) ;
- cr2->SetBinContent(i,fit2->GetParameter(5)) ;
- cr2->SetBinError (i,fit2->GetParError(5)) ;
- hpcopy->SetAxisRange(0.04,0.3,"X") ;
- hpcopy->Draw() ;
- label ->Draw();
- if(c2)
- c2->Update() ;
- else
- c1->Update() ;
-// if(getchar()=='q')return ;
-
- fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4),fit1->GetParameter(5));
- fbg2->SetParameters(fit2->GetParameter(3),fit2->GetParameter(4),fit2->GetParameter(5),
- fit2->GetParameter(6));
-
- Double_t intRangeMin = PeakPosition(pt,centrality)-3.*PeakWidth(pt,centrality) ;
- Double_t intRangeMax = PeakPosition(pt,centrality)+3.*PeakWidth(pt,centrality) ;
- // printf("pt=%f, %f < M < %f\n",pt,intRangeMin,intRangeMax);
- Int_t intBinMin = hp->GetXaxis()->FindBin(intRangeMin) ;
- Int_t intBinMax = hp->GetXaxis()->FindBin(intRangeMax) ;
- Double_t errStat = hpm->Integral(intBinMin,intBinMax);
-
- hpm ->Multiply(fbg1) ;
- hpm2->Multiply(fbg2) ;
- hp ->Add(hpm ,-1.) ;
- hp2 ->Add(hpm2,-1.) ;
-
- if(i<=15)
- c3->cd(i) ;
- else{
- if(c4==0){
- c4 = new TCanvas("mggFit2_Signal","mggFit2",10,10,1200,800) ;
- c4->Divide(nPadX,nPadY) ;
- }
- c4->cd(i-15) ;
- }
-
- Int_t binPi0 = hp->FindBin(kMean);
- fgs->SetParameters(hp->Integral(binPi0-1,binPi0+1)/3.,fit1->GetParameter(1),fit1->GetParameter(2)) ;
- fgs->SetParLimits(0,0.000,1.e+5) ;
- fgs->SetParLimits(1,0.120,0.145) ;
- fgs->SetParLimits(2,0.002,0.015) ;
- hp->Fit(fgs,"Q","",rangeMin,rangeMax) ;
- hp->SetMaximum(hp2->GetMaximum()*1.4) ;
- hp->SetMinimum(hp2->GetMinimum()*1.1) ;
- hp->SetMarkerStyle(20) ;
- hp->SetMarkerSize(0.7) ;
- mr1->SetBinContent(i,fgs->GetParameter(1)) ;
- mr1->SetBinError (i,fgs->GetParError(1) ) ;
- sr1->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;
- sr1->SetBinError (i,fgs->GetParError(2) ) ;
-
- Double_t y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;
- nr1->SetBinContent(i,y) ;
- Double_t ey=0 ;
- if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){
- Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;
- Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;
- ey=y*TMath::Sqrt(en*en+es*es) ;
- }
- nr1->SetBinError(i,ey) ;
-
- Double_t npiInt = hp->Integral(intBinMin,intBinMax) ;
- Double_t norm = fbg1->GetParameter(0) ;
- Double_t normErr= fbg1->GetParError(0) ;
- if(npiInt>0.){
- nr1int->SetBinContent(i,npiInt) ;
- nr1int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
- }
- hp2->SetAxisRange(0.04,0.3,"X") ;
- hp2->SetMaximum(hp2->GetMaximum()*1.4) ;
- hp2->SetMinimum(hp2->GetMinimum()*1.1) ;
- hp2->SetMarkerStyle(20) ;
- hp2->SetMarkerSize(0.7) ;
- hp2->Fit(fgs,"Q","",rangeMin,rangeMax) ;
- mr2->SetBinContent(i,fgs->GetParameter(1)) ;
- mr2->SetBinError (i,fgs->GetParError(1)) ;
- sr2->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;
- sr2->SetBinError (i,fgs->GetParError(2)) ;
- y=fgs->Integral(intRangeMin,intRangeMax)/hp->GetXaxis()->GetBinWidth(1) ;
- nr2->SetBinContent(i,y) ;
- ey=0 ;
- if(fgs->GetParameter(0)!=0. && fgs->GetParameter(2)!=0.){
- Double_t en=fgs->GetParError(0)/fgs->GetParameter(0) ;
- Double_t es=fgs->GetParError(2)/fgs->GetParameter(2) ;
- ey=y*TMath::Sqrt(en*en+es*es) ;
- }
- nr2->SetBinError(i,ey) ;
- npiInt=hp2->Integral(intBinMin,intBinMax) ;
- norm=fbg2->GetParameter(0) ;
- normErr=fbg2->GetParError(0) ;
- if(npiInt>0.){
- nr2int->SetBinContent(i,npiInt) ;
- nr2int->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
- }
- hp2->SetTitle(key) ;
- hp2->Draw() ;
- label ->Draw();
- if(c4)
- c4->Update() ;
- else
- c3->Update() ;
-
- delete hp ;
-// delete hp2 ;
-// delete hpcopy ;
- delete hpm ;
- delete hpm2 ;
- }
- char name[55] ;
-
- sprintf(name,"Pi0_ratio_cent%d%s.eps" ,centrality,cModule) ;
- c1->Print(name) ;
- sprintf(name,"Pi0_signal_cent%d%s.eps",centrality,cModule) ;
- c3->Print(name) ;
-
- //Normalize by the number of events
- Int_t cMin,cMax;
- if (centrality == 0) {
- cMin=1;
- cMax=20;
- }
- else if (centrality == 1) {
- cMin=21;
- cMax=80;
- }
- Double_t nevents = hCentrality1->Integral(cMin,cMax);
- nr1 ->Scale(1./nevents) ;
- nr1int->Scale(1./nevents) ;
- nr2 ->Scale(1./nevents) ;
- nr2int->Scale(1./nevents) ;
-
- printf("\t==============================\n");
- printf("\t| N events = %i8 |\n",nevents);
- printf("\t==============================\n");
-
- TFile fout("LHC10h_Pi0_FitResult.root","update");
- mr1->Write(0,TObject::kOverwrite) ;
- sr1->Write(0,TObject::kOverwrite) ;
- ar1->Write(0,TObject::kOverwrite) ;
- br1->Write(0,TObject::kOverwrite) ;
- nr1->Write(0,TObject::kOverwrite) ;
- nr1int->Write(0,TObject::kOverwrite) ;
- ar2->Write(0,TObject::kOverwrite) ;
- br2->Write(0,TObject::kOverwrite) ;
- cr2->Write(0,TObject::kOverwrite) ;
- mr2->Write(0,TObject::kOverwrite) ;
- sr2->Write(0,TObject::kOverwrite) ;
- nr2->Write(0,TObject::kOverwrite) ;
- nr2int->Write(0,TObject::kOverwrite) ;
- fout.Close() ;
-
-}
-
-//-----------------------------------------------------------------------------
-Double_t PeakPosition(Double_t pt, Int_t centrality = 0){
- //Fit to the measured peak position
- Double_t pi0Peak = 0.135;
- if (centrality == 0) pi0Peak=0.1413;
- else if (centrality == 1) pi0Peak=0.1380;
- else if (centrality == 2) pi0Peak=0.1367;
- else if (centrality == 3) pi0Peak=0.1359;
- return pi0Peak;
- // return 4.99292e-003*exp(-pt*9.32300e-001)+1.34944e-001 ;
-}
-//-----------------------------------------------------------------------------
-Double_t PeakWidth(Double_t pt, Int_t centrality = 0){
- //fit to the measured peak width
- Double_t pi0Sigma = 0.07;
- if (centrality == 0) pi0Sigma=0.0084;
- else if (centrality == 1) pi0Sigma=0.0071;
- else if (centrality == 2) pi0Sigma=0.0068;
- else if (centrality == 3) pi0Sigma=0.0064;
- return pi0Sigma;
-
- // Double_t a=0.0068 ;
- // Double_t b=0.0025 ;
- // Double_t c=0.000319 ;
- // return TMath::Sqrt(a*a+b*b/pt/pt+c*c*pt*pt) ;
-}
-
-//-----------------------------------------------------------------------------
-Double_t CB(Double_t * x, Double_t * par){
- //Parameterization of Real/Mixed ratio
- Double_t m=par[1] ;
- Double_t s=par[2] ;
- Double_t dx=(x[0]-m)/s ;
- Double_t dMpi0 = x[0]-kMean;
- Double_t peak = par[0] * TMath::Exp(-dx*dx/2.);
- Double_t bg = par[3] + par[4]*dMpi0 + par[5]*dMpi0*dMpi0;
- return peak+bg ;
-}
-
-//-----------------------------------------------------------------------------
-Double_t CB2(Double_t * x, Double_t * par){
- //Another parameterization of Real/Mixed ratio7TeV/README
- Double_t m=par[1] ;
- Double_t s=par[2] ;
- Double_t dx=(x[0]-m)/s ;
- Double_t dMpi0 = x[0]-kMean;
- Double_t peak = par[0] * TMath::Exp(-dx*dx/2.);
- Double_t bg = par[3] + par[4]*dMpi0 + par[5]*dMpi0*dMpi0 + par[6]*dMpi0*dMpi0*dMpi0;
- return peak+bg ;
-}
-//-----------------------------------------------------------------------------
-Double_t CBs(Double_t * x, Double_t * par){
- //Parameterizatin of signal
- Double_t m=par[1] ;
- Double_t s=par[2] ;
- Double_t dx=(x[0]-m)/s ;
- Double_t peak = par[0] * TMath::Exp(-dx*dx/2.);
- return peak ;
-}
-//-----------------------------------------------------------------------------
-Double_t BG1(Double_t * x, Double_t * par){
- //Normalizatino of Mixed
- return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean) ;
-}
-//-----------------------------------------------------------------------------
-Double_t BG2(Double_t * x, Double_t * par){
- //Another normalization of Mixed
- return par[0]+par[1]*(x[0]-kMean)+par[2]*(x[0]-kMean)*(x[0]-kMean)+par[3]*(x[0]-kMean)*(x[0]-kMean)*(x[0]-kMean) ;
-}
-
-
-//-----------------------------------------------------------------------------
-PPRstyle()
-{
-
- //////////////////////////////////////////////////////////////////////
- //
- // ROOT style macro for the TRD TDR
- //
- //////////////////////////////////////////////////////////////////////
-
- gStyle->SetPalette(1);
- gStyle->SetCanvasBorderMode(-1);
- gStyle->SetCanvasBorderSize(1);
- gStyle->SetCanvasColor(10);
-
- gStyle->SetFrameFillColor(10);
- gStyle->SetFrameBorderSize(1);
- gStyle->SetFrameBorderMode(-1);
- gStyle->SetFrameLineWidth(1.2);
- gStyle->SetFrameLineColor(1);
-
- gStyle->SetHistFillColor(0);
- gStyle->SetHistLineWidth(1);
- gStyle->SetHistLineColor(1);
-
- gStyle->SetPadColor(10);
- gStyle->SetPadBorderSize(1);
- gStyle->SetPadBorderMode(-1);
-
- gStyle->SetStatColor(10);
- gStyle->SetTitleColor(kBlack,"X");
- gStyle->SetTitleColor(kBlack,"Y");
-
- gStyle->SetLabelSize(0.04,"X");
- gStyle->SetLabelSize(0.04,"Y");
- gStyle->SetLabelSize(0.04,"Z");
- gStyle->SetTitleSize(0.04,"X");
- gStyle->SetTitleSize(0.04,"Y");
- gStyle->SetTitleSize(0.04,"Z");
- gStyle->SetTitleFont(42,"X");
- gStyle->SetTitleFont(42,"Y");
- gStyle->SetTitleFont(42,"X");
- gStyle->SetLabelFont(42,"X");
- gStyle->SetLabelFont(42,"Y");
- gStyle->SetLabelFont(42,"Z");
- gStyle->SetStatFont(42);
-
- gStyle->SetTitleOffset(1.0,"X");
- gStyle->SetTitleOffset(1.4,"Y");
-
- gStyle->SetFillColor(kWhite);
- gStyle->SetTitleFillColor(kWhite);
-
- gStyle->SetOptDate(0);
- gStyle->SetOptTitle(1);
- gStyle->SetOptStat(0);
- gStyle->SetOptFit(0);
-
-}
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff