ITS/PlotITSHits.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include<Riostream.h>
   3 #include<TROOT.h>
   4 #include<TArrayI.h>
   5 #include<TBranch.h>
   6 #include<TCanvas.h>
   7 #include<TClassTable.h>
   8 #include<TClonesArray.h>
   9 #include<TFile.h>
  10 #include<TStyle.h>
  11 #include<TH1.h>
  12 #include<TH2.h>
  13 #include<TLatex.h>
  14 #include <TInterpreter.h>
  15 #include <TGeoManager.h>
  16 #include<TObject.h>
  17 #include<TObjArray.h>
  18 #include<TTree.h>
  19 #include<TNtuple.h>
  20 #include<TParticle.h>
  21 #include "AliStack.h"
  22 #include "AliRun.h"
  23 #include "AliCDBManager.h"
  24 #include "AliGeomManager.h"
  25 #include "AliITS.h"
  26 #include "AliITSgeomTGeo.h"
  27 #include "AliITSDetTypeRec.h"
  28 #include "AliITSRecPoint.h"
  29 #include "AliITSRecPoint.h"
  30 #include "AliITSdigit.h"
  31 #include "AliITSdigitSSD.h"
  32 #include "AliITShit.h"
  33 #include "AliITSmodule.h"
  34 #include "AliITSsegmentation.h"
  35 #include "AliITSsegmentationSPD.h"
  36 #include "AliITSsegmentationSDD.h"
  37 #include "AliITSsegmentationSSD.h"
  38 #include "AliRunLoader.h"
  39 #include "AliITSLoader.h"
  40 #include "AliHeader.h"
  41 #endif
  42
  43 /*  $Id$    */
  44
  45 // macro to display the coordindates (local+global) and the energy deposit
  46 // of ITS hits
  47
  48
  49 void PlotITSHits() {
  50
  51
  52   AliCDBManager* man = AliCDBManager::Instance();
  53   if (!man->IsDefaultStorageSet()) {
  54     printf("Setting a local default storage and run number 0\n");
  55     man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  56     man->SetRun(0);
  57   }
  58   else {
  59     printf("Using deafult storage \n");
  60   }
  61   // retrives geometry
  62   if(!gGeoManager){
  63     AliGeomManager::LoadGeometry("geometry.root");
  64   }
  65   AliGeomManager::ApplyAlignObjsFromCDB("ITS");
  66
  67   Int_t totmod=AliITSgeomTGeo::GetNModules();
  68   Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
  69   Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
  70
  71   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  72   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
  73
  74   TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
  75   TH1F** hmod=new TH1F*[6];
  76   TH1F** hxl=new TH1F*[6];
  77   TH1F** hzl=new TH1F*[6];
  78   TH1F** hxg=new TH1F*[6];
  79   TH1F** hyg=new TH1F*[6];
  80   TH1F** hzg=new TH1F*[6];
  81   TH1F** hr=new TH1F*[6];
  82   TH1F** hphi=new TH1F*[6];
  83   TH2F** hzphi=new TH2F*[6];
  84   TH1F** hq=new TH1F*[6];
  85
  86   Char_t name[10];
  87   for(Int_t iLay=0;iLay<6;iLay++){
  88     sprintf(name,"hmod%d",iLay+1);
  89     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
  90     hmod[iLay]->GetXaxis()->SetTitle("Module");
  91     hmod[iLay]->GetXaxis()->CenterTitle();
  92     sprintf(name,"hxloc%d",iLay+1);
  93     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
  94     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
  95     hxl[iLay]->GetXaxis()->CenterTitle();
  96     sprintf(name,"hzloc%d",iLay+1);
  97     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
  98     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
  99     hzl[iLay]->GetXaxis()->CenterTitle();
 100     sprintf(name,"hxgl%d",iLay+1);
 101     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
 102     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
 103     hxg[iLay]->GetXaxis()->CenterTitle();
 104     sprintf(name,"hygl%d",iLay+1);
 105     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
 106     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
 107     hyg[iLay]->GetXaxis()->CenterTitle();
 108     sprintf(name,"hzgl%d",iLay+1);
 109     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
 110     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
 111     hzg[iLay]->GetXaxis()->CenterTitle();
 112     sprintf(name,"hr%d",iLay+1);
 113     hr[iLay]=new TH1F(name,"",100,0.,50.);
 114     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
 115     hr[iLay]->GetXaxis()->CenterTitle();
 116     sprintf(name,"hphi%d",iLay+1);
 117     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
 118     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 119     hphi[iLay]->GetXaxis()->CenterTitle();
 120     sprintf(name,"hq%d",iLay+1);
 121     hq[iLay]=new TH1F(name,"",100,0.,300.);
 122     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
 123     hq[iLay]->GetXaxis()->CenterTitle();
 124     sprintf(name,"hzphi%d",iLay+1);
 125     hzphi[iLay]=new TH2F(name,Form("Layer %d",iLay+1),50,-TMath::Pi(),TMath::Pi(),50,-zlim[iLay],zlim[iLay]);
 126     hzphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 127     hzphi[iLay]->GetYaxis()->SetTitle("Zglob (cm)");
 128     hzphi[iLay]->SetStats(0);
 129   }
 130
 131
 132   AliRunLoader* rl = AliRunLoader::Open("galice.root");
 133   if (rl == 0x0){
 134     cerr<<"Can not open session RL=NULL"<< endl;
 135     return;
 136   }
 137   Int_t retval = rl->LoadgAlice();
 138   if (retval){
 139     cerr<<"LoadgAlice returned error"<<endl;
 140     return;
 141   }
 142   gAlice=rl->GetAliRun();
 143
 144   retval = rl->LoadHeader();
 145   if (retval){
 146     cerr<<"LoadHeader returned error"<<endl;
 147     return;
 148   }
 149
 150   retval = rl->LoadKinematics();
 151   if (retval){
 152     cerr<<"LoadKinematics returned error"<<endl;
 153     return;
 154   }
 155
 156   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
 157   if(!ITSloader){
 158     cerr<<"ITS loader not found"<<endl;
 159     return;
 160   }
 161
 162   ITSloader->LoadHits("read");
 163   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 164   ITS->SetTreeAddress();
 165   Int_t totev=rl->GetNumberOfEvents();
 166
 167   for(Int_t iev=0; iev<totev; iev++){
 168     rl->GetEvent(iev);
 169
 170     // HITS
 171     TTree *TH = ITSloader->TreeH();
 172     printf("Event %d  Tracks %d\n",iev,(Int_t)TH->GetEntries());
 173
 174     // ITS
 175     Int_t nmodules;
 176     ITS->InitModules(-1,nmodules);
 177     ITS->FillModules(0,0,nmodules," "," ");
 178
 179
 180     for (Int_t mod=modmin; mod<=modmax; mod++){
 181       Int_t lay,lad,det;
 182       AliITSgeomTGeo::GetModuleId(mod,lay,lad,det);
 183       lay--;
 184
 185       // Hits
 186       AliITSmodule *modu = ITS->GetModule(mod);
 187       TObjArray *arrHits = modu->GetHits();
 188       Int_t nhits = arrHits->GetEntriesFast();
 189       for (Int_t iHit=0;iHit<nhits;iHit++) {
 190         AliITShit *hit = (AliITShit*) arrHits->At(iHit);
 191         Int_t iMod=hit->GetModule();
 192         hlayer->Fill(lay);
 193         Double_t xl,yl,zl,xl0,yl0,zl0;
 194         Double_t xg,yg,zg,tof,xg0,yg0,zg0,tof0;
 195         Float_t hitloc[3],hitglo[3];
 196         hit->GetPositionL(xl,yl,zl,tof);
 197         hit->GetPositionL0(xl0,yl0,zl0,tof0);
 198         hit->GetPositionG(xg,yg,zg,tof);
 199         hit->GetPositionG0(xg0,yg0,zg0,tof0);
 200         Double_t hitlen=TMath::Abs(yl-yl0);
 201         if(hitlen<0.005) continue; // remove hits "shorter" than 50 um
 202         if(lay > 1 && hitlen<0.025) continue; // remove hits "shorter" than 250 um in SDD,SSD
 203         hitloc[0]=0.5*(xl+xl0);
 204         hitloc[1]=0.5*(yl+yl0);
 205         hitloc[2]=0.5*(zl+zl0);
 206         hitglo[0]=0.5*(xg+xg0);
 207         hitglo[1]=0.5*(yg+yg0);
 208         hitglo[2]=0.5*(zg+zg0);
 209         Float_t edep=hit->GetIonization()*1000000;
 210         Float_t rad=TMath::Sqrt(hitglo[0]*hitglo[0]+hitglo[1]*hitglo[1]);
 211         Float_t phi=TMath::ATan2(hitglo[1],hitglo[0]);
 212         hmod[lay]->Fill(iMod);
 213         hzl[lay]->Fill(hitloc[2]);
 214         hxl[lay]->Fill(hitloc[0]);
 215         hzg[lay]->Fill(hitglo[2]);
 216         hyg[lay]->Fill(hitglo[1]);
 217         hxg[lay]->Fill(hitglo[0]);
 218         hr[lay]->Fill(rad);
 219         hphi[lay]->Fill(phi);
 220         hq[lay]->Fill(edep);
 221         hzphi[lay]->Fill(phi,hitglo[2]);
 222       }
 223     }
 224   }
 225
 226   gStyle->SetOptStat(10);
 227   gStyle->SetPadBottomMargin(0.14);
 228
 229
 230   TCanvas **c=new TCanvas*[6];
 231   Char_t ctit[30];
 232   for(Int_t iLay=0;iLay<6;iLay++){
 233     sprintf(name,"can%d",iLay+1);
 234     sprintf(ctit,"Layer %d",iLay+1);
 235     c[iLay]=new TCanvas(name,ctit,1200,900);
 236     c[iLay]->Divide(3,3,0.001,0.001);
 237     c[iLay]->cd(1);
 238     hmod[iLay]->Draw();
 239     c[iLay]->cd(2);
 240     hxl[iLay]->Draw();
 241     c[iLay]->cd(3);
 242     hzl[iLay]->Draw();
 243     c[iLay]->cd(4);
 244     hxg[iLay]->Draw();
 245     c[iLay]->cd(5);
 246     hyg[iLay]->Draw();
 247     c[iLay]->cd(6);
 248     hzg[iLay]->Draw();
 249     c[iLay]->cd(7);
 250     hr[iLay]->Draw();
 251     c[iLay]->cd(8);
 252     hphi[iLay]->Draw();
 253     c[iLay]->cd(9);
 254     hq[iLay]->Draw();
 255   }
 256
 257   gStyle->SetPalette(1);
 258   TLatex* tstat=new TLatex();
 259   tstat->SetNDC();
 260   TCanvas* cspd=new TCanvas("cspd","SPD",1000,600);
 261   cspd->Divide(2,1);
 262   cspd->cd(1);
 263   hzphi[0]->Draw("colz");
 264   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[0]->GetEntries())));
 265   cspd->cd(2);
 266   hzphi[1]->Draw("colz");
 267   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[1]->GetEntries())));
 268
 269   TCanvas* csdd=new TCanvas("csdd","SDD",1000,600);
 270   csdd->Divide(2,1);
 271   csdd->cd(1);
 272   hzphi[2]->Draw("colz");
 273   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[2]->GetEntries())));
 274   csdd->cd(2);
 275   hzphi[3]->Draw("colz");
 276   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[3]->GetEntries())));
 277
 278   TCanvas* cssd=new TCanvas("cssd","SSD",1000,600);
 279   cssd->Divide(2,1);
 280   cssd->cd(1);
 281   hzphi[4]->Draw("colz");
 282   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[4]->GetEntries())));
 283   cssd->cd(2);
 284   hzphi[5]->Draw("colz");
 285   tstat->DrawLatex(0.6,0.92,Form("# Clusters = %d",int(hzphi[5]->GetEntries())));
 286
 287 }