ITS/macrosSDD/PlotOutputQAtrainSDD.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 #include "AliITSgeomTGeo.h"
  17 #endif
  18
  19 Double_t LangausFun(Double_t *x, Double_t *par);
  20
  21
  22 void PlotOutputQAtrainSDD(TString option="local",
  23                           Int_t nRun=0,
  24                           TString period="LHC11a",
  25                           TString qaTrain="",
  26                           TString fileName="QAresults.root"){
  27
  28   gStyle->SetOptStat(0);
  29   //  gStyle->SetOptTitle(0);
  30   TFile *f;
  31   TString path;
  32   Int_t year=2011;
  33   if(period.Contains("LHC10")) year=2010;
  34   else if(period.Contains("LHC09")) year=2009;
  35
  36   if(option.Contains("local")){
  37     f=new TFile(fileName.Data());
  38     printf("Opened file %s\n",f->GetName());
  39   }else{
  40     TGrid::Connect("alien:");
  41     path=Form("/alice/data/%d/%s/%09d/ESDs/",year,period.Data(),nRun);
  42     printf("search in path %s\n",path.Data());
  43     if(!gGrid||!gGrid->IsConnected()) {
  44       printf("gGrid not found! exit macro\n");
  45       return;
  46     }
  47     TGridResult *gr = gGrid->Query(path,fileName);
  48     Int_t nFiles = gr->GetEntries();
  49     if (nFiles < 1){
  50       printf("QA file for run %d not found\n",nRun);
  51       return;
  52     }
  53     printf("================>%d files found\n", nFiles);
  54     if(qaTrain.Length()>0){
  55       Int_t found=kFALSE;
  56       for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
  57         fileName=Form("%s",gr->GetKey(iFil,"turl"));
  58         TString isMerged=fileName;
  59         isMerged.Remove(isMerged.Sizeof()-16);
  60         isMerged.Remove(0,isMerged.Sizeof()-5);
  61         if(!isMerged.Contains("QA")) continue;
  62         if(fileName.Contains(qaTrain.Data())){
  63           found=kTRUE;
  64           break;
  65         }
  66       }
  67       if(!found){
  68         printf(" File from %s train not found\n",qaTrain.Data());
  69         return;
  70       }
  71     }else{
  72       Int_t theFile=0;
  73       Int_t maxVer=0;
  74       if (nFiles > 1){
  75         for (Int_t iFil = 0; iFil <nFiles ; iFil++) {
  76           fileName=Form("%s",gr->GetKey(iFil,"turl"));
  77           TString isMerged=fileName;
  78           isMerged.Remove(isMerged.Sizeof()-16);
  79           isMerged.Remove(0,isMerged.Sizeof()-5);
  80           if(!isMerged.Contains("QA")) continue;
  81           TString cutFilename=fileName.ReplaceAll("/QAresults.root","");
  82           Int_t last=cutFilename.Sizeof()-3;
  83           cutFilename.Remove(0,last);
  84           Int_t qaver=atoi(cutFilename.Data());
  85           if(qaver>maxVer){
  86             maxVer=qaver;
  87             theFile=iFil;
  88           }
  89         }
  90       }
  91       fileName=Form("%s",gr->GetKey(theFile,"turl"));
  92     }
  93     f=TFile::Open(fileName.Data());
  94   }
  95
  96   TDirectoryFile* df=(TDirectoryFile*)f->Get("SDD_Performance");
  97   if(!df){
  98     printf("SDD_Performance MISSING -> Exit\n");
  99     return;
 100   }
 101   TList* l=(TList*)df->Get("coutputRP");
 102   if(!df){
 103     printf("coutputRP TList MISSING -> Exit\n");
 104     return;
 105   }
 106
 107   TH1F* hcllay=(TH1F*)l->FindObject("hCluInLay");
 108   TH1F* hapmod=(TH1F*)l->FindObject("hAllPmod");
 109   TH1F* hgpmod=(TH1F*)l->FindObject("hGoodPmod");
 110   TH1F* hmpmod=(TH1F*)l->FindObject("hMissPmod");
 111   TH1F* hbrmod=(TH1F*)l->FindObject("hBadRegmod");
 112   TH1F* hskmod=(TH1F*)l->FindObject("hSkippedmod");
 113   TH1F* hoamod=(TH1F*)l->FindObject("hOutAccmod");
 114   TH1F* hnrmod=(TH1F*)l->FindObject("hNoRefitmod");
 115
 116   //  TH1F* hapxl=(TH1F*)l->FindObject("hAllPxloc");
 117   TH1F* hgpxl=(TH1F*)l->FindObject("hGoodPxloc");
 118   TH1F* hmpxl=(TH1F*)l->FindObject("hMissPxloc");
 119   TH1F* hbrxl=(TH1F*)l->FindObject("hBadRegxloc");
 120   //  TH1F* hapzl=(TH1F*)l->FindObject("hAllPzloc");
 121   TH1F* hgpzl=(TH1F*)l->FindObject("hGoodPzloc");
 122   TH1F* hmpzl=(TH1F*)l->FindObject("hMissPzloc");
 123   TH1F* hbrzl=(TH1F*)l->FindObject("hBadRegzloc");
 124
 125   TH2F* hClSizAn=(TH2F*)l->FindObject("hCluSizAn");
 126   TH2F* hClSizTb=(TH2F*)l->FindObject("hCluSizTb");
 127
 128   TH2F* hdedx3=(TH2F*)l->FindObject("hdEdxL3VsP");
 129   TH2F* hdedx4=(TH2F*)l->FindObject("hdEdxL4VsP");
 130   TH2F* hdedxmod=(TH2F*)l->FindObject("hdEdxVsMod");
 131
 132
 133   TH1F* hmodR=(TH1F*)l->FindObject("hRPMod");
 134   TH1F* hmodT=(TH1F*)l->FindObject("hTPMod");
 135   TH1F* hgamod=(TH1F*)l->FindObject("hGAMod");
 136
 137   TH2F* h2dmodR3=new TH2F("h2dmodR3","Rec Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
 138   TH2F* h2dmodR4=new TH2F("h2dmodR4","Rec Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
 139   TH2F* h2dmodT3=new TH2F("h2dmodT3","Track Points, Layer 3",6,0.5,6.5,14,0.5,14.5);
 140   TH2F* h2dmodT4=new TH2F("h2dmodT4","Track Points, Layer 4",8,0.5,8.5,22,0.5,22.5);
 141   TH2F* h2dmodR3N=new TH2F("h2dmodR3N","Rec Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
 142   TH2F* h2dmodR4N=new TH2F("h2dmodR4N","Rec Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
 143   TH2F* h2dmodT3N=new TH2F("h2dmodT3N","Track Points/GoodAnode/Event, Layer 3",6,0.5,6.5,14,0.5,14.5);
 144   TH2F* h2dmodT4N=new TH2F("h2dmodT4N","Track Points/GoodAnode/Event, Layer 4",8,0.5,8.5,22,0.5,22.5);
 145   TH1F* hmodRN=new TH1F("hmodRN","Normalized Rec Points per Module",260,239.5,499.5);
 146   TH1F* hmodTN=new TH1F("hmodTN","Normalized Track Points per Module",260,239.5,499.5);
 147
 148   TH1F* hev=(TH1F*)l->FindObject("hNEvents");
 149   Int_t nTotEvents=hev->GetBinContent(2);
 150   Int_t nTrigEvents=hev->GetBinContent(3);
 151   Int_t nEvents=nTotEvents;
 152   printf("---- Statistics ----\n");
 153   printf("Number of Events = %d\n",nTotEvents);
 154   if(nTrigEvents>0){
 155     printf("Number of Triggered Events = %d\n",nTrigEvents);
 156     nEvents=nTrigEvents;
 157   }else{
 158     printf("No request on the trigger done when running the task\n");
 159   }
 160   Int_t bestMod=0;
 161   for(Int_t iMod=0; iMod<260;iMod++){
 162     Int_t gda=(Int_t)hgamod->GetBinContent(iMod+1);
 163     if(gda>bestMod) bestMod=gda;
 164   }
 165   Int_t nChunks=1;
 166   if(bestMod>512){
 167     nChunks=(Int_t)(bestMod/512.+0.5);
 168   }
 169   printf("Chunks merged = %d\n",nChunks);
 170   hgamod->Scale(1./nChunks);
 171   TCanvas* cgan=new TCanvas("cgan","Good Anodes");
 172   cgan->SetTickx();
 173   cgan->SetTicky();
 174   hgamod->SetMarkerStyle(20);
 175   hgamod->SetMarkerSize(0.6);
 176   hgamod->SetMinimum(-10.);
 177   hgamod->Draw("P");
 178   hgamod->GetXaxis()->SetTitle("SDD Module Id");
 179   hgamod->GetYaxis()->SetTitle("Number of good anodes");
 180   cgan->Update();
 181
 182   printf("---- Modules with > 2%% of bad anodes ----\n");
 183   for(Int_t iMod=0; iMod<260; iMod++){
 184     Int_t idMod=iMod+240;
 185     Float_t rps=hmodR->GetBinContent(iMod+1);
 186     Float_t tps=hmodT->GetBinContent(iMod+1);
 187     Float_t ga=hgamod->GetBinContent(iMod+1);
 188     if(ga<500){
 189       printf("Module %d - Good Anodes = %d\n",idMod,(Int_t)ga);
 190     }
 191     Float_t rpsN=0.;
 192     Float_t tpsN=0.;
 193     Float_t erpsN=0.;
 194     Float_t etpsN=0.;
 195     if(ga>0 && nEvents>0){
 196       rpsN=rps/ga/(Float_t)nEvents;
 197       tpsN=tps/ga/(Float_t)nEvents;
 198       erpsN=TMath::Sqrt(rps)/ga/(Float_t)nEvents;
 199       etpsN=TMath::Sqrt(tps)/ga/(Float_t)nEvents;
 200     }
 201     hmodRN->SetBinContent(iMod+1,rpsN);
 202     hmodTN->SetBinContent(iMod+1,tpsN);
 203     hmodRN->SetBinError(iMod+1,erpsN);
 204     hmodTN->SetBinError(iMod+1,etpsN);
 205     Int_t iLay,iLad,iDet;
 206     AliITSgeomTGeo::GetModuleId(idMod,iLay,iLad,iDet);
 207     if(iLay==3){
 208       h2dmodR3->SetBinContent(iDet,iLad,rps);
 209       h2dmodT3->SetBinContent(iDet,iLad,tps);
 210       h2dmodR3N->SetBinContent(iDet,iLad,rpsN);
 211       h2dmodT3N->SetBinContent(iDet,iLad,tpsN);
 212     }
 213     else if(iLay==4){
 214       h2dmodR4->SetBinContent(iDet,iLad,rps);
 215       h2dmodT4->SetBinContent(iDet,iLad,tps);
 216       h2dmodR4N->SetBinContent(iDet,iLad,rpsN);
 217       h2dmodT4N->SetBinContent(iDet,iLad,tpsN);
 218     }
 219   }
 220   if(nEvents<1) return;
 221
 222   gStyle->SetPalette(1);
 223
 224   if(hmodR->GetEntries()>0){
 225     TCanvas* cmodR=new TCanvas("cmodR","RecPoint Occup",1200,1200);
 226     cmodR->Divide(2,3);
 227     cmodR->cd(1);
 228     gPad->SetLeftMargin(0.14);
 229     hmodR->Draw();
 230     hmodR->GetXaxis()->SetTitle("SDD Module Id");
 231     hmodR->GetYaxis()->SetTitle("RecPoints");
 232     hmodR->GetYaxis()->SetTitleOffset(1.55);
 233     cmodR->cd(2);
 234     gPad->SetLeftMargin(0.14);
 235     hmodRN->Draw("E");
 236     hmodRN->GetXaxis()->SetTitle("SDD Module Id");
 237     hmodRN->GetYaxis()->SetTitle("RecPoints/GoodAnode/Event");
 238     hmodRN->GetYaxis()->SetTitleOffset(1.55);
 239     cmodR->cd(3);
 240     gPad->SetLeftMargin(0.14);
 241     h2dmodR3->Draw("colz");
 242     h2dmodR3->GetXaxis()->SetTitle("Detector");
 243     h2dmodR3->GetYaxis()->SetTitle("Ladder");
 244     cmodR->cd(4);
 245     gPad->SetLeftMargin(0.14);
 246     h2dmodR3N->Draw("colz");
 247     h2dmodR3N->GetXaxis()->SetTitle("Detector");
 248     h2dmodR3N->GetYaxis()->SetTitle("Ladder");
 249     cmodR->cd(5);
 250     gPad->SetLeftMargin(0.14);
 251     h2dmodR4->Draw("colz");
 252     h2dmodR4->GetXaxis()->SetTitle("Detector");
 253     h2dmodR4->GetYaxis()->SetTitle("Ladder");
 254     cmodR->cd(6);
 255     gPad->SetLeftMargin(0.14);
 256     gPad->SetLeftMargin(0.14);
 257     h2dmodR4N->Draw("colz");
 258     h2dmodR4N->GetXaxis()->SetTitle("Detector");
 259     h2dmodR4N->GetYaxis()->SetTitle("Ladder");
 260     cmodR->Update();
 261   }
 262
 263
 264   if(hmodT->GetEntries()>0){
 265     TCanvas* cmodT=new TCanvas("cmodT","TrackPoint Occup",1200,1200);
 266     cmodT->Divide(2,3);
 267     cmodT->cd(1);
 268     hmodT->Draw();
 269     hmodT->GetXaxis()->SetTitle("SDD Module Id");
 270     hmodT->GetYaxis()->SetTitle("TrackPoints");
 271     hmodT->GetYaxis()->SetTitleOffset(1.4);
 272     cmodT->cd(2);
 273     gPad->SetLeftMargin(0.14);
 274     hmodTN->Draw("E");
 275     hmodTN->GetXaxis()->SetTitle("SDD Module Id");
 276     hmodTN->GetYaxis()->SetTitle("TrackPoints");
 277     hmodTN->GetYaxis()->SetTitleOffset(1.4);
 278     cmodT->cd(3);
 279     gPad->SetLeftMargin(0.14);
 280     h2dmodT3->Draw("colz");
 281     h2dmodT3->GetXaxis()->SetTitle("Detector");
 282     h2dmodT3->GetYaxis()->SetTitle("Ladder");
 283     cmodT->cd(4);
 284     gPad->SetLeftMargin(0.14);
 285     h2dmodT3N->Draw("colz");
 286     h2dmodT3N->GetXaxis()->SetTitle("Detector");
 287     h2dmodT3N->GetYaxis()->SetTitle("Ladder");
 288     cmodT->cd(5);
 289     gPad->SetLeftMargin(0.14);
 290     h2dmodT4->Draw("colz");
 291     h2dmodT4->GetXaxis()->SetTitle("Detector");
 292     h2dmodT4->GetYaxis()->SetTitle("Ladder");
 293     cmodT->cd(6);
 294     gPad->SetLeftMargin(0.14);
 295     h2dmodT4N->Draw("colz");
 296     h2dmodT4N->GetXaxis()->SetTitle("Detector");
 297     h2dmodT4N->GetYaxis()->SetTitle("Ladder");
 298     cmodT->Update();
 299   }
 300
 301   TH1F* htplad3=(TH1F*)l->FindObject("hTPLad3");
 302   TH1F* htplad4=(TH1F*)l->FindObject("hTPLad4");
 303   TH1F* hgalad3=(TH1F*)l->FindObject("hGALad3");
 304   TH1F* hgalad4=(TH1F*)l->FindObject("hGALad4");
 305   TH1F* hnormOcc3=new TH1F("hnormOcc3","",14,-0.5,13.5);
 306   TH1F* hnormOcc4=new TH1F("hnormOcc4","",22,-0.5,21.5);
 307   Bool_t tpok=kFALSE;
 308   for(Int_t ilad=0;ilad<14;ilad++){
 309     Float_t occ=0.;
 310     Float_t eocc=0.;
 311     Int_t gd3=hgalad3->GetBinContent(ilad+1);
 312     if(gd3>0){
 313       occ=(Float_t)htplad3->GetBinContent(ilad+1)/(Float_t)gd3/(Float_t)nEvents;
 314       eocc=TMath::Sqrt((Float_t)htplad3->GetBinContent(ilad+1))/(Float_t)gd3/(Float_t)nEvents;
 315     }
 316     hnormOcc3->SetBinContent(ilad+1,occ);
 317     hnormOcc3->SetBinError(ilad+1,eocc);
 318   }
 319   for(Int_t ilad=0;ilad<22;ilad++){
 320     Float_t occ=0.;
 321     Float_t eocc=0.;
 322     Int_t gd4=hgalad4->GetBinContent(ilad+1);
 323     if(gd4>0){
 324       occ=(Float_t)htplad4->GetBinContent(ilad+1)/(Float_t)gd4/(Float_t)nEvents;
 325       eocc=TMath::Sqrt((Float_t)htplad4->GetBinContent(ilad+1))/(Float_t)gd4/(Float_t)nEvents;
 326     }
 327     hnormOcc4->SetBinContent(ilad+1,occ);
 328     hnormOcc4->SetBinError(ilad+1,eocc);
 329   }
 330
 331
 332   if(tpok){
 333     TCanvas* cn0=new TCanvas("cn0","Normalized Ladder Occupancy",1400,600);
 334     cn0->Divide(2,1);
 335     cn0->cd(1);
 336     gPad->SetLeftMargin(0.14);
 337     hnormOcc3->Draw();
 338     hnormOcc3->GetXaxis()->SetTitle("Ladder number (layer 3)");
 339     hnormOcc3->GetYaxis()->SetTitle("TrackPoints/GoodAnodes/Events");
 340     hnormOcc3->GetYaxis()->SetTitleOffset(1.35);
 341     cn0->cd(2);
 342     gPad->SetLeftMargin(0.14);
 343     hnormOcc4->Draw();
 344     hnormOcc4->GetXaxis()->SetTitle("Ladder number (layer 4)");
 345     hnormOcc4->GetYaxis()->SetTitle("TrackPoints/GoodAnode/Events");
 346     hnormOcc4->GetYaxis()->SetTitleOffset(1.35);
 347     cn0->Update();
 348   }
 349
 350   if(hcllay){
 351     Double_t norm=hcllay->GetBinContent(1);
 352     if(norm>0.){
 353       hcllay->Scale(1./norm);
 354       hcllay->SetTitle("");
 355       hcllay->GetXaxis()->SetRange(2,7);
 356       hcllay->SetMinimum(0.);
 357       hcllay->SetMaximum(1.1);
 358       hcllay->SetMarkerStyle(23);
 359       TCanvas* ceffL=new TCanvas("ceffL","PointPerLayer",1000,600);
 360       ceffL->SetGridy();
 361       hcllay->Draw();
 362       hcllay->GetXaxis()->SetTitle("Layer");
 363       hcllay->GetYaxis()->SetTitle("Fraction of tracks with point in layer");
 364       ceffL->Update();
 365     }
 366   }
 367
 368   hgpmod->SetTitle("");
 369   TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
 370   hgpmod->Draw();
 371   hgpmod->GetXaxis()->SetTitle("SDD Module Id");
 372   hgpmod->GetYaxis()->SetTitle("Number of tracks");
 373   hmpmod->SetLineColor(2);
 374   hmpmod->SetMarkerColor(2);
 375   hmpmod->SetMarkerStyle(22);
 376   hmpmod->SetMarkerSize(0.5);
 377   hmpmod->Draw("psame");
 378   hbrmod->SetLineColor(kGreen+1);
 379   hbrmod->SetMarkerColor(kGreen+1);
 380   hbrmod->SetMarkerStyle(20);
 381   hbrmod->SetMarkerSize(0.5);
 382   hbrmod->Draw("same");
 383   hskmod->SetLineColor(kYellow);
 384   hskmod->Draw("same");
 385   hoamod->SetLineColor(4);
 386   hoamod->Draw("same");
 387   hnrmod->SetLineColor(6);
 388   hnrmod->Draw("same");
 389   TLatex* t1=new TLatex(0.7,0.85,"Good Point");
 390   t1->SetNDC();
 391   t1->SetTextColor(1);
 392   t1->Draw();
 393   TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
 394   t2->SetNDC();
 395   t2->SetTextColor(2);
 396   t2->Draw();
 397   TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
 398   t3->SetNDC();
 399   t3->SetTextColor(kGreen+1);
 400   t3->Draw();
 401   ceff0->Update();
 402
 403   TH1F* heff=new TH1F("heff","",260,239.5,499.5);
 404   for(Int_t imod=0; imod<260;imod++){
 405     Float_t numer=hgpmod->GetBinContent(imod+1)+hbrmod->GetBinContent(imod+1)+hoamod->GetBinContent(imod+1)+hnrmod->GetBinContent(imod+1);
 406     Float_t denom=hapmod->GetBinContent(imod+1);
 407     Float_t eff=0.;
 408     Float_t erreff=0.;
 409     if(denom>0){
 410       eff=numer/denom;
 411       erreff=TMath::Sqrt(eff*(1-eff)/denom);
 412     }
 413     heff->SetBinContent(imod+1,eff);
 414     heff->SetBinError(imod+1,erreff);
 415   }
 416
 417   printf("---- Modules with efficiency < 90%% ----\n");
 418   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,600);
 419   heff->Draw();
 420   heff->GetXaxis()->SetTitle("SDD Module Id");
 421   heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
 422   for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
 423     Float_t e=heff->GetBinContent(ibin);
 424     if(e<0.9){
 425       Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
 426       Int_t lay,lad,det;
 427       AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
 428       printf("Module %d - Layer %d Ladder %2d Det %d  -   Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
 429     }
 430   }
 431   ceff1->Update();
 432
 433   if(hgpxl){
 434     hgpxl->SetTitle("");
 435     hgpzl->SetTitle("");
 436     TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
 437     ceff2->Divide(2,1);
 438     ceff2->cd(1);
 439     hgpxl->Draw();
 440     hgpxl->GetXaxis()->SetTitle("Xlocal (cm)");
 441     hgpxl->GetYaxis()->SetTitle("Number of tracks");
 442     hmpxl->SetLineColor(2);
 443     hmpxl->SetMarkerColor(2);
 444     hmpxl->SetMarkerStyle(22);
 445     hmpxl->SetMarkerSize(0.5);
 446     hmpxl->Draw("psame");
 447     hbrxl->SetLineColor(kGreen+1);
 448     hbrxl->SetMarkerColor(kGreen+1);
 449     hbrxl->SetMarkerStyle(20);
 450     hbrxl->SetMarkerSize(0.5);
 451     hbrxl->Draw("same");
 452     t1->Draw();
 453     t2->Draw();
 454     t3->Draw();
 455     ceff2->cd(2);
 456     hgpzl->Draw();
 457     hgpzl->GetXaxis()->SetTitle("Zlocal (cm)");
 458     hgpzl->GetYaxis()->SetTitle("Number of tracks");
 459     hmpzl->SetLineColor(2);
 460     hmpzl->SetMarkerColor(2);
 461     hmpzl->SetMarkerStyle(22);
 462     hmpzl->SetMarkerSize(0.5);
 463     hmpzl->Draw("psame");
 464     hbrzl->SetLineColor(kGreen+1);
 465     hbrzl->SetMarkerColor(kGreen+1);
 466     hbrzl->SetMarkerStyle(20);
 467     hbrzl->SetMarkerSize(0.5);
 468     hbrzl->Draw("same");
 469     t1->Draw();
 470     t2->Draw();
 471     t3->Draw();
 472     ceff2->Update();
 473   }
 474
 475
 476   if(hClSizAn && hClSizTb){
 477     TCanvas* ccs=new TCanvas("ccs","Cluster Size",1200,600);
 478     ccs->Divide(2,1);
 479     ccs->cd(1);
 480     gPad->SetLogz();
 481     gPad->SetRightMargin(0.12);
 482     hClSizAn->Draw("colz");
 483     hClSizAn->GetXaxis()->SetTitle("Drift Time (ns)");
 484     hClSizAn->GetYaxis()->SetTitle("Cluster Size - Anodes");
 485     ccs->cd(2);
 486     gPad->SetLogz();
 487     gPad->SetRightMargin(0.12);
 488     hClSizTb->Draw("colz");
 489     hClSizTb->GetXaxis()->SetTitle("Drift Time (ns)");
 490     hClSizTb->GetYaxis()->SetTitle("Cluster Size - Time Bins");
 491     ccs->Update();
 492   }
 493
 494   TH1F* htimR=(TH1F*)l->FindObject("hDrTimRP");
 495   TH1F* htimT=(TH1F*)l->FindObject("hDrTimTPAll");
 496   TH1F* htimTe=(TH1F*)l->FindObject("hDrTimTPExtra");
 497   TH1F* htimTne=(TH1F*)l->FindObject("hDrTimTPNoExtra");
 498   htimR->Rebin(4);
 499   htimT->Rebin(4);
 500   htimTe->Rebin(4);
 501   htimTne->Rebin(4);
 502   htimR->SetLineWidth(2);
 503   htimT->SetLineWidth(2);
 504   htimTe->SetLineWidth(2);
 505   htimTne->SetLineWidth(2);
 506
 507   TCanvas* ctim=new TCanvas("ctim","DriftTime",1400,600);
 508   ctim->Divide(2,1);
 509   ctim->cd(1);
 510   htimR->Draw();
 511   htimR->GetYaxis()->SetTitleOffset(1.2);
 512   htimR->GetXaxis()->SetTitle("Drift Time (ns)");
 513   htimR->GetYaxis()->SetTitle("RecPoints");
 514   ctim->cd(2);
 515   htimT->Draw();
 516   htimTe->SetLineColor(2);
 517   htimTe->Draw("same");
 518   htimTne->SetLineColor(4);
 519   htimTne->Draw("same");
 520   htimT->GetXaxis()->SetTitle("Drift Time (ns)");
 521   htimT->GetYaxis()->SetTitle("TrackPoints");
 522   htimT->GetYaxis()->SetTitleOffset(1.2);
 523   TLatex* ta=new TLatex(0.5,0.85,"All Clusters");
 524   ta->SetNDC();
 525   ta->SetTextColor(1);
 526   ta->Draw();
 527   TLatex* te=new TLatex(0.5,0.8,"Extra Clusters");
 528   te->SetNDC();
 529   te->SetTextColor(2);
 530   te->Draw();
 531   TLatex* tn=new TLatex(0.5,0.75,"Non-Extra Clusters");
 532   tn->SetNDC();
 533   tn->SetTextColor(4);
 534   tn->Draw();
 535   ctim->Update();
 536
 537   TCanvas* cdedx=new TCanvas("cdedx","dedx",1400,600);
 538   cdedx->Divide(3,1);
 539   cdedx->cd(1);
 540   gPad->SetLogz();
 541   hdedx3->Draw("col");
 542   hdedx3->GetXaxis()->SetTitle("P (GeV/c)");
 543   hdedx3->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 3");
 544   hdedx3->GetYaxis()->SetTitleOffset(1.25);
 545   cdedx->cd(2);
 546   gPad->SetLogz();
 547   hdedx4->Draw("col");
 548   hdedx4->GetXaxis()->SetTitle("P (GeV/c)");
 549   hdedx4->GetYaxis()->SetTitle("dE/dx (keV/300 #mum) Layer 4");
 550   hdedx4->GetYaxis()->SetTitleOffset(1.25);
 551   cdedx->cd(3);
 552   gPad->SetLogz();
 553   hdedxmod->Draw("col");
 554   hdedxmod->GetXaxis()->SetTitle("SDD Module Id");
 555   hdedxmod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 556   hdedxmod->GetYaxis()->SetTitleOffset(1.25);
 557   cdedx->Update();
 558
 559   printf("---- dE/dx vs.DriftTime ----\n");
 560   TCanvas* csig=new TCanvas("csig","dedx vs. DriftTime",1000,700);
 561   csig->Divide(4,2);
 562   TH1F* hSigTim[8];
 563   TGraphErrors* gmpv=new TGraphErrors(0);
 564   TGraphErrors* gsigg=new TGraphErrors(0);
 565   TGraphErrors* gsigl=new TGraphErrors(0);
 566   gmpv->SetTitle("");
 567   gsigg->SetTitle("");
 568   gsigl->SetTitle("");
 569   Int_t iPoint=0;
 570   TF1 *lfun = new TF1("LangausFun",LangausFun,50.,300.,4);
 571   for(Int_t it=0; it<8; it++){
 572     hSigTim[it]=(TH1F*)l->FindObject(Form("hSigTimeInt%d",it));
 573     csig->cd(it+1);
 574     hSigTim[it]->Draw();
 575     if(hSigTim[it]->GetEntries()>200){
 576       lfun->SetLineWidth(2);
 577       lfun->SetParameter(0,5.);
 578       lfun->SetParameter(1,80.);
 579       lfun->SetParameter(2,hSigTim[it]->GetEntries()/10.);
 580       lfun->SetParameter(3,10.);
 581       lfun->SetParLimits(3,0.,20);
 582
 583       hSigTim[it]->Fit("LangausFun","QLR");
 584       hSigTim[it]->GetXaxis()->SetTitle(Form("dE/dx, time interval %d",it+1));
 585       hSigTim[it]->GetYaxis()->SetTitle("Events");
 586       Float_t mpv=lfun->GetParameter(1);
 587       Float_t empv=lfun->GetParError(1);
 588       Float_t sig=lfun->GetParameter(3);
 589       Float_t esig=lfun->GetParError(3);
 590       Float_t sigl=lfun->GetParameter(0);
 591       Float_t esigl=lfun->GetParError(0);
 592       gmpv->SetPoint(iPoint,(Float_t)it,mpv);
 593       gmpv->SetPointError(iPoint,0.,empv);
 594       gsigg->SetPoint(iPoint,(Float_t)it,sig);
 595       gsigg->SetPointError(iPoint,0.,esig);
 596       gsigl->SetPoint(iPoint,(Float_t)it,sigl);
 597       gsigl->SetPointError(iPoint,0.,esigl);
 598       ++iPoint;
 599       gPad->Update();
 600       printf("Bin %d - MPV=%.3f  \t SigmaLandau=%.3f  \t SigmaGaus=%.3f\n",it,mpv,sigl,sig);
 601     }
 602   }
 603
 604   TCanvas* cpars=new TCanvas("cpars","Params",800,900);
 605   cpars->Divide(1,3,0.01,0.);
 606   cpars->cd(1);
 607   gPad->SetLeftMargin(0.14);
 608   gPad->SetFrameLineWidth(2);
 609   gPad->SetTickx();
 610   gPad->SetTicky();
 611   gmpv->SetMarkerStyle(20);
 612   //  gmpv->SetMinimum(0);
 613   //  gmpv->SetMaximum(120);
 614   gmpv->GetXaxis()->SetLimits(-0.2,6.8);
 615   gmpv->Draw("AP");
 616   //  gmpv->GetXaxis()->SetTitle("Drift Time interval number");
 617   gmpv->GetYaxis()->SetTitle("Landau MPV (keV)");
 618   gmpv->GetXaxis()->SetTitleSize(0.05);
 619   gmpv->GetYaxis()->SetTitleSize(0.05);
 620   gmpv->GetYaxis()->SetTitleOffset(1.2);
 621   cpars->cd(2);
 622   gPad->SetLeftMargin(0.14);
 623   gPad->SetFrameLineWidth(2);
 624   gPad->SetTickx();
 625   gPad->SetTicky();
 626   gsigl->SetMarkerStyle(20);
 627   gsigl->GetXaxis()->SetLimits(-0.2,6.8);
 628   gsigl->Draw("AP");
 629   //  gsigl->GetXaxis()->SetTitle("Drift Time interval number");
 630   gsigl->GetYaxis()->SetTitle("#sigma_{Landau} (keV)");
 631   gsigl->GetXaxis()->SetTitleSize(0.05);
 632   gsigl->GetYaxis()->SetTitleSize(0.05);
 633   gsigl->GetYaxis()->SetTitleOffset(1.2);
 634   cpars->cd(3);
 635   gPad->SetLeftMargin(0.14);
 636   gPad->SetFrameLineWidth(2);
 637   gPad->SetTickx();
 638   gPad->SetTicky();
 639   gsigg->SetMarkerStyle(20);
 640   gsigg->GetXaxis()->SetLimits(-0.2,6.8);
 641   gsigg->Draw("AP");
 642   gsigg->GetXaxis()->SetTitle("Drift Time interval number");
 643   gsigg->GetYaxis()->SetTitle("#sigma_{Gauss} (keV)");
 644   gsigg->GetXaxis()->SetTitleSize(0.05);
 645   gsigg->GetYaxis()->SetTitleSize(0.05);
 646   gsigg->GetYaxis()->SetTitleOffset(1.2);
 647   cpars->Update();
 648
 649 }
 650
 651
 652 Double_t LangausFun(Double_t *x, Double_t *par) {
 653
 654   //Fit parameters:
 655   //par[0]=Width (scale) parameter of Landau density
 656   //par[1]=Most Probable (MP, location) parameter of Landau density
 657   //par[2]=Total area (integral -inf to inf, normalization constant)
 658   //par[3]=Width (sigma) of convoluted Gaussian function
 659   //
 660   //In the Landau distribution (represented by the CERNLIB approximation),
 661   //the maximum is located at x=-0.22278298 with the location parameter=0.
 662   //This shift is corrected within this function, so that the actual
 663   //maximum is identical to the MP parameter.
 664
 665   // Numeric constants
 666   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 667   Double_t mpshift  = -0.22278298;       // Landau maximum location
 668
 669   // Control constants
 670   Double_t np = 100.0;      // number of convolution steps
 671   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 672
 673   // Variables
 674   Double_t xx;
 675   Double_t mpc;
 676   Double_t fland;
 677   Double_t sum = 0.0;
 678   Double_t xlow,xupp;
 679   Double_t step;
 680   Double_t i;
 681
 682
 683   // MP shift correction
 684   mpc = par[1] - mpshift * par[0];
 685
 686   // Range of convolution integral
 687   xlow = x[0] - sc * par[3];
 688   xupp = x[0] + sc * par[3];
 689
 690   step = (xupp-xlow) / np;
 691
 692   // Convolution integral of Landau and Gaussian by sum
 693   for(i=1.0; i<=np/2; i++) {
 694     xx = xlow + (i-.5) * step;
 695     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 696     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 697
 698     xx = xupp - (i-.5) * step;
 699     fland = TMath::Landau(xx,mpc,par[0]) / par[0];
 700     sum += fland * TMath::Gaus(x[0],xx,par[3]);
 701   }
 702
 703   return (par[2] * step * sum * invsq2pi / par[3]);
 704 }
 705