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 "AliITSgeomTGeo.h"
  14 #include "AliITSRecPoint.h"
  15 #include "AliRun.h"
  16 #endif
  17
  18 Int_t ShowITSRecPoints(Int_t nevfordisp=0){
  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
  41   AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  42   if(!ITSloader){
  43     cerr<<"ITS loader not found"<<endl;
  44     return -1;
  45   }
  46   ITSloader->LoadRecPoints("read");
  47
  48   Float_t cluglo[3]={0.,0.,0.};
  49   AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  50
  51   Int_t totmod=AliITSgeomTGeo::GetNModules();
  52   Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
  53   Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
  54
  55   Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  56   Float_t zlim[6]={15.,15.,22.,30.,45.,55.};
  57
  58   TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
  59   TH1F** hmod=new TH1F*[6];
  60   TH1F** hxl=new TH1F*[6];
  61   TH1F** hzl=new TH1F*[6];
  62   TH1F** hxg=new TH1F*[6];
  63   TH1F** hyg=new TH1F*[6];
  64   TH1F** hzg=new TH1F*[6];
  65   TH1F** hr=new TH1F*[6];
  66   TH1F** hphi=new TH1F*[6];
  67   TH1F** hq=new TH1F*[6];
  68
  69   Char_t name[10];
  70   for(Int_t iLay=0;iLay<6;iLay++){
  71     sprintf(name,"hmod%d",iLay+1);
  72     hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
  73     hmod[iLay]->GetXaxis()->SetTitle("Module");
  74     hmod[iLay]->GetXaxis()->CenterTitle();
  75     sprintf(name,"hxloc%d",iLay+1);
  76     hxl[iLay]=new TH1F(name,"",100,-4.,4.);
  77     hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
  78     hxl[iLay]->GetXaxis()->CenterTitle();
  79     sprintf(name,"hzloc%d",iLay+1);
  80     hzl[iLay]=new TH1F(name,"",100,-4.,4.);
  81     hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
  82     hzl[iLay]->GetXaxis()->CenterTitle();
  83     sprintf(name,"hxgl%d",iLay+1);
  84     hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
  85     hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
  86     hxg[iLay]->GetXaxis()->CenterTitle();
  87     sprintf(name,"hygl%d",iLay+1);
  88     hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
  89     hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
  90     hyg[iLay]->GetXaxis()->CenterTitle();
  91     sprintf(name,"hzgl%d",iLay+1);
  92     hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
  93     hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
  94     hzg[iLay]->GetXaxis()->CenterTitle();
  95     sprintf(name,"hr%d",iLay+1);
  96     hr[iLay]=new TH1F(name,"",100,0.,50.);
  97     hr[iLay]->GetXaxis()->SetTitle("r (cm)");
  98     hr[iLay]->GetXaxis()->CenterTitle();
  99     sprintf(name,"hphi%d",iLay+1);
 100     hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());
 101     hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
 102     hphi[iLay]->GetXaxis()->CenterTitle();
 103     sprintf(name,"hq%d",iLay+1);
 104     hq[iLay]=new TH1F(name,"",100,0.,300.);
 105     hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
 106     hq[iLay]->GetXaxis()->CenterTitle();
 107   }
 108
 109   TGraph *gptsXY=new TGraph(0);
 110   TGraph *gptsRZ=new TGraph(0);
 111
 112   Int_t totev=rl->GetNumberOfEvents();
 113   printf("Total Number of events = %d\n",totev);
 114
 115   for(Int_t iev=0;iev<totev;iev++){
 116     rl->GetEvent(iev);
 117     TTree *TR = ITSloader->TreeR();
 118     TClonesArray *ITSrec  = detTypeRec->RecPoints();
 119     TBranch *branch = 0;
 120     if(TR && ITSrec){
 121       branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
 122       if(branch)branch->SetAddress(&ITSrec);
 123     }
 124     if(iev%100==0) printf("Event #%d\n",iev);
 125
 126
 127     Int_t ipt=0;
 128     for (Int_t mod=modmin; mod<=modmax; mod++){
 129       detTypeRec->ResetRecPoints();
 130       branch->GetEvent(mod);
 131       Int_t nrecp = ITSrec->GetEntries();
 132       if(nrecp>0){
 133         for(Int_t irec=0;irec<nrecp;irec++) {
 134           AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
 135           Int_t lay=recp->GetLayer();
 136           hlayer->Fill(lay);
 137           recp->GetGlobalXYZ(cluglo);
 138           Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]);
 139           Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
 140           if(iev==nevfordisp){
 141             gptsXY->SetPoint(ipt,cluglo[0],cluglo[1]);
 142             if(cluglo[1]>0) gptsRZ->SetPoint(ipt,cluglo[2],rad);
 143             else gptsRZ->SetPoint(ipt,cluglo[2],-rad);
 144             ipt++;
 145           }
 146           hmod[lay]->Fill(mod);
 147           hzl[lay]->Fill(recp->GetDetLocalZ());
 148           hxl[lay]->Fill(recp->GetDetLocalX());
 149           hzg[lay]->Fill(cluglo[2]);
 150           hyg[lay]->Fill(cluglo[1]);
 151           hxg[lay]->Fill(cluglo[0]);
 152           hr[lay]->Fill(rad);
 153           hphi[lay]->Fill(phi);
 154           hq[lay]->Fill(recp->GetQ());
 155         }
 156       }
 157     }
 158   }
 159
 160   TCanvas **c=new TCanvas*[6];
 161   Char_t ctit[30];
 162   for(Int_t iLay=0;iLay<6;iLay++){
 163     sprintf(name,"can%d",iLay+1);
 164     sprintf(ctit,"Layer %d",iLay+1);
 165     c[iLay]=new TCanvas(name,ctit,1200,900);
 166     c[iLay]->Divide(3,3,0.001,0.001);
 167     c[iLay]->cd(1);
 168     hmod[iLay]->Draw();
 169     c[iLay]->cd(2);
 170     hxl[iLay]->Draw();
 171     c[iLay]->cd(3);
 172     hzl[iLay]->Draw();
 173     c[iLay]->cd(4);
 174     hxg[iLay]->Draw();
 175     c[iLay]->cd(5);
 176     hyg[iLay]->Draw();
 177     c[iLay]->cd(6);
 178     hzg[iLay]->Draw();
 179     c[iLay]->cd(7);
 180     hr[iLay]->Draw();
 181     c[iLay]->cd(8);
 182     hphi[iLay]->Draw();
 183     c[iLay]->cd(9);
 184     hq[iLay]->Draw();
 185   }
 186
 187   TCanvas *cev0;
 188   sprintf(ctit,"Event %d XY",nevfordisp);
 189   cev0=new TCanvas("cev0",ctit,600,600);
 190   if(gptsXY->GetN()>0){
 191     gptsXY->SetMarkerStyle(7);
 192     gptsXY->SetTitle(0);
 193     gptsXY->Draw("AP");
 194     gptsXY->GetXaxis()->SetTitle("Xglob");
 195     gptsXY->GetYaxis()->SetTitle("Yglob");
 196    }
 197
 198   TCanvas *cev1;
 199   sprintf(ctit,"Event %d Zr",nevfordisp);
 200   cev1=new TCanvas("cev1",ctit,600,600);
 201   if(gptsRZ->GetN()>0){
 202     gptsRZ->SetMarkerStyle(7);
 203     gptsRZ->SetTitle(0);
 204     gptsRZ->Draw("AP");
 205     gptsRZ->GetXaxis()->SetTitle("Zglob");
 206     gptsRZ->GetYaxis()->SetTitle("Radius");
 207   }
 208
 209   return 0;
 210 }
 211
 212