PWGGA/PHOSTasks/PHOSCalib/CompWidth.C

   1 void CompWidth(){
   2
   3   Int_t nIter=18;
   4   TFile *f[20] ;
   5   TH3F * hReM1[20] ;
   6   TH3F * hMiM1[20] ;
   7   TH3F * hReM2[20] ;
   8   TH3F * hMiM2[20] ;
   9   TH3F * hReM3[20] ;
  10   TH3F * hMiM3[20] ;
  11   char key[155] ;
  12
  13   for(Int_t i=0;i<nIter;i++){
  14     sprintf(key,"histos_pass%d.root",i) ;
  15     f[i] = TFile::Open(key) ;
  16     hReM1[i] = (TH3F*)f[i]->Get("hMggDispM1") ;
  17     sprintf(key,"hMggDispM1_iter%d",i) ;
  18     hReM1[i]->SetName(key) ;
  19     hMiM1[i] = (TH3F*)f[i]->Get("hMiMggDispM1") ;
  20     sprintf(key,"hMiMggDispM1_iter%d",i) ;
  21     hMiM1[i]->SetName(key) ;
  22
  23     hReM2[i] = (TH3F*)f[i]->Get("hMggDispM2") ;
  24     sprintf(key,"hMggDispM2_iter%d",i) ;
  25     hReM2[i]->SetName(key) ;
  26     hMiM2[i] = (TH3F*)f[i]->Get("hMiMggDispM2") ;
  27     sprintf(key,"hMiMggDispM2_iter%d",i) ;
  28     hMiM2[i]->SetName(key) ;
  29
  30     hReM3[i] = (TH3F*)f[i]->Get("hMggDispM3") ;
  31     sprintf(key,"hMggDispM3_iter%d",i) ;
  32     hReM3[i]->SetName(key) ;
  33     hMiM3[i] = (TH3F*)f[i]->Get("hMiMggDispM3") ;
  34     sprintf(key,"hMiMggDispM3_iter%d",i) ;
  35     hMiM3[i]->SetName(key) ;
  36   }
  37
  38   TH1D * hW1 = new TH1D("WidthM1","Module 1",nIter,-0.5,nIter-0.5) ;
  39   TH1D * hW2 = new TH1D("WidthM2","Module 1",nIter,-0.5,nIter-0.5) ;
  40   TH1D * hW3 = new TH1D("WidthM3","Module 1",nIter,-0.5,nIter-0.5) ;
  41
  42 //  TF1 * gs = new TF1("gs","[0]*exp(-(x-[1])*(x-[1])/2./[2]/[2])+[3]+[4]*x",0.,1.) ;
  43   TF1 * gs = new TF1("gs","[0]*exp(-(x-[1])*(x-[1])/2./[2]/[2])+[3]",0.,1.) ;
  44   TCanvas * c = new TCanvas("c","c") ;
  45
  46   TH1D * hD[100] ;
  47
  48   for(Int_t i=0;i<nIter;i++){
  49     TH1D * tmp = hReM1[i]->ProjectionZ("a") ;
  50     TH1D * tmpMi = hMiM1[i]->ProjectionZ("b") ;
  51     tmp->Sumw2() ;
  52     //Normalize
  53     Double_t nMi=tmpMi->Integral(70,120) ;
  54     Double_t nRe=tmp->Integral(70,120) ;
  55     if(nMi>0.)
  56       tmpMi->Scale(nRe/nMi) ;
  57     tmp->Add(tmpMi,-1.) ;
  58     gs->SetParameters(3500.,0.135,0.008,0.,0.) ;
  59     //    gs->SetParLimits(0,0.,1.e+6) ;
  60     gs->SetParLimits(1,0.09,0.2) ;
  61     gs->SetParLimits(2,0.003,0.015) ;
  62
  63     c->cd() ;
  64 //    tmp->Fit(gs,"q","",0.05,0.220) ;
  65     tmp->GetXaxis()->SetRangeUser(0.,0.3) ;
  66     tmp->SetXTitle("m_{#gamma#gamma} (GeV/c^{2})") ;
  67     tmp->SetTitle("Module 2") ;
  68     tmp->Fit(gs,"","",0.125,0.145) ;
  69     hD[i]=(TH1D*)tmp->Clone(Form("Copy%d",i)) ;
  70 //    c->Update() ; if(getchar()=='q')return  ;
  71 //    if(i==12) return ;
  72     //    mass->SetBinContent(i,j,gs->GetParameter(1)) ;
  73     //  mass->SetBinError(i,j,gs->GetParError(1)) ;
  74     hW1->SetBinContent(i+1,gs->GetParameter(2)) ;
  75     hW1->SetBinError(i+1,gs->GetParError(2)) ;
  76     delete tmp ;
  77     delete tmpMi ;
  78
  79     tmp = hReM2[i]->ProjectionZ("a") ;
  80     tmpMi = hMiM2[i]->ProjectionZ("b") ;
  81     tmp->Sumw2() ;
  82     //Normalize
  83     nMi=tmpMi->Integral(70,120) ;
  84     nRe=tmp->Integral(70,120) ;
  85     if(nMi>0.)
  86       tmpMi->Scale(nRe/nMi) ;
  87     tmp->Add(tmpMi,-1.) ;
  88     gs->SetParameters(10000.,0.150,0.008,0.,0.) ;
  89     //    gs->SetParLimits(0,0.,100.) ;
  90     gs->SetParLimits(1,0.09,0.2) ;
  91     gs->SetParLimits(2,0.003,0.015) ;
  92
  93     tmp->Fit(gs,"q","",0.125,0.145) ;
  94     //    mass->SetBinContent(i,j,gs->GetParameter(1)) ;
  95     //  mass->SetBinError(i,j,gs->GetParError(1)) ;
  96     hW2->SetBinContent(i+1,gs->GetParameter(2)) ;
  97     hW2->SetBinError(i+1,gs->GetParError(2)) ;
  98     delete tmp ;
  99     delete tmpMi ;
 100
 101     tmp = hReM3[i]->ProjectionZ("a") ;
 102     tmpMi = hMiM3[i]->ProjectionZ("b") ;
 103     tmp->Sumw2() ;
 104     //Normalize
 105     Double_t nMi=tmpMi->Integral(70,120) ;
 106     Double_t nRe=tmp->Integral(70,120) ;
 107     if(nMi>0.)
 108       tmpMi->Scale(nRe/nMi) ;
 109     tmp->Add(tmpMi,-1.) ;
 110     gs->SetParameters(10000.,0.150,0.008,0.,0.) ;
 111     //    gs->SetParLimits(0,0.,100.) ;
 112     gs->SetParLimits(1,0.09,0.2) ;
 113     gs->SetParLimits(2,0.003,0.015) ;
 114
 115 //    tmp->Fit(gs,"q","",0.05,0.220) ;
 116     tmp->Fit(gs,"q","",0.125,0.145) ;
 117     //    mass->SetBinContent(i,j,gs->GetParameter(1)) ;
 118     //  mass->SetBinError(i,j,gs->GetParError(1)) ;
 119     hW3->SetBinContent(i+1,gs->GetParameter(2)) ;
 120     hW3->SetBinError(i+1,gs->GetParError(2)) ;
 121     delete tmp ;
 122     delete tmpMi ;
 123
 124   }
 125
 126   hW1->SetMarkerStyle(20) ;
 127   hW1->SetMarkerColor(2) ;
 128   hW2->SetMarkerStyle(21) ;
 129   hW2->SetMarkerColor(4) ;
 130   hW3->SetMarkerStyle(24) ;
 131   hW3->SetMarkerColor(8) ;
 132
 133   hW1->SetXTitle("Iteration") ;
 134   hW1->SetYTitle("#sigma (GeV/c^{2})") ;
 135   hW1->Draw() ;
 136   hW2->Draw("same") ;
 137   hW3->Draw("same") ;
 138
 139   TLegend * l = new TLegend(0.7,0.7,0.85,0.85) ;
 140   l->AddEntry(hW3,"Module 2","p") ;
 141   l->AddEntry(hW2,"Module 3","p") ;
 142   l->AddEntry(hW1,"Module 4","p") ;
 143   l->Draw() ;
 144
 145   TCanvas * cMinv = new TCanvas("cMinv") ;
 146   hD[0]->SetMarkerStyle(20) ;
 147   hD[0]->SetMarkerSize(0.8) ;
 148   hD[0]->Draw() ;
 149   for(Int_t i=1; i<nIter; i++){
 150    hD[i]->SetMarkerStyle(20+i) ;
 151    hD[i]->SetMarkerSize(0.8) ;
 152    hD[i]->Draw("same") ;
 153   }
 154
 155 }