ITS/GetGainModuleLevel.C

   1 #include "TMath.h"
   2 #include "TList.h"
   3 #include "TH2F.h"
   4 #include "TH1F.h"
   5 #include "TF1.h"
   6 #include "TH1D.h"
   7 #include "TProfile.h"
   8 #include "TFile.h"
   9 #include "TCanvas.h"
  10 #include <Riostream.h>
  11 #include <TROOT.h>
  12 #include <TString.h>
  13 #include "TStyle.h"
  14 #include "TGrid.h"
  15
  16
  17 class TList;
  18 class TH2F;
  19 class TH1F;
  20 class TH1D;
  21 class TProfile;
  22 class TFile;
  23 class TCanvas;
  24 class TStyle;
  25 class TF1;
  26 class TGrid;
  27
  28
  29
  30
  31
  32 Double_t convolution(Double_t* x,Double_t *par)
  33 {
  34           //Fit parameters:
  35    //par[0]=Width (scale) parameter of Landau density
  36    //par[1]=Most Probable (MP, location) parameter of Landau density
  37    //par[2]=Total area (integral -inf to inf, normalization constant)
  38    //par[3]=Width (sigma) of convoluted Gaussian function
  39    //
  40    //In the Landau distribution (represented by the CERNLIB approximation),
  41    //the maximum is located at x=-0.22278298 with the location parameter=0.
  42    //This shift is corrected within this function, so that the actual
  43    //maximum is identical to the MP parameter.
  44
  45       // Numeric constants
  46       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  47       Double_t mpshift  = -0.22278298;       // Landau maximum location
  48
  49       // Control constants
  50       Double_t np = 200.0;      // number of convolution steps
  51       Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
  52
  53       // Variables
  54       Double_t xx;
  55       Double_t mpc;
  56       Double_t fland;
  57       Double_t sum = 0.0;
  58       Double_t xlow,xupp;
  59       Double_t step;
  60       Double_t i;
  61
  62
  63       // MP shift correction
  64       mpc = par[1] - mpshift * par[0];
  65
  66       // Range of convolution integral
  67       xlow = x[0] - sc * par[3];
  68       xupp = x[0] + sc * par[3];
  69
  70       step = (xupp-xlow) / np;
  71
  72       // Convolution integral of Landau and Gaussian by sum
  73       for(i=1.0; i<=np/2; i++)
  74       {
  75                  xx = xlow + (i-.5) * step;
  76                 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
  77                 sum += fland * TMath::Gaus(x[0],xx,par[3]);
  78
  79                 xx = xupp - (i-.5) * step;
  80                 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
  81                 sum += fland * TMath::Gaus(x[0],xx,par[3]);
  82       }
  83
  84       return (par[2] * step * sum * invsq2pi / par[3]);
  85 }
  86
  87
  88
  89
  90 void GetGainModuleLevel(TString filename,Bool_t normal=1,Int_t ntofit=500, Bool_t grid=0)
  91 {
  92         gROOT->SetStyle("Plain");
  93         gStyle->SetPalette(1,0);
  94         gStyle->SetOptStat(111111);
  95         gStyle->SetOptFit(1);
  96         if(grid)
  97                 TGrid::Connect("alien://");
  98
  99         TFile* file_data=TFile::Open(filename);
 100         if(!file_data)
 101                 return;
 102         TList *listin=0x0;
 103         listin=(TList*)file_data->Get("output");
 104         if(!listin)
 105                 listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/output");
 106         if(!listin)
 107                 listin=(TList*)file_data->Get("PWGPPdEdxSSDQA/SSDdEdxQA");
 108         if(!listin)
 109                 listin=(TList*)file_data->Get("SSDdEdxQA");
 110         if(!listin)
 111                 return;
 112         TH2F* fHistQ=0x0;
 113         if (normal)
 114                 fHistQ=(TH2F*)listin ->FindObject("QACharge");
 115         else
 116                 fHistQ=(TH2F*)listin ->FindObject("QAChargeCorrected");
 117         if(!fHistQ)
 118                 return;
 119         TH2F* fHistCR=(TH2F*)listin ->FindObject("QAChargeRatio");
 120         if(!fHistCR)
 121                 return;
 122         TList *listout1=new TList();
 123
 124         TList *listout2=new TList();
 125
 126         TH1F* fHistMPVs=new TH1F("HistMPVS","HistMPVs;MPV;N",75,70,95);
 127         fHistMPVs->SetDirectory(0);
 128         listout2->Add(fHistMPVs);
 129
 130         TH1F* fHistSL=new TH1F("HistSL","#sigma_{Landau};#sigma_{Landau};N",40,0,16);
 131         fHistSL->SetDirectory(0);
 132         listout2->Add(fHistSL);
 133
 134         TH1F* fHistSG=new TH1F("HistSG","#sigma_{Gaus};#sigma_{Gaus};N",40,0,16);
 135         fHistSG->SetDirectory(0);
 136         listout2->Add(fHistSG);
 137
 138         TH1F* fHistCRmean=new TH1F("HistCRmean","HistCRmean;CRmean;N",200,-1,1);
 139         fHistCRmean->SetDirectory(0);
 140         listout2->Add(fHistCRmean);
 141
 142         TH1F* fHistCRRMS=new TH1F("HistCRRMS","HistCRRMS;CRRMS;N",100,0,1);
 143         fHistCRRMS->SetDirectory(0);
 144         listout2->Add(fHistCRRMS);
 145
 146         TH1F* fHistGainP=new TH1F("HistGainP","HistGainP;CRGainPcorr;N",120,0.5,2.0);
 147         fHistGainP->SetDirectory(0);
 148         listout2->Add(fHistGainP);
 149
 150         TH1F* fHistGainN=new TH1F("HistGainN","HistGainN;CRGainNcorr;N",120,0.5,2.0);
 151         fHistGainN->SetDirectory(0);
 152         listout2->Add(fHistGainN);
 153
 154         TH1F *fMPVGraph = new TH1F("MPVgraph","MPVgraph;Module number;MPV",1698,-0.5,1697.5);
 155         fMPVGraph->SetMarkerColor(kRed);
 156         fMPVGraph->SetMarkerStyle(22);
 157         listout2->Add(fMPVGraph);
 158
 159         TH1F *fCRmeanGraph = new TH1F("CRmeangraph","CRmeangraph;Module number;MPV",1698,-0.5,1697.5);
 160         fCRmeanGraph->SetMarkerColor(kBlue);
 161         fCRmeanGraph->SetMarkerStyle(23);
 162         listout2->Add(fCRmeanGraph);
 163
 164         Float_t gainP[1698];
 165         Float_t gainN[1698];
 166         Float_t mpv[1698];
 167         Int_t flag[1698];
 168
 169         ofstream outfiletxt;
 170         outfiletxt.open("gain.txt");
 171          outfiletxt.width(10) ;
 172         outfiletxt.setf(outfiletxt.left);
 173         outfiletxt<<"MODULE"<<"\t";
 174         outfiletxt.width(10);
 175         outfiletxt.setf(outfiletxt.left);
 176         outfiletxt<<"FLAG"<<"\t";
 177         outfiletxt.width(10) ;
 178         outfiletxt.setf(outfiletxt.left);
 179         outfiletxt<<"GainPcorr"<<"\t";
 180         outfiletxt.width(10) ;
 181         outfiletxt.setf(outfiletxt.left);
 182         outfiletxt<<"GainNcorr"<<"\t";
 183         outfiletxt.width(10) ;
 184         outfiletxt.setf(outfiletxt.left);
 185         outfiletxt<<"MPV"<<endl;
 186
 187
 188         ofstream outfiletxtbad;
 189         outfiletxtbad.open("badModules.txt");
 190
 191
 192
 193
 194         for (int i =0;i<1698;i++)
 195         {
 196                 cout<<i<<endl;
 197                 TString tmpQ("Q");
 198                 tmpQ+=i;
 199                 TString tmpCR("CR");
 200                 tmpCR+=i;
 201                 TH1D* fHist1DCR= fHistCR->ProjectionY(tmpCR,i+1,i+1);
 202                 Double_t mean=fHist1DCR->GetMean();
 203                 if(!(TMath::Abs(mean)<1.0)||fHist1DCR->GetEntries()<10)
 204                 {
 205                         flag[i]=-2;
 206                         gainN[i]=1.0;
 207                         gainP[i]=1.0;
 208                         mpv[i]=1.0;
 209                         continue;
 210                 }
 211                 fHistCRmean->Fill(mean);
 212                 fHistCRRMS->Fill(fHist1DCR->GetRMS());
 213                 gainN[i]=1.0/(1.0+mean);
 214                 gainP[i]=1.0/(1.0-mean);
 215                 fHistGainP->Fill(gainP[i]);
 216                 fHistGainN->Fill(gainN[i]);
 217                 fCRmeanGraph->SetBinContent(i+1,mean);
 218                 fCRmeanGraph->SetBinError(i+1,fHist1DCR->GetRMS());
 219
 220
 221                 TH1D* fHist1DQ=fHistQ->ProjectionY(tmpQ,i+1,i+1);
 222                  fHist1DQ->SetDirectory(0);
 223                  listout1->Add(fHist1DQ);
 224                 if(fHist1DQ->GetEntries()<ntofit)
 225                 {
 226                         flag[i]=-1;
 227                         mpv[i]=1.0;
 228                          outfiletxtbad<<"Low statistic \t module= "<<i<<" netries="<<fHist1DQ->GetEntries()<<endl;
 229                         continue;
 230                 }
 231                 else
 232                 {
 233                         tmpQ+="fit";
 234                         Float_t range=fHist1DQ->GetBinCenter(fHist1DQ->GetMaximumBin());
 235                         TF1 *f1 = new TF1(tmpQ,convolution,range*0.45,range*3.0,4);
 236                         f1->SetParameters(7.0,range,1.0,5.5);
 237                         Float_t normalization=fHist1DQ->GetEntries()*fHist1DQ->GetXaxis()->GetBinWidth(2)/f1->Integral(range*0.45,range*3.0);
 238                         f1->SetParameters(7.0,range,normalization,5.5);
 239                         //f1->SetParameters(7.0,range,fHist1DQ->GetMaximum(),5.5);
 240                         f1->SetParNames("sigma Landau","MPV","N","sigma Gaus");
 241                         f1->SetParLimits(0,2.0,100.0);
 242                         f1->SetParLimits(3,0.0,100.0);
 243                         if(fHist1DQ->Fit(tmpQ,"BRQ")==0)
 244                         {
 245                                 mpv[i]=f1->GetParameter(1);
 246                                 fHistMPVs->Fill(mpv[i]);
 247                                 fHistSL->Fill(f1->GetParameter(0));
 248                                 fHistSG->Fill(f1->GetParameter(3));
 249                                 flag[i]=1;
 250                                 fMPVGraph->SetBinContent(i+1,f1->GetParameter(1));
 251                                 fMPVGraph->SetBinError(i+1,f1->GetParError(1));
 252                                 if(mpv[i]<75.0)
 253                                 {
 254                                         outfiletxtbad<<"MPV lower than 75 \t module="<<i<<endl;
 255                                         flag[i]=0;
 256                                 }
 257                                 if(mpv[i]>100.0)
 258                                 {
 259                                         outfiletxtbad<<"MPV higher than 100 \t module="<<i<<endl;
 260                                         flag[i]=0;
 261
 262                                 }
 263                                 if(f1->GetParError(1)>1.0)
 264                                 {
 265                                         outfiletxtbad<<"MPV high error on MPV  \t module="<<i<<endl;
 266                                         flag[i]=0;
 267                                 }
 268                         }
 269                         else
 270                         {
 271                                 mpv[i]=1;
 272                                 flag[i]=0;
 273                                  outfiletxtbad<<"BAD FIT \t module="<<i<<endl;
 274                                 //continue;
 275                         }
 276                 }
 277         }
 278
 279         for (int i=0;i<1698;i++)
 280         {
 281                 outfiletxt.setf(outfiletxt.scientific);
 282                 outfiletxt.precision(2);
 283                  outfiletxt.width(10) ;
 284                 outfiletxt.setf(outfiletxt.left);
 285                 outfiletxt<<i<<"\t";
 286                  outfiletxt.width(10) ;
 287                 outfiletxt.setf(outfiletxt.left);
 288                 outfiletxt<<flag[i]<<"\t";
 289                  outfiletxt.width(10) ;
 290                 outfiletxt.setf(outfiletxt.left);
 291                 outfiletxt<<gainP[i]<<"\t";
 292                  outfiletxt.width(10) ;
 293                 outfiletxt.setf(outfiletxt.left);
 294                 outfiletxt<<gainN[i]<<"\t";
 295                  outfiletxt.width(10) ;
 296                 outfiletxt.setf(outfiletxt.left);
 297                 outfiletxt<<mpv[i]<<endl;
 298         }
 299
 300         TCanvas *c1 = new TCanvas("1","1",1200,800);
 301         c1->Divide(2,1);
 302         c1->cd(1);
 303         fHistQ->Draw("colz");
 304         c1->cd(2);
 305         fHistCR->Draw("colz");
 306
 307
 308
 309         TFile* fout1=TFile::Open("gain_all_fits.root","recreate");
 310         listout1->Write("output",TObject::kSingleKey);
 311         fout1->Close();
 312
 313         TFile* fout2=TFile::Open("gain_control_plots.root","recreate");
 314         listout2->Write("output",TObject::kSingleKey);
 315         fout2->Close();
 316
 317
 318
 319 }