ITS/macrosSDD/MakeSDDADCCalib.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TH1F.h>
   3 #include <TH2F.h>
   4 #include <TF1.h>
   5 #include <TPad.h>
   6 #include <TGraphErrors.h>
   7 #include <TROOT.h>
   8 #include <TFile.h>
   9 #include <TTree.h>
  10 #include <TGrid.h>
  11 #include <TGridResult.h>
  12 #include <TMath.h>
  13 #include <TCanvas.h>
  14 #include <TStyle.h>
  15 #include <TLatex.h>
  16 #endif
  17
  18 // For each module, finds the dEdx distribution in 8 drift time intervals. Fits each distribution with LanGaus. For each module builds the histogram of
  19 // Most Probable Value vs drift time and fits it with pol1. Stores the fit parameters vs module number in histograms hmpvModpar0 and hmpvModpar1.
  20
  21 Double_t LangausFun(Double_t *x, Double_t *par);
  22
  23
  24 void MakeSDDADCCalib(TString foldname = ".",TString filename ="QAresults_barrel.root")
  25 {
  26
  27   TFile *fin = new TFile(Form("%s/%s",foldname.Data(),filename.Data()));
  28   TString lname=Form("clistITSAlignQA");
  29   TDirectory *dirFile=(TDirectory*)fin->Get("ITSAlignQA");
  30   TList *cOutput=0x0;
  31   cOutput = (TList*)dirFile->Get(lname.Data());
  32
  33   const Int_t nDrTimeBin=8;
  34   Float_t drTimeLim[nDrTimeBin+1]={0,800,1600,2400,3200,4000,4800,5600,6400};
  35
  36   if(!cOutput){
  37     Printf("E: Cannot open TList %s",lname.Data());
  38     return;
  39   }
  40
  41   TH2F *hdEdxvsDrTime[260];
  42   TF1 *lfun[260];
  43   TH1F *hmpv = new TH1F("hmpv","hmpv",nDrTimeBin,drTimeLim);
  44   TH1F *hsig = new TH1F("hsig","hsig",nDrTimeBin,drTimeLim);
  45   TH1F *hsigl = new TH1F("hsigl","hsigl",nDrTimeBin,drTimeLim);
  46   TH1F *hmpvModpar0 = new TH1F("hmpvModpar0","hmpvModpar0",280,229.5,509.5);
  47   TH1F *hmpvModpar1 = new TH1F("hmpvModpar1","hmpvModpar1",280,229.5,509.5);
  48   TH1F *hsigModpar0 = new TH1F("hsigModpar0","hsigModpar0",280,229.5,509.5);
  49   TH1F *hsigModpar1 = new TH1F("hsigModpar1","hsigModpar1",280,229.5,509.5);
  50   TH1F *hsiglModpar0 = new TH1F("hsiglModpar0","hsiglModpar0",280,229.5,509.5);
  51   TH1F *hsiglModpar1 = new TH1F("hsiglModpar1","hsiglModpar1",280,229.5,509.5);
  52
  53   TString fout="SDDADCCalibResults.root";
  54   TFile *out=new TFile(Form("%s/%s",foldname.Data(),fout.Data()),"recreate");
  55
  56   TCanvas *chdEdxproj=new TCanvas("chdEdxproj","chdEdxproj",1000,800);
  57   chdEdxproj->Divide(4,2);
  58   TCanvas* cmod=new TCanvas("cmod","module",600,900);
  59   cmod->Divide(1,3);
  60   for(Int_t ihist=0;ihist<260;ihist++){//loop on modules
  61   // for(Int_t ihist=0;ihist<10;ihist++){//only 10 to test
  62     Int_t imod=ihist+240;
  63     hdEdxvsDrTime[ihist]=(TH2F*)cOutput->FindObject(Form("hSDDdEdxvsDrTime%i",imod));
  64     chdEdxproj->Clear("D");
  65     //    chdEdxproj->Divide(4,2);
  66     printf("Mod. # %i \n",imod);
  67     hmpv->SetTitle(Form("MPV Mod. # %i",imod));
  68     hmpv->SetName(Form("MPVModule%i",imod));
  69     hsig->SetTitle(Form("Gauss sigma Mod. # %i",imod));
  70     hsigl->SetTitle(Form("Landau sigma Mod. # %i",imod));
  71     hmpv->GetXaxis()->SetTitle("Drift time (ns)");
  72     hsig->GetXaxis()->SetTitle("Drift time (ns)");
  73     hsigl->GetXaxis()->SetTitle("Drift time (ns)");
  74     hmpv->GetYaxis()->SetTitle("MPV (keV)");
  75     hsig->GetYaxis()->SetTitle("Gaussian #sigma (keV)");
  76     hsigl->GetYaxis()->SetTitle("Landau #sigma (keV)");
  77     TH1F *hdEdxproj[nDrTimeBin];
  78     for(Int_t idEdx=0;idEdx<nDrTimeBin;idEdx++){//loop on DrTime Bin
  79       hdEdxproj[idEdx]=(TH1F*)hdEdxvsDrTime[ihist]->ProjectionY(Form("%i",idEdx),hdEdxvsDrTime[ihist]->GetXaxis()->FindBin(drTimeLim[idEdx]),hdEdxvsDrTime[ihist]->GetXaxis()->FindBin(drTimeLim[idEdx+1]));
  80       hdEdxproj[idEdx]->SetTitle(hdEdxvsDrTime[ihist]->GetTitle());
  81       lfun[ihist] = new TF1(Form("LangausFun%d",imod),LangausFun,50.,300.,4); //Langaus fit on a DrTime slice
  82       lfun[ihist]->SetLineWidth(2);
  83       lfun[ihist]->SetParameter(0,5.);
  84       lfun[ihist]->SetParameter(1,80.);
  85       lfun[ihist]->SetParameter(2,hdEdxproj[idEdx]->GetEntries()/10.);
  86       lfun[ihist]->SetParameter(3,10.);
  87       lfun[ihist]->SetParLimits(3,0.,20);
  88       hdEdxproj[idEdx]->Fit(lfun[ihist],"NQLR");
  89       hdEdxproj[idEdx]->GetXaxis()->SetTitle(Form("dE/dx, time interval %d",idEdx));
  90       hdEdxproj[idEdx]->GetYaxis()->SetTitle("Events");
  91       Float_t mpv=lfun[ihist]->GetParameter(1);
  92       Float_t empv=lfun[ihist]->GetParError(1);
  93       Float_t sig=lfun[ihist]->GetParameter(3);
  94       Float_t esig=lfun[ihist]->GetParError(3);
  95       Float_t sigl=lfun[ihist]->GetParameter(0);
  96       Float_t esigl=lfun[ihist]->GetParError(0);
  97       //filling histos mpv vs DrTime for each module
  98       hmpv->Fill(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1]),mpv);
  99       hmpv->SetBinError(hmpv->FindBin(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1])),empv);
 100       hsig->Fill(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1]),sig);
 101       hsig->SetBinError(hsig->FindBin(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1])),esig);
 102       hsigl->Fill(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1]),sigl);
 103       hsigl->SetBinError(hsigl->FindBin(0.5*(drTimeLim[idEdx]+drTimeLim[idEdx+1])),esigl);
 104       chdEdxproj->cd(idEdx+1);
 105       hdEdxproj[idEdx]->Draw();
 106       lfun[ihist]->Draw("same");
 107       chdEdxproj->Update();
 108     }//end loop on DrTime slices
 109     TF1 *pol1mpv = new TF1("pol1mpv","pol1mpv",0,6400);
 110     if(imod!=469)hmpv->Fit(pol1mpv,"NQLR","",1000,5400);
 111     else hmpv->Fit(pol1mpv,"NQLR","",1000,3000); //Mod 469 only one part of the dr Time region is full
 112     if(imod==376){//Hard coded, bad module
 113       pol1mpv->SetParameter(0,84);
 114       pol1mpv->SetParameter(1,0);
 115     }
 116     hmpv->Write();
 117     hmpvModpar1->Fill(imod,pol1mpv->GetParameter(1));
 118     hmpvModpar1->SetBinError(hmpvModpar1->FindBin(imod),pol1mpv->GetParError(1));
 119     Printf("Par0_mpv:%f Err_Par0_mpv:%f    Par1_mpv:%f Err_Par1_mpv:%f",pol1mpv->GetParameter(0),pol1mpv->GetParError(0),pol1mpv->GetParameter(1),pol1mpv->GetParError(1));
 120     hmpvModpar0->Fill(imod,pol1mpv->GetParameter(0));
 121     hmpvModpar0->SetBinError(hmpvModpar0->FindBin(imod),pol1mpv->GetParError(0));
 122     cmod->cd(1);
 123     hmpv->Draw();
 124     pol1mpv->Draw("SAME");
 125     cmod->Update();
 126     hmpv->Reset("M");
 127
 128     TF1 *pol1sig = new TF1("pol1sig","pol1sig",0,6400);
 129     hsig->Fit(pol1sig,"NQLR");
 130     hsigModpar0->Fill(imod,pol1sig->GetParameter(0));
 131     hsigModpar0->SetBinError(hsigModpar0->FindBin(imod),pol1sig->GetParError(0));
 132     hsigModpar1->Fill(imod,pol1sig->GetParameter(1));
 133     hsigModpar1->SetBinError(hsigModpar1->FindBin(imod),pol1sig->GetParError(1));
 134     cmod->cd(2);
 135     hsig->Draw();
 136     pol1sig->Draw("SAME");
 137     cmod->Update();
 138     hsig->Reset("M");
 139
 140     TF1 *pol1sigl = new TF1("pol1sigl","pol1sigl",0,6400);
 141     hsigl->Fit(pol1sigl,"NQLR");
 142     hsiglModpar0->Fill(imod,pol1sigl->GetParameter(0));
 143     hsiglModpar0->SetBinError(hsiglModpar0->FindBin(imod),pol1sigl->GetParError(0));
 144     hsiglModpar1->Fill(imod,pol1sigl->GetParameter(1));
 145     hsiglModpar1->SetBinError(hsiglModpar1->FindBin(imod),pol1sigl->GetParError(1));
 146     cmod->cd(3);
 147     hsigl->Draw();
 148     pol1sigl->Draw("SAME");
 149     cmod->Update();
 150     hsigl->Reset("M");
 151
 152   }//end loop on modules
 153
 154   hmpvModpar0->GetXaxis()->SetTitle("SDD Module");
 155   hmpvModpar1->GetXaxis()->SetTitle("SDD Module");
 156   hsigModpar0->GetXaxis()->SetTitle("SDD Module");
 157   hsigModpar1->GetXaxis()->SetTitle("SDD Module");
 158   hsiglModpar0->GetXaxis()->SetTitle("SDD Module");
 159   hsiglModpar1->GetXaxis()->SetTitle("SDD Module");
 160   hmpvModpar0->GetYaxis()->SetTitle("par0 MPV");
 161   hmpvModpar1->GetYaxis()->SetTitle("par1 MPV");
 162   hsigModpar0->GetYaxis()->SetTitle("par0 sigma Gaus");
 163   hsigModpar1->GetYaxis()->SetTitle("par1 sigma Gaus");
 164   hsiglModpar0->GetYaxis()->SetTitle("par0 sigma Landau");
 165   hsiglModpar1->GetYaxis()->SetTitle("par1 sigma Landau");
 166   Double_t offset=1.3;
 167   hmpvModpar0->GetYaxis()->SetTitleOffset(offset);
 168   hmpvModpar1->GetYaxis()->SetTitleOffset(offset);
 169   hsigModpar0->GetYaxis()->SetTitleOffset(offset);
 170   hsigModpar1->GetYaxis()->SetTitleOffset(offset);
 171   hsiglModpar0->GetYaxis()->SetTitleOffset(offset);
 172   hsiglModpar1->GetYaxis()->SetTitleOffset(offset);
 173
 174   TCanvas *chmpvModpar0=new TCanvas("chmpvModpar0","chmpvModpar0",1000,800);
 175   chmpvModpar0->cd();
 176   hmpvModpar0->Draw();
 177   TCanvas *chmpvModpar1=new TCanvas("chmpvModpar1","chmpvModpar1",1000,800);
 178   chmpvModpar1->cd();
 179   hmpvModpar1->Draw();
 180   TCanvas *chsigModpar0=new TCanvas("chsigModpar0","chsigModpar0",1000,800);
 181   chsigModpar0->cd();
 182   hsigModpar0->Draw();
 183   TCanvas *chsigModpar1=new TCanvas("chsigModpar1","chsigModpar1",1000,800);
 184   chsigModpar1->cd();
 185   hsigModpar1->Draw();
 186   TCanvas *chsiglModpar0=new TCanvas("chsiglModpar0","chsiglModpar0",1000,800);
 187   chsiglModpar0->cd();
 188   hsiglModpar0->SetMinimum(0);
 189   hsiglModpar0->SetMaximum(20);
 190   hsiglModpar0->Draw();
 191   TCanvas *chsiglModpar1=new TCanvas("chsiglModpar1","chsiglModpar1",1000,800);
 192   chsiglModpar1->cd();
 193   hsiglModpar1->Draw();
 194
 195   hmpvModpar0->Write();
 196   hmpvModpar1->Write();
 197   hsigModpar0->Write();
 198   hsigModpar1->Write();
 199   hsiglModpar0->Write();
 200   hsiglModpar1->Write();
 201   out->Close();
 202   delete out;
 203
 204 } //end main
 205
 206
 207 Double_t LangausFun(Double_t *x, Double_t *par) {
 208
 209   //Fit parameters:
 210   //par[0]=Width (scale) parameter of Landau density
 211   //par[1]=Most Probable (MP, location) parameter of Landau density
 212   //par[2]=Total area (integral -inf to inf, normalization constant)
 213   //par[3]=Width (sigma) of convoluted Gaussian function
 214   //
 215   //In the Landau distribution (represented by the CERNLIB approximation),
 216   //the maximum is located at x=-0.22278298 with the location parameter=0.
 217   //This shift is corrected within this function, so that the actual
 218   //maximum is identical to the MP parameter.
 219
 220   // Numeric constants
 221   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 222   Double_t mpshift  = -0.22278298;       // Landau maximum location
 223
 224   // Control constants
 225   Double_t np = 100.0;      // number of convolution steps
 226   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 227
 228   // Variables
 229   Double_t xx;
 230   Double_t mpc;
 231   Double_t fland;
 232   Double_t sum = 0.0;
 233   Double_t xlow,xupp;
 234   Double_t step;
 235   Double_t i;
 236
 237
 238   // MP shift correction
 239   mpc = par[1] - mpshift * par[0];
 240
 241   // Range of convolution integral
 242   xlow = x[0] - sc * par[3];
 243   xupp = x[0] + sc * par[3];
 244
 245   step = (xupp-xlow) / np;
 246
 247   // Convolution integral of Landau and Gaussian by sum
 248   for(i=1.0; i<=np/2; i++) {
 249     xx = xlow + (i-.5) * step;
 250     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 251     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 252
 253     xx = xupp - (i-.5) * step;
 254     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 255     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 256   }
 257
 258   return (par[2] * step * sum * invsq2pi / par[3]);
 259 }