TPC/TPCupgrade/macros/fitResolution.C

   1 void EmpiricalGEMQFit(const char * input ="GEMScansIKF.root" ){
   2   //
   3   // Empirical fit of the relative Q resolution for the iron source with 4 GEM layers
   4   // Fit:
   5   // Assumption :
   6   // 1.) Effective electron transparency proportiaonal to the effective ion transparency
   7   // 2.) RMS of the effecitive gain can be expressed as linear function of U1/U3,U2/U4 and U3/U4
   8   //
   9   // Agreement with data within the expect error estimate -relative agreement
  10   // RMS (sigma_{meas}/sigma_{fit}) ~ 3%
  11   //
  12   TFile *fgem = TFile::Open(input);
  13   TTree * tree= (TTree*)fgem->Get("NeCO2N2");
  14   tree->SetMarkerStyle(25);
  15   TStatToolkit toolkit;
  16   Double_t chi20=0;
  17   Int_t    npoints=0;
  18   Int_t npointsMax=10000;
  19   TVectorD param0,param1,param2,param3;
  20   TMatrixD covar0,covar1,covar2,covar3;
  21
  22   tree->SetAlias("UGEMA","(UGEM1+UGEM2+UGEM3+UGEM4)");
  23   TString fstringFast="";
  24   fstringFast+="1/IB++";                // fraction of the
  25   fstringFast+="(UGEM1/UGEM4)/IB++";    // fraction of the gain
  26   fstringFast+="(UGEM2/UGEM4)/IB++";    // fraction of the gain
  27   fstringFast+="(UGEM3/UGEM4)/IB++";    // fraction of the gain
  28   //
  29   TCut cutFit="ET2Scan+ET3Scan==0";
  30   TString *strResolFit = TStatToolkit::FitPlane(tree,"Sigma^2:Sigma^2", fstringFast.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 0);
  31   strResolFit->Tokenize("++")->Print();
  32   tree->SetAlias("fitSigma2",strResolFit->Data());
  33   //
  34   gStyle->SetOptTitle(0);
  35   //
  36   TCanvas *canvas = new TCanvas("canvasEmp","canvasEmp",600,500);
  37   canvas->SetBottomMargin(0.15);
  38   canvas->SetRightMargin(0.1);
  39   canvas->SetTicks(1,1);
  40   tree->SetMarkerSize(0.7);
  41   {
  42     tree->Draw("Sigma:sqrt(fitSigma2):sqrt(1/IB)",cutFit,"colz");
  43     TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
  44     htemp->GetXaxis()->SetTitle("#sigma_{fit}(%)");
  45     htemp->GetYaxis()->SetTitle("#sigma_{meas}(%)");
  46     htemp->GetZaxis()->SetTitle("#sqrt{IBF}");
  47     htemp->SetTitle("Fe resolution");
  48     htemp->GetXaxis()->SetLimits(8,20);
  49     htemp->GetYaxis()->SetLimits(8,20);
  50     htemp->Draw("colz");
  51     TLatex latex;
  52     latex.DrawLatex(8.5,18.5,"#sigma=#sqrt{p_{0}+p_{1}#frac{1}{IB}+p_{2}#frac{U1}{U4xIB}+p_{3}#frac{U2}{U4xIB} + p_{4}#frac{U3}{U4xIB}}");
  53     latex.DrawLatex(8.5,17,TString::Format("p_{0}=%1.f",param0[0]));
  54     latex.DrawLatex(8.5,16,TString::Format("p_{1}=%1.f",param0[1]));
  55     latex.DrawLatex(8.5,15,TString::Format("p_{2}=%1.f",param0[2]));
  56     latex.DrawLatex(8.5,14,TString::Format("p_{3}=%1.f",param0[3]));
  57     latex.DrawLatex(8.5,13,TString::Format("p_{4}=%1.f",param0[4]));
  58   }
  59
  60   canvas->SaveAs("canvasFEResolutionFit.pdf");
  61   canvas->SaveAs("canvasFEResolutionFit.png");
  62 }
  63
  64
  65 void EmpiricalGEMQFit(){
  66   //
  67   // Empirical fit of the relative Q resolution for the iron source with 4 GEM layers
  68   // Fit:
  69   // Assumption :
  70   // 1.) Effective electron transparency proportiaonal to the effective ion transparency
  71   // 2.) RMS of the effecitive gain can be expressed as linear function of U1/U3,U2/U4 and U3/U4
  72   //
  73   // Agreement with data within the expect error estimate -relative agreement
  74   // RMS (sigma_{meas}/sigma_{fit}) ~ 3%
  75   //
  76   TFile *fgem = TFile::Open("ETscansIKF.root");
  77   TTree * tree= (TTree*)fgem->Get("NeCO2N2");
  78   tree->SetMarkerStyle(25);
  79   TStatToolkit toolkit;
  80   Double_t chi20=0;
  81   Int_t    npoints=0;
  82   Int_t npointsMax=10000;
  83   TVectorD param0,param1,param2,param3;
  84   TMatrixD covar0,covar1,covar2,covar3;
  85   //
  86   TString fstringFast="";
  87   fstringFast+="1/IB++";                // fraction of the
  88   fstringFast+="(ET1/1000)/IB++";    // fraction of the gain
  89   fstringFast+="(ET2/1000)/IB++";    // fraction of the gain
  90   fstringFast+="(ET3/1000)/IB++";    // fraction of the gain
  91   fstringFast+="(ET4/1000)/IB++";    // fraction of the gain
  92   //fstringFast+="(UGEM1/UGEM4)/IB++";    // fraction of the gain
  93   //
  94   TCut cutFit="ET1<5500";
  95   TString *strResolFit = TStatToolkit::FitPlane(tree,"Sigma^2:Sigma^2", fstringFast.Data(),cutFit, chi20,npoints,param0,covar0,-1,0, npointsMax, 0);
  96   strResolFit->Tokenize("++")->Print();
  97   tree->SetAlias("fitSigma2",strResolFit->Data());
  98   //
  99   gStyle->SetOptTitle(0);
 100   //
 101   TCanvas *canvas = new TCanvas("canvasEmp","canvasEmp",600,500);
 102   canvas->SetBottomMargin(0.15);
 103   canvas->SetRightMargin(0.1);
 104   canvas->SetTicks(1,1);
 105   tree->SetMarkerSize(0.7);
 106   {
 107     tree->Draw("Sigma:sqrt(fitSigma2):sqrt(1/IB)",cutFit,"colz");
 108     TH2F *htemp = (TH2F*)gPad->GetPrimitive("htemp");
 109     htemp->GetXaxis()->SetTitle("#sigma_{fit}(%)");
 110     htemp->GetYaxis()->SetTitle("#sigma_{meas}(%)");
 111     htemp->GetZaxis()->SetTitle("#sqrt{IBF}");
 112     htemp->SetTitle("Fe resolution");
 113     htemp->GetXaxis()->SetLimits(8,20);
 114     htemp->GetYaxis()->SetLimits(8,20);
 115     htemp->Draw("colz");
 116     TLatex latex;
 117     latex.DrawLatex(8.5,18.5,"#sigma=#sqrt{p_{0}+p_{1}#frac{1}{IB}+p_{2}#frac{ET1}{1000xIB}+p_{3}#frac{ET2}{1000xIB} + p_{4}#frac{ET3}{1000xIB}} +  p_{4}#frac{ET4}{1000xIB}} ");
 118     latex.DrawLatex(8.5,17,TString::Format("p_{0}=%1.f",param0[0]));
 119     latex.DrawLatex(8.5,16,TString::Format("p_{1}=%1.f",param0[1]));
 120     latex.DrawLatex(8.5,15,TString::Format("p_{2}=%1.f",param0[2]));
 121     latex.DrawLatex(8.5,14,TString::Format("p_{3}=%1.f",param0[3]));
 122     latex.DrawLatex(8.5,13,TString::Format("p_{4}=%1.f",param0[4]));
 123     latex.DrawLatex(8.5,12,TString::Format("p_{5}=%1.f",param0[5]));
 124   }
 125
 126   canvas->SaveAs("canvasFEResolutionFitET.pdf");
 127   canvas->SaveAs("canvasFEResolutionFitET.png");
 128 }
 129
 130
 131 void FitMCParam(){
 132   //
 133   // Fit the parmaterizations
 134   //
 135   TChain *chain  = AliXRDPROOFtoolkit::MakeChain("outscan_all.list","d",0,100000);
 136
 137
 138 }