ITS/ShowITSRecPoints.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TCanvas.h>
   3 #include <TClassTable.h>
   4 #include <TGraph.h>
   5 #include <TGraph2D.h>
   6 #include <TGeoManager.h>
   7 #include <TH1.h>
   8 #include <TH2.h>
   9 #include <TLatex.h>
  10 #include <TStyle.h>
  11 #include <TInterpreter.h>
  12 #include "AliGeomManager.h"
  13 #include "AliHeader.h"
  14 #include "AliITS.h"
  15 #include "AliITSDetTypeRec.h"
  16 #include "AliITSgeomTGeo.h"
  17 #include "AliITSRecPoint.h"
  18 #include "AliRun.h"
  19 #endif
  20
  21 Int_t ShowITSRecPoints(Int_t nevfordisp=0){
  22   ///////////////////////////////////////////////////////////////////////
  23   // Macro to check clusters in the 6 ITS layers                       //
  24   // Provides:                                                         //
  25   //  6 canvases with 9 plots each (1 canvas for each layer)           //
  26   //  3 canvases with cluster XY coordinates for the first 3 events    //
  27   ///////////////////////////////////////////////////////////////////////
  28
  29   if (gClassTable->GetID("AliRun") < 0) {
  30     gInterpreter->ExecuteMacro("loadlibs.C");
  31   }
  32
  33   // retrives geometry
  34   if(!gGeoManager){
  35     AliGeomManager::LoadGeometry("geometry.root");
  36   }
  37
  38   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  39   if (rl == 0x0){
  40     cerr<<"Can not open session RL=NULL"<< endl;
  41     return -1;
  42   }
  43
  44   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  45   if(!ITSloader){
  46     cerr<<"ITS loader not found"<<endl;
  47     return -1;
  48   }
  49   ITSloader->LoadRecPoints("read");
  50
  51   Float_t cluglo[3]={0.,0.,0.};
  52   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  53
  54   Int_t totmod=AliITSgeomTGeo::GetNModules();
  55   Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
  56   Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
  57
  58   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  59   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
  60
  61   TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
  62   TH1F** hmod=new TH1F*[6];
  63   TH1F** hxl=new TH1F*[6];
  64   TH1F** hzl=new TH1F*[6];
  65   TH1F** hxg=new TH1F*[6];
  66   TH1F** hyg=new TH1F*[6];
  67   TH1F** hzg=new TH1F*[6];
  68   TH1F** hr=new TH1F*[6];
  69   TH1F** hphi=new TH1F*[6];
  70   TH2F** hzphi=new TH2F*[6];
  71   TH1F** hq=new TH1F*[6];
  72   TH1F** hdrtime=new TH1F*[6];
  73
  74   Char_t name[10];
  75   for(Int_t iLay=0;iLay<6;iLay++){
  76     sprintf(name,"hmod%d",iLay+1);
  77     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
  78     hmod[iLay]->GetXaxis()->SetTitle("Module");
  79     hmod[iLay]->GetXaxis()->CenterTitle();
  80     sprintf(name,"hxloc%d",iLay+1);
  81     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
  82     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
  83     hxl[iLay]->GetXaxis()->CenterTitle();
  84     sprintf(name,"hzloc%d",iLay+1);
  85     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
  86     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
  87     hzl[iLay]->GetXaxis()->CenterTitle();
  88     sprintf(name,"hxgl%d",iLay+1);
  89     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
  90     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
  91     hxg[iLay]->GetXaxis()->CenterTitle();
  92     sprintf(name,"hygl%d",iLay+1);
  93     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
  94     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
  95     hyg[iLay]->GetXaxis()->CenterTitle();
  96     sprintf(name,"hzgl%d",iLay+1);
  97     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
  98     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
  99     hzg[iLay]->GetXaxis()->CenterTitle();
 100     sprintf(name,"hr%d",iLay+1);
 101     hr[iLay]=new TH1F(name,"",100,0.,50.);
 102     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
 103     hr[iLay]->GetXaxis()->CenterTitle();
 104     sprintf(name,"hphi%d",iLay+1);
 105     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
 106     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 107     hphi[iLay]->GetXaxis()->CenterTitle();
 108     sprintf(name,"hq%d",iLay+1);
 109     hq[iLay]=new TH1F(name,"",100,0.,300.);
 110     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
 111     hq[iLay]->GetXaxis()->CenterTitle();
 112     sprintf(name,"hdtim%d",iLay+1);
 113     hdrtime[iLay]=new TH1F(name,"",100,0.,7000.);
 114     hdrtime[iLay]->GetXaxis()->SetTitle("Drift Time (ns)");
 115     hdrtime[iLay]->GetXaxis()->CenterTitle();
 116     sprintf(name,"hzphi%d",iLay+1);
 117     hzphi[iLay]=new TH2F(name,Form("Layer %d",iLay+1),50,-TMath::Pi(),TMath::Pi(),50,-zlim[iLay],zlim[iLay]);
 118     hzphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 119     hzphi[iLay]->GetYaxis()->SetTitle("Zglob (cm)");
 120     hzphi[iLay]->SetStats(0);
 121   }
 122
 123   TGraph *gptsXY=new TGraph(0);
 124   TGraph *gptsRZ=new TGraph(0);
 125
 126   Int_t totev=rl->GetNumberOfEvents();
 127   printf("Total Number of events = %d\n",totev);
 128
 129   for(Int_t iev=0;iev<totev;iev++){
 130     rl->GetEvent(iev);
 131     TTree *TR = ITSloader->TreeR();
 132     TClonesArray *ITSrec  = detTypeRec->RecPoints();
 133     TBranch *branch = 0;
 134     if(TR && ITSrec){
 135       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 136       if(branch)branch->SetAddress(&ITSrec);
 137     }
 138     if(iev%100==0) printf("Event #%d\n",iev);
 139
 140
 141     Int_t ipt=0;
 142     for (Int_t mod=modmin; mod<=modmax; mod++){
 143       detTypeRec->ResetRecPoints();
 144       branch->GetEvent(mod);
 145       Int_t nrecp = ITSrec->GetEntries();
 146       if(nrecp>0){
 147         for(Int_t irec=0;irec<nrecp;irec++) {
 148           AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 149           Int_t lay=recp->GetLayer();
 150           hlayer->Fill(lay);
 151           recp->GetGlobalXYZ(cluglo);
 152           Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]);
 153           Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
 154           if(iev==nevfordisp){
 155             gptsXY->SetPoint(ipt,cluglo[0],cluglo[1]);
 156             if(cluglo[1]>0) gptsRZ->SetPoint(ipt,cluglo[2],rad);
 157             else gptsRZ->SetPoint(ipt,cluglo[2],-rad);
 158             ipt++;
 159           }
 160           hmod[lay]->Fill(mod);
 161           hzl[lay]->Fill(recp->GetDetLocalZ());
 162           hxl[lay]->Fill(recp->GetDetLocalX());
 163           hzg[lay]->Fill(cluglo[2]);
 164           hyg[lay]->Fill(cluglo[1]);
 165           hxg[lay]->Fill(cluglo[0]);
 166           hr[lay]->Fill(rad);
 167           hphi[lay]->Fill(phi);
 168           hq[lay]->Fill(recp->GetQ());
 169           hdrtime[lay]->Fill(recp->GetDriftTime());
 170           hzphi[lay]->Fill(phi,cluglo[2]);
 171         }
 172       }
 173     }
 174   }
 175
 176   TCanvas **c=new TCanvas*[6];
 177   Char_t ctit[30];
 178   for(Int_t iLay=0;iLay<6;iLay++){
 179     sprintf(name,"can%d",iLay+1);
 180     sprintf(ctit,"Layer %d",iLay+1);
 181     c[iLay]=new TCanvas(name,ctit,1200,900);
 182     c[iLay]->Divide(3,3,0.001,0.001);
 183     c[iLay]->cd(1);
 184     hmod[iLay]->Draw();
 185     c[iLay]->cd(2);
 186     hxl[iLay]->Draw();
 187     c[iLay]->cd(3);
 188     hzl[iLay]->Draw();
 189     c[iLay]->cd(4);
 190     hxg[iLay]->Draw();
 191     c[iLay]->cd(5);
 192     hyg[iLay]->Draw();
 193     c[iLay]->cd(6);
 194     hzg[iLay]->Draw();
 195     c[iLay]->cd(7);
 196     hr[iLay]->Draw();
 197     c[iLay]->cd(8);
 198     hphi[iLay]->Draw();
 199     c[iLay]->cd(9);
 200     hq[iLay]->Draw();
 201   }
 202
 203   TCanvas* cdrtim=new TCanvas("cdrtim","SDD drift time",1000,600);
 204   cdrtim->Divide(2,1);
 205   cdrtim->cd(1);
 206   hdrtime[2]->Draw();
 207   cdrtim->cd(2);
 208   hdrtime[3]->Draw();
 209
 210   gStyle->SetPalette(1);
 211   TLatex* tstat=new TLatex();
 212   tstat->SetNDC();
 213   TCanvas* cspd=new TCanvas("cspd","SPD",1000,600);
 214   cspd->Divide(2,1);
 215   cspd->cd(1);
 216   hzphi[0]->Draw("colz");
 217   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[0]->GetEntries())));
 218   cspd->cd(2);
 219   hzphi[1]->Draw("colz");
 220   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[1]->GetEntries())));
 221
 222   TCanvas* csdd=new TCanvas("csdd","SDD",1000,600);
 223   csdd->Divide(2,1);
 224   csdd->cd(1);
 225   hzphi[2]->Draw("colz");
 226   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[2]->GetEntries())));
 227   csdd->cd(2);
 228   hzphi[3]->Draw("colz");
 229   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[3]->GetEntries())));
 230
 231   TCanvas* cssd=new TCanvas("cssd","SSD",1000,600);
 232   cssd->Divide(2,1);
 233   cssd->cd(1);
 234   hzphi[4]->Draw("colz");
 235   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[4]->GetEntries())));
 236   cssd->cd(2);
 237   hzphi[5]->Draw("colz");
 238   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[5]->GetEntries())));
 239
 240   TCanvas *cev0;
 241   sprintf(ctit,"Event %d XY",nevfordisp);
 242   cev0=new TCanvas("cev0",ctit,600,600);
 243   if(gptsXY->GetN()>0){
 244     gptsXY->SetMarkerStyle(7);
 245     gptsXY->SetTitle(0);
 246     gptsXY->Draw("AP");
 247     gptsXY->GetXaxis()->SetTitle("Xglob");
 248     gptsXY->GetYaxis()->SetTitle("Yglob");
 249    }
 250
 251   TCanvas *cev1;
 252   sprintf(ctit,"Event %d Zr",nevfordisp);
 253   cev1=new TCanvas("cev1",ctit,600,600);
 254   if(gptsRZ->GetN()>0){
 255     gptsRZ->SetMarkerStyle(7);
 256     gptsRZ->SetTitle(0);
 257     gptsRZ->Draw("AP");
 258     gptsRZ->GetXaxis()->SetTitle("Zglob");
 259     gptsRZ->GetYaxis()->SetTitle("Radius");
 260   }
 261
 262   return 0;
 263 }
 264
 265