ITS/macrosSDD/CheckSDDInESD.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 <TMath.h>
  11 #include <TCanvas.h>
  12 #include <TStyle.h>
  13 #include <TLatex.h>
  14 #include "AliITSgeomTGeo.h"
  15 #include "AliESDEvent.h"
  16 #endif
  17
  18 enum {kAll, kTPCITS, kITSsa, kITSpureSA};
  19
  20 void CheckSDDInESD(TString filename="AliESDs.root", Int_t optTracks=kAll){
  21
  22
  23   TFile* esdFile = TFile::Open(filename.Data());
  24   if (!esdFile || !esdFile->IsOpen()) {
  25     printf("Error in opening ESD file");
  26     return;
  27   }
  28
  29   AliESDEvent * esd = new AliESDEvent;
  30   TTree* tree = (TTree*) esdFile->Get("esdTree");
  31   if (!tree) {
  32     printf("Error: no ESD tree found");
  33     return;
  34   }
  35   esd->ReadFromTree(tree);
  36   TH1F* hpt=new TH1F("hpt","",100,0.,10.);
  37   TH1F* hphi=new TH1F("hphi","",100,-1,1);
  38   TH1F* hlam=new TH1F("hlam","",100,-2.,2.);
  39   TH1F* halpha=new TH1F("halpha","",100,-7,7);
  40   TH1F* hitscl=new TH1F("hitscl","",7,-0.5,6.5);
  41   TH1F* htpccl=new TH1F("htpccl","",200,-0.5,199.5);
  42   TH1F* hitsmap=new TH1F("hitsmap","",64,-0.5,63.5);
  43   TH1F* hclulay=new TH1F("hclulay","",7,-1.5,5.5);
  44
  45   TH1F* hvx=new TH1F("hvx","",100,-1.,1.);
  46   TH1F* hvy=new TH1F("hvy","",100,-1.,1.);
  47   TH1F* hvz=new TH1F("hvz","",100,-20.,20.);
  48   TH1F* hdedx3=new TH1F("hdedx3","",100,0.,300.);
  49   TH1F* hdedx4=new TH1F("hdedx4","",100,0.,300.);
  50   TH1F* hdedx5=new TH1F("hdedx5","",100,0.,300.);
  51   TH1F* hdedx6=new TH1F("hdedx6","",100,0.,300.);
  52   TH1F* hStatus=new TH1F("hStatus","",11,-1.5,9.5);
  53
  54
  55   // -- Local coordinates
  56
  57
  58   // -- Module histos
  59
  60   TH1F* hAllPMod  = new TH1F("hAllPmod","Crossing Tracks vs. Module",260,239.5,499.5);
  61   TH1F* hGoodPMod  = new TH1F("hGoodPmod","PointsAssocToTrack per Module",260,239.5,499.5);
  62   TH1F* hBadRegMod  = new TH1F("hBadRegmod","Tracks in BadRegion per Module",260,239.5,499.5);
  63   TH1F* hMissPMod  = new TH1F("hMissPmod","Missing Points per Module",260,239.5,499.5);
  64   TH1F* hSkippedMod  = new TH1F("hSkippedmod","Tracks in Skipped Module",260,239.5,499.5);
  65   TH1F* hOutAccMod  = new TH1F("hOutAccmod","Tracks outside zAcc per Module",260,239.5,499.5);
  66   TH1F* hNoRefitMod  = new TH1F("hNoRefitmod","Points rejected in refit per Module",260,239.5,499.5);
  67
  68   TH1F* hAllPXloc  = new TH1F("hAllPxloc","Crossing Tracks vs. Xloc",75, -3.75, 3.75);
  69   TH1F* hGoodPXloc  = new TH1F("hGoodPxloc","PointsAssocToTrack vs. Xloc",75, -3.75, 3.75);
  70   TH1F* hBadRegXloc  = new TH1F("hBadRegxloc","Tracks in BadRegion vs. Xloc",75, -3.75, 3.75);
  71   TH1F* hMissPXloc  = new TH1F("hMissPxloc","Missing Points vs. Xloc",75, -3.75, 3.75);
  72   TH1F* hAllPZloc  = new TH1F("hAllPzloc","Crossing Tracks vs. Zloc",77, -3.85, 3.85);
  73   TH1F* hGoodPZloc  = new TH1F("hGoodPzloc","PointsAssocToTrack vs. Zloc",77, -3.85, 3.85);
  74   TH1F* hBadRegZloc  = new TH1F("hBadRegzloc","Tracks in BadRegion vs. Zloc",77, -3.85, 3.85);
  75   TH1F* hMissPZloc  = new TH1F("hMissPzloc","Missing Points vs. Zloc",77, -3.85, 3.85);
  76   TH2F* hdEdxVsMod=new TH2F("hdEdxVsMod","dE/dx vs. mod",260,239.5,499.5,100,0.,500.);
  77
  78   gStyle->SetPalette(1);
  79
  80
  81   for (Int_t iEvent = 0; iEvent < tree->GetEntries(); iEvent++) {
  82     tree->GetEvent(iEvent);
  83     if (!esd) {
  84       printf("Error: no ESD object found for event %d", iEvent);
  85       return;
  86     }
  87     cout<<"-------- Event "<<iEvent<<endl;
  88     printf(" Tracks # = %d\n",esd->GetNumberOfTracks());
  89     const AliESDVertex *spdv=esd->GetVertex();
  90     printf(" SPD Primary Vertex in %f %f %f with %d contributors\n",spdv->GetX(),spdv->GetY(),spdv->GetZ(),spdv->GetNContributors());
  91     const AliESDVertex *trkv=esd->GetPrimaryVertex();
  92     printf(" Track Primary Vertex with %d contributors\n",trkv->GetNContributors());
  93     if(spdv->IsFromVertexer3D()){
  94       hvx->Fill(spdv->GetX());
  95       hvy->Fill(spdv->GetY());
  96       hvz->Fill(spdv->GetZ());
  97     }
  98     Double_t itss[4];
  99     for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
 100       AliESDtrack* track = esd->GetTrack(iTrack);
 101       Int_t nITSclus=track->GetNcls(0);
 102       UChar_t clumap=track->GetITSClusterMap();
 103       Int_t nPointsForPid=0;
 104       for(Int_t i=2; i<6; i++){
 105         if(clumap&(1<<i)) ++nPointsForPid;
 106       }
 107       //      track->PropagateTo(4.,5.);
 108       htpccl->Fill(track->GetNcls(1));
 109       Int_t status=track->GetStatus();
 110       Bool_t tpcin=0;
 111       hStatus->Fill(-1.);
 112       if(status & AliESDtrack::kTPCin){
 113         tpcin=1;
 114         hStatus->Fill(0.);
 115       }
 116       if(status & AliESDtrack::kTPCout){
 117         hStatus->Fill(1.);
 118       }
 119       if(status & AliESDtrack::kTPCrefit){
 120         hStatus->Fill(2.);
 121       }
 122       Bool_t itsin=0;
 123       if(status & AliESDtrack::kITSin){
 124         itsin=1;
 125         hStatus->Fill(3.);
 126       }
 127       if(status & AliESDtrack::kITSout){
 128         hStatus->Fill(4.);
 129       }
 130       if(status & AliESDtrack::kITSrefit){
 131         hStatus->Fill(5.);
 132       }
 133       if(!tpcin && itsin){
 134         hStatus->Fill(6.);
 135       }
 136       if(status & AliESDtrack::kITSpureSA){
 137         hStatus->Fill(7.);
 138       }
 139
 140       if(status & AliESDtrack::kITSrefit){
 141         if((optTracks==kTPCITS) && !(status & AliESDtrack::kTPCin)) continue;
 142         if((optTracks==kITSsa) && (status & AliESDtrack::kTPCin)) continue;
 143         if((optTracks==kITSsa) && (status & AliESDtrack::kITSpureSA)) continue;
 144         if((optTracks==kITSpureSA) && (status & AliESDtrack::kITSpureSA)) continue;
 145
 146          track->GetITSdEdxSamples(itss);
 147         //      printf("Track %d (label %d) in ITS with %d clusters clumap %d pointspid= %d\n",iTrack,track->GetLabel(),nITSclus,clumap,nPointsForPid);
 148         //printf("   dedx=%f %f %f %f\n",itss[0],itss[1],itss[2],itss[3]);
 149         hitscl->Fill(nITSclus);
 150         hdedx3->Fill(itss[0]);
 151         hdedx4->Fill(itss[1]);
 152         hdedx5->Fill(itss[2]);
 153         hdedx6->Fill(itss[3]);
 154         hitsmap->Fill(clumap);
 155         hclulay->Fill(-1.);
 156         for(Int_t iLay=0;iLay<6;iLay++){
 157           if(clumap&1<<iLay) hclulay->Fill(iLay);
 158         }
 159         hpt->Fill(track->Pt());
 160         hphi->Fill(TMath::ASin(track->GetSnp()));
 161         hlam->Fill(TMath::ATan(track->GetTgl()));
 162         halpha->Fill(track->GetAlpha());
 163         Int_t iMod,status;
 164         Float_t xloc,zloc;
 165         for(Int_t iLay=2; iLay<=3; iLay++){
 166           Bool_t ok=track->GetITSModuleIndexInfo(iLay,iMod,status,xloc,zloc);
 167           if(ok){
 168             iMod+=240;
 169             hAllPMod->Fill(iMod);
 170             hAllPXloc->Fill(xloc);
 171             hAllPZloc->Fill(zloc);
 172             if(status==1){
 173               hGoodPMod->Fill(iMod);
 174               hGoodPXloc->Fill(xloc);
 175               hGoodPZloc->Fill(zloc);
 176               if(track->Pt()>1.) hdEdxVsMod->Fill(iMod,itss[iLay-2]);
 177             }
 178             else if(status==2){
 179               hBadRegMod->Fill(iMod);
 180               hBadRegXloc->Fill(xloc);
 181               hBadRegZloc->Fill(zloc);
 182             }
 183             else if(status==3) hSkippedMod->Fill(iMod);
 184             else if(status==4) hOutAccMod->Fill(iMod);
 185             else if(status==5){
 186               hMissPMod->Fill(iMod);
 187               hMissPXloc->Fill(xloc);
 188               hMissPZloc->Fill(zloc);
 189             }
 190             else if(status==6) hNoRefitMod->Fill(iMod);
 191           }
 192         }
 193       }
 194     }
 195   }
 196   Float_t norm=hclulay->GetBinContent(1);
 197   if(norm<1.) norm=1.;
 198   hclulay->Scale(1./norm);
 199   gStyle->SetLineWidth(2);
 200
 201   TCanvas* c1=new TCanvas("c1","Track quantities",900,900);
 202   c1->Divide(2,2);
 203   c1->cd(1);
 204   htpccl->Draw();
 205   htpccl->GetXaxis()->SetTitle("Clusters in TPC ");
 206   c1->cd(2);
 207   hitscl->Draw();
 208   hitscl->GetXaxis()->SetTitle("Clusters in ITS ");
 209   c1->cd(3);
 210   hclulay->Draw();
 211   hclulay->GetXaxis()->SetRange(2,7);
 212   hclulay->GetXaxis()->SetTitle("# ITS Layer");
 213   hclulay->GetYaxis()->SetTitle("Fraction of tracks with point in Layer x");
 214   c1->cd(4);
 215
 216   TCanvas* c2=new TCanvas("c2","dedx per Layer",900,900);
 217   c2->Divide(2,2);
 218   c2->cd(1);
 219   hdedx3->Draw();
 220   hdedx3->GetXaxis()->SetTitle("dE/dx Lay3");
 221   c2->cd(2);
 222   hdedx4->Draw();
 223   hdedx4->GetXaxis()->SetTitle("dE/dx Lay4");
 224   c2->cd(3);
 225   hdedx5->Draw();
 226   hdedx5->GetXaxis()->SetTitle("dE/dx Lay5");
 227   c2->cd(4);
 228   hdedx6->Draw();
 229   hdedx6->GetXaxis()->SetTitle("dE/dx Lay6");
 230
 231   hdEdxVsMod->SetStats(0);
 232   TCanvas* cdedx=new TCanvas("cdedx","dedx SDD",1400,600);
 233   cdedx->SetLogz();
 234   hdEdxVsMod->Draw("col");
 235   hdEdxVsMod->GetXaxis()->SetTitle("SDD Module Id");
 236   hdEdxVsMod->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 237   hdEdxVsMod->GetYaxis()->SetTitleOffset(1.25);
 238
 239
 240
 241   TCanvas* cv=new TCanvas("cv","Vertex",600,900);
 242   cv->Divide(1,3);
 243   cv->cd(1);
 244   hvx->Draw();
 245   hvx->GetXaxis()->SetTitle("Xv (cm)");
 246   cv->cd(2);
 247   hvy->Draw();
 248   hvy->GetXaxis()->SetTitle("Yv (cm)");
 249   cv->cd(3);
 250   hvz->Draw();
 251   hvz->GetXaxis()->SetTitle("Xv (cm)");
 252
 253   hGoodPMod->SetStats(0);
 254   hGoodPMod->SetTitle("");
 255   TCanvas* ceff0=new TCanvas("ceff0","ModuleIndexInfo",1000,600);
 256   hGoodPMod->Draw("e");
 257   hGoodPMod->GetXaxis()->SetTitle("SDD Module Id");
 258   hGoodPMod->GetYaxis()->SetTitle("Number of tracks");
 259   hMissPMod->SetLineColor(2);
 260   hMissPMod->SetMarkerColor(2);
 261   hMissPMod->SetMarkerStyle(22);
 262   hMissPMod->SetMarkerSize(0.5);
 263   hMissPMod->Draw("psame");
 264   hBadRegMod->SetLineColor(kGreen+1);
 265   hBadRegMod->SetMarkerColor(kGreen+1);
 266   hBadRegMod->SetMarkerStyle(20);
 267   hBadRegMod->SetMarkerSize(0.5);
 268   hBadRegMod->Draw("esame");
 269   hSkippedMod->SetLineColor(kYellow);
 270   hSkippedMod->Draw("esame");
 271   hOutAccMod->SetLineColor(4);
 272   hOutAccMod->Draw("esame");
 273   hNoRefitMod->SetLineColor(6);
 274   hNoRefitMod->Draw("esame");
 275   TLatex* t1=new TLatex(0.7,0.85,"Good Point");
 276   t1->SetNDC();
 277   t1->SetTextColor(1);
 278   t1->Draw();
 279   TLatex* t2=new TLatex(0.7,0.8,"Missing Point");
 280   t2->SetNDC();
 281   t2->SetTextColor(2);
 282   t2->Draw();
 283   TLatex* t3=new TLatex(0.7,0.75,"Bad Region");
 284   t3->SetNDC();
 285   t3->SetTextColor(kGreen+1);
 286   t3->Draw();
 287   ceff0->Update();
 288
 289   TH1F* heff=new TH1F("heff","",260,239.5,499.5);
 290   for(Int_t imod=0; imod<260;imod++){
 291     Float_t numer=hGoodPMod->GetBinContent(imod+1)+hBadRegMod->GetBinContent(imod+1)+hOutAccMod->GetBinContent(imod+1)+hNoRefitMod->GetBinContent(imod+1);
 292     Float_t denom=hAllPMod->GetBinContent(imod+1);
 293     Float_t eff=0.;
 294     Float_t erreff=0.;
 295     if(denom>0){
 296       eff=numer/denom;
 297       erreff=TMath::Sqrt(eff*(1-eff)/denom);
 298     }
 299     heff->SetBinContent(imod+1,eff);
 300     heff->SetBinError(imod+1,erreff);
 301   }
 302
 303   printf("---- Modules with efficiency < 90%% ----\n");
 304   heff->SetStats(0);
 305   TCanvas* ceff1=new TCanvas("ceff1","Efficiency",1000,600);
 306   heff->Draw();
 307   heff->GetXaxis()->SetTitle("SDD Module Id");
 308   heff->GetYaxis()->SetTitle("Fraction of tracks with point in good region");
 309   for(Int_t ibin=1; ibin<=heff->GetNbinsX(); ibin++){
 310     Float_t e=heff->GetBinContent(ibin);
 311     if(e<0.9){
 312       Int_t iMod=(Int_t)heff->GetBinCenter(ibin);
 313       Int_t lay,lad,det;
 314       AliITSgeomTGeo::GetModuleId(iMod,lay,lad,det);
 315       printf("Module %d - Layer %d Ladder %2d Det %d  -   Eff. %.3f\n",iMod,lay,lad,det,heff->GetBinContent(ibin));
 316     }
 317   }
 318   ceff1->Update();
 319
 320
 321
 322   hGoodPXloc->SetTitle("");
 323   hGoodPZloc->SetTitle("");
 324   hGoodPXloc->SetStats(0);
 325   hGoodPZloc->SetStats(0);
 326   hGoodPXloc->SetMinimum(0);
 327   hGoodPZloc->SetMinimum(0);
 328   TCanvas* ceff2=new TCanvas("ceff2","LocalCoord",1000,600);
 329   ceff2->Divide(2,1);
 330   ceff2->cd(1);
 331   hGoodPXloc->Draw("e");
 332   hGoodPXloc->GetXaxis()->SetTitle("Xlocal (cm)");
 333   hGoodPXloc->GetYaxis()->SetTitle("Number of tracks");
 334   hMissPXloc->SetLineColor(2);
 335   hMissPXloc->SetMarkerColor(2);
 336   hMissPXloc->SetMarkerStyle(22);
 337   hMissPXloc->SetMarkerSize(0.5);
 338   hMissPXloc->Draw("psame");
 339   hBadRegXloc->SetLineColor(kGreen+1);
 340   hBadRegXloc->SetMarkerColor(kGreen+1);
 341   hBadRegXloc->SetMarkerStyle(20);
 342   hBadRegXloc->SetMarkerSize(0.5);
 343   hBadRegXloc->Draw("psame");
 344   t1->Draw();
 345   t2->Draw();
 346   t3->Draw();
 347   ceff2->cd(2);
 348   hGoodPZloc->Draw("e");
 349   hGoodPZloc->GetXaxis()->SetTitle("Zlocal (cm)");
 350   hGoodPZloc->GetYaxis()->SetTitle("Number of tracks");
 351   hMissPZloc->SetLineColor(2);
 352   hMissPZloc->SetMarkerColor(2);
 353   hMissPZloc->SetMarkerStyle(22);
 354   hMissPZloc->SetMarkerSize(0.5);
 355   hMissPZloc->Draw("psame");
 356   hBadRegZloc->SetLineColor(kGreen+1);
 357   hBadRegZloc->SetMarkerColor(kGreen+1);
 358   hBadRegZloc->SetMarkerStyle(20);
 359   hBadRegZloc->SetMarkerSize(0.5);
 360   hBadRegZloc->Draw("psame");
 361   t1->Draw();
 362   t2->Draw();
 363   t3->Draw();
 364   ceff2->Update();
 365
 366
 367
 368 }
 369
 370
 371
 372
 373