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