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 <TInterpreter.h>
   9 #include "AliCDBManager.h"
  10 #include "AliGeomManager.h"
  11 #include "AliHeader.h"
  12 #include "AliITS.h"
  13 #include "AliITSDetTypeRec.h"
  14 #include "AliITSgeom.h"
  15 #include "AliITSRecPoint.h"
  16 #include "AliRun.h"
  17 #endif
  18
  19 Int_t ShowITSRecPoints(){
  20   ///////////////////////////////////////////////////////////////////////
  21   // Macro to check clusters in the 6 ITS layers                       //
  22   // Provides:                                                         //
  23   //  6 canvases with 9 plots each (1 canvas for each layer)           //
  24   //  3 canvases with cluster XY coordinates for the first 3 events    //
  25   ///////////////////////////////////////////////////////////////////////
  26
  27   if (gClassTable->GetID("AliRun") < 0) {
  28     gInterpreter->ExecuteMacro("loadlibs.C");
  29   }
  30   else {
  31     if(gAlice){
  32       delete gAlice->GetRunLoader();
  33       delete gAlice;
  34       gAlice=0;
  35     }
  36   }
  37   // Set OCDB if needed
  38   AliCDBManager* man = AliCDBManager::Instance();
  39   if (!man->IsDefaultStorageSet()) {
  40     printf("Setting a local default storage and run number 0\n");
  41     man->SetDefaultStorage("local://$ALICE_ROOT");
  42     man->SetRun(0);
  43   }
  44   else {
  45     printf("Using deafult storage \n");
  46   }
  47
  48   // retrives geometry
  49   if(!gGeoManager){
  50     AliGeomManager::LoadGeometry("geometry.root");
  51   }
  52
  53   AliRunLoader* rl = AliRunLoader::Open("galice.root");
  54   if (rl == 0x0){
  55     cerr<<"Can not open session RL=NULL"<< endl;
  56     return -1;
  57   }
  58   Int_t retval = rl->LoadgAlice();
  59   if (retval){
  60     cerr<<"LoadgAlice returned error"<<endl;
  61     return -1;
  62   }
  63   gAlice=rl->GetAliRun();
  64
  65   retval = rl->LoadHeader();
  66   if (retval){
  67     cerr<<"LoadHeader returned error"<<endl;
  68     return -1;
  69   }
  70
  71
  72   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  73   if(!ITSloader){
  74     cerr<<"ITS loader not found"<<endl;
  75     return -1;
  76   }
  77   ITSloader->LoadRecPoints("read");
  78
  79   Float_t cluglo[3]={0.,0.,0.};
  80   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  81   ITS->SetTreeAddress();
  82   AliITSgeom *geom = ITS->GetITSgeom();
  83   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  84   detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
  85   detTypeRec->SetDefaults();
  86
  87   Int_t modmin=geom->GetStartDet(0);
  88   Int_t modmax=geom->GetLastDet(2);
  89   Int_t totmod=modmax-modmin;
  90   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  91   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
  92
  93   TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
  94   TH1F** hmod=new TH1F*[6];
  95   TH1F** hxl=new TH1F*[6];
  96   TH1F** hzl=new TH1F*[6];
  97   TH1F** hxg=new TH1F*[6];
  98   TH1F** hyg=new TH1F*[6];
  99   TH1F** hzg=new TH1F*[6];
 100   TH1F** hr=new TH1F*[6];
 101   TH1F** hphi=new TH1F*[6];
 102   TH1F** hq=new TH1F*[6];
 103
 104   Char_t name[10];
 105   for(Int_t iLay=0;iLay<6;iLay++){
 106     sprintf(name,"hmod%d",iLay+1);
 107     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
 108     hmod[iLay]->GetXaxis()->SetTitle("Module");
 109     hmod[iLay]->GetXaxis()->CenterTitle();
 110     sprintf(name,"hxloc%d",iLay+1);
 111     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
 112     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
 113     hxl[iLay]->GetXaxis()->CenterTitle();
 114     sprintf(name,"hzloc%d",iLay+1);
 115     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
 116     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
 117     hzl[iLay]->GetXaxis()->CenterTitle();
 118     sprintf(name,"hxgl%d",iLay+1);
 119     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
 120     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
 121     hxg[iLay]->GetXaxis()->CenterTitle();
 122     sprintf(name,"hygl%d",iLay+1);
 123     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
 124     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
 125     hyg[iLay]->GetXaxis()->CenterTitle();
 126     sprintf(name,"hzgl%d",iLay+1);
 127     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
 128     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
 129     hzg[iLay]->GetXaxis()->CenterTitle();
 130     sprintf(name,"hr%d",iLay+1);
 131     hr[iLay]=new TH1F(name,"",100,0.,50.);
 132     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
 133     hr[iLay]->GetXaxis()->CenterTitle();
 134     sprintf(name,"hphi%d",iLay+1);
 135     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
 136     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 137     hphi[iLay]->GetXaxis()->CenterTitle();
 138     sprintf(name,"hq%d",iLay+1);
 139     hq[iLay]=new TH1F(name,"",100,0.,300.);
 140     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
 141     hq[iLay]->GetXaxis()->CenterTitle();
 142   }
 143
 144   TGraph **gpts=new TGraph*[3];
 145   TGraph2D **gpts3d=new TGraph2D*[3];
 146   for(Int_t i=0;i<3;i++){
 147     gpts[i]=new TGraph(0);
 148     gpts3d[i]=new TGraph2D(0);
 149   }
 150   Int_t totev=rl->GetNumberOfEvents();
 151   printf("Total Number of events = %d\n",totev);
 152
 153   for(Int_t iev=0;iev<totev;iev++){
 154     rl->GetEvent(iev);
 155     TTree *TR = ITSloader->TreeR();
 156     TClonesArray *ITSrec  = detTypeRec->RecPoints();
 157     TBranch *branch = 0;
 158     if(TR && ITSrec){
 159       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 160       if(branch)branch->SetAddress(&ITSrec);
 161     }
 162     Int_t nparticles = rl->GetHeader()->GetNtrack();
 163     cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;
 164
 165
 166     Int_t ipt=0;
 167     for(Int_t subd=0;subd<3;subd++){
 168
 169       Int_t first = geom->GetStartDet(subd);
 170       Int_t last = geom->GetLastDet(subd);
 171
 172       for (Int_t mod=first; mod<=last; mod++){
 173         detTypeRec->ResetRecPoints();
 174         branch->GetEvent(mod);
 175         Int_t nrecp = ITSrec->GetEntries();
 176         if(nrecp>0){
 177           for(Int_t irec=0;irec<nrecp;irec++) {
 178             AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 179             Int_t lay=recp->GetLayer();
 180             hlayer->Fill(lay);
 181             recp->GetGlobalXYZ(cluglo);
 182             Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]);
 183             Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
 184             if(iev<3){
 185               gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
 186               gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
 187               ipt++;
 188             }
 189             hmod[lay]->Fill(mod);
 190             hzl[lay]->Fill(recp->GetDetLocalZ());
 191             hxl[lay]->Fill(recp->GetDetLocalX());
 192             hzg[lay]->Fill(cluglo[2]);
 193             hyg[lay]->Fill(cluglo[1]);
 194             hxg[lay]->Fill(cluglo[0]);
 195             hr[lay]->Fill(rad);
 196             hphi[lay]->Fill(phi);
 197             hq[lay]->Fill(recp->GetQ());
 198           }
 199         }
 200       }
 201     }
 202   }
 203   TCanvas **c=new TCanvas*[6];
 204   Char_t ctit[30];
 205   for(Int_t iLay=0;iLay<6;iLay++){
 206     sprintf(name,"can%d",iLay+1);
 207     sprintf(ctit,"Layer %d",iLay+1);
 208     c[iLay]=new TCanvas(name,ctit,1200,900);
 209     c[iLay]->Divide(3,3,0.001,0.001);
 210     c[iLay]->cd(1);
 211     hmod[iLay]->Draw();
 212     c[iLay]->cd(2);
 213     hxl[iLay]->Draw();
 214     c[iLay]->cd(3);
 215     hzl[iLay]->Draw();
 216     c[iLay]->cd(4);
 217     hxg[iLay]->Draw();
 218     c[iLay]->cd(5);
 219     hyg[iLay]->Draw();
 220     c[iLay]->cd(6);
 221     hzg[iLay]->Draw();
 222     c[iLay]->cd(7);
 223     hr[iLay]->Draw();
 224     c[iLay]->cd(8);
 225     hphi[iLay]->Draw();
 226     c[iLay]->cd(9);
 227     hq[iLay]->Draw();
 228   }
 229
 230   TCanvas *cev0;
 231   cev0=new TCanvas("cev0","Event 0",600,600);
 232   gpts[0]->SetMarkerStyle(7);
 233   gpts[0]->SetTitle(0);
 234   gpts[0]->Draw("AP");
 235
 236   TCanvas *cev1;
 237   cev1=new TCanvas("cev1","Event 1",600,600);
 238   gpts[1]->SetMarkerStyle(7);
 239   gpts[1]->SetTitle(0);
 240   gpts[1]->Draw("AP");
 241
 242   TCanvas *cev2;
 243   cev2=new TCanvas("cev2","Event 2",600,600);
 244   gpts[2]->SetMarkerStyle(7);
 245   gpts[2]->SetTitle(0);
 246   gpts[2]->Draw("AP");
 247
 248
 249   return 0;
 250 }
 251
 252