[u/mrichter/AliRoot.git] / ITS / ShowITSRecPoints.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TCanvas.h>
#include <TClassTable.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TGeoManager.h>
#include <TH1.h>
#include <TInterpreter.h>
#include "AliCDBManager.h"
#include "AliGeomManager.h"
#include "AliHeader.h"
#include "AliITS.h"
#include "AliITSDetTypeRec.h"
#include "AliITSgeom.h"
#include "AliITSRecPoint.h"
#include "AliRun.h"
#endif

Int_t ShowITSRecPoints(){
  ///////////////////////////////////////////////////////////////////////
  // Macro to check clusters in the 6 ITS layers                       //
  // Provides:                                                         //
  //  6 canvases with 9 plots each (1 canvas for each layer)           //
  //  3 canvases with cluster XY coordinates for the first 3 events    //
  ///////////////////////////////////////////////////////////////////////

  if (gClassTable->GetID("AliRun") < 0) {
    gInterpreter->ExecuteMacro("loadlibs.C");
  }
  else { 
    if(gAlice){
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice=0;
    }
  }
  // Set OCDB if needed
  AliCDBManager* man = AliCDBManager::Instance();
  if (!man->IsDefaultStorageSet()) {
    printf("Setting a local default storage and run number 0\n");
    man->SetDefaultStorage("local://$ALICE_ROOT");
    man->SetRun(0);
  }
  else {
    printf("Using deafult storage \n");
  }
 
  // retrives geometry 
  if(!gGeoManager){
    AliGeomManager::LoadGeometry("geometry.root");
  }

  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  if (rl == 0x0){
    cerr<<"Can not open session RL=NULL"<< endl;
    return -1;
  }
  Int_t retval = rl->LoadgAlice();
  if (retval){
    cerr<<"LoadgAlice returned error"<<endl;
    return -1;
  }
  gAlice=rl->GetAliRun();

  retval = rl->LoadHeader();
  if (retval){
    cerr<<"LoadHeader returned error"<<endl;
    return -1;
  }


  AliITSLoader* ITSloader =  (AliITSLoader*) rl->GetLoader("ITSLoader");
  if(!ITSloader){
    cerr<<"ITS loader not found"<<endl;
    return -1;
  }
  ITSloader->LoadRecPoints("read");

  Float_t cluglo[3]={0.,0.,0.}; 
  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
  ITS->SetTreeAddress();
  AliITSgeom *geom = ITS->GetITSgeom();
  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
  detTypeRec->SetITSgeom(ITSloader->GetITSgeom());
  detTypeRec->SetDefaults();

  Int_t modmin=geom->GetStartDet(0);
  Int_t modmax=geom->GetLastDet(2);
  Int_t totmod=modmax-modmin;
  Float_t xlim[6]={4.5,7.5,16.,26.,40.,45.};
  Float_t zlim[6]={15.,15.,22.,30.,45.,55.};

  TH1F* hlayer=new TH1F("hlayer","",6,0.5,6.5);
  TH1F** hmod=new TH1F*[6];
  TH1F** hxl=new TH1F*[6];
  TH1F** hzl=new TH1F*[6];
  TH1F** hxg=new TH1F*[6];
  TH1F** hyg=new TH1F*[6];
  TH1F** hzg=new TH1F*[6];
  TH1F** hr=new TH1F*[6];
  TH1F** hphi=new TH1F*[6];
  TH1F** hq=new TH1F*[6];

  Char_t name[10];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"hmod%d",iLay+1);
    hmod[iLay]=new TH1F(name,"",totmod,modmin-0.5,modmax+0.5);
    hmod[iLay]->GetXaxis()->SetTitle("Module");
    hmod[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxloc%d",iLay+1);
    hxl[iLay]=new TH1F(name,"",100,-4.,4.);
    hxl[iLay]->GetXaxis()->SetTitle("Xloc (cm)");
    hxl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hzloc%d",iLay+1);
    hzl[iLay]=new TH1F(name,"",100,-4.,4.);
    hzl[iLay]->GetXaxis()->SetTitle("Zloc (cm)");
    hzl[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hxgl%d",iLay+1);
    hxg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hxg[iLay]->GetXaxis()->SetTitle("Xglob (cm)");
    hxg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hygl%d",iLay+1);
    hyg[iLay]=new TH1F(name,"",100,-xlim[iLay],xlim[iLay]);
    hyg[iLay]->GetXaxis()->SetTitle("Yglob (cm)");
    hyg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hzgl%d",iLay+1);
    hzg[iLay]=new TH1F(name,"",100,-zlim[iLay],zlim[iLay]);
    hzg[iLay]->GetXaxis()->SetTitle("Zglob (cm)");
    hzg[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hr%d",iLay+1);
    hr[iLay]=new TH1F(name,"",100,0.,50.);
    hr[iLay]->GetXaxis()->SetTitle("r (cm)");
    hr[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hphi%d",iLay+1);
    hphi[iLay]=new TH1F(name,"",100,-TMath::Pi(),TMath::Pi());    
    hphi[iLay]->GetXaxis()->SetTitle("#varphi (rad)");
    hphi[iLay]->GetXaxis()->CenterTitle();
    sprintf(name,"hq%d",iLay+1);
    hq[iLay]=new TH1F(name,"",100,0.,300.);    
    hq[iLay]->GetXaxis()->SetTitle("Charge (keV)");
    hq[iLay]->GetXaxis()->CenterTitle();
  }

  TGraph **gpts=new TGraph*[3];
  TGraph2D **gpts3d=new TGraph2D*[3];
  for(Int_t i=0;i<3;i++){
    gpts[i]=new TGraph(0);
    gpts3d[i]=new TGraph2D(0);
  }
  Int_t totev=rl->GetNumberOfEvents();
  printf("Total Number of events = %d\n",totev);

  for(Int_t iev=0;iev<totev;iev++){
    rl->GetEvent(iev);
    TTree *TR = ITSloader->TreeR();
    TClonesArray *ITSrec  = detTypeRec->RecPoints();
    TBranch *branch = 0;
    if(TR && ITSrec){
      branch = ITSloader->TreeR()->GetBranch("ITSRecPoints");
      if(branch)branch->SetAddress(&ITSrec);
    }
    Int_t nparticles = rl->GetHeader()->GetNtrack();
    cout<<"Event #"<<iev<<"   #Particles="<<nparticles<<endl;


    Int_t ipt=0;
    for(Int_t subd=0;subd<3;subd++){

      Int_t first = geom->GetStartDet(subd);
      Int_t last = geom->GetLastDet(subd);

      for (Int_t mod=first; mod<=last; mod++){
	detTypeRec->ResetRecPoints();
	branch->GetEvent(mod);
	Int_t nrecp = ITSrec->GetEntries();
	if(nrecp>0){
	  for(Int_t irec=0;irec<nrecp;irec++) {
	    AliITSRecPoint *recp = (AliITSRecPoint*)ITSrec->At(irec);
	    Int_t lay=recp->GetLayer();
	    hlayer->Fill(lay);
	    recp->GetGlobalXYZ(cluglo);
	    Float_t rad=TMath::Sqrt(cluglo[0]*cluglo[0]+cluglo[1]*cluglo[1]); 
	    Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
	    if(iev<3){
 	      gpts[iev]->SetPoint(ipt,cluglo[0],cluglo[1]);
 	      gpts3d[iev]->SetPoint(ipt,cluglo[0],cluglo[1],cluglo[2]);
	      ipt++;
	    }
	    hmod[lay]->Fill(mod);
	    hzl[lay]->Fill(recp->GetDetLocalZ());
	    hxl[lay]->Fill(recp->GetDetLocalX());
	    hzg[lay]->Fill(cluglo[2]);
	    hyg[lay]->Fill(cluglo[1]);
	    hxg[lay]->Fill(cluglo[0]);
	    hr[lay]->Fill(rad);
	    hphi[lay]->Fill(phi);
	    hq[lay]->Fill(recp->GetQ());
	  }
	}
      }
    }
  }
  TCanvas **c=new TCanvas*[6];
  Char_t ctit[30];
  for(Int_t iLay=0;iLay<6;iLay++){
    sprintf(name,"can%d",iLay+1);
    sprintf(ctit,"Layer %d",iLay+1);
    c[iLay]=new TCanvas(name,ctit,1200,900);
    c[iLay]->Divide(3,3,0.001,0.001);
    c[iLay]->cd(1);
    hmod[iLay]->Draw();
    c[iLay]->cd(2);
    hxl[iLay]->Draw();
    c[iLay]->cd(3);
    hzl[iLay]->Draw();
    c[iLay]->cd(4);
    hxg[iLay]->Draw();
    c[iLay]->cd(5);
    hyg[iLay]->Draw();
    c[iLay]->cd(6);
    hzg[iLay]->Draw();
    c[iLay]->cd(7);
    hr[iLay]->Draw();
    c[iLay]->cd(8);
    hphi[iLay]->Draw();    
    c[iLay]->cd(9);
    hq[iLay]->Draw();    
  }

  TCanvas *cev0;
  cev0=new TCanvas("cev0","Event 0",600,600);
  gpts[0]->SetMarkerStyle(7);
  gpts[0]->SetTitle(0);
  gpts[0]->Draw("AP");

  TCanvas *cev1;
  cev1=new TCanvas("cev1","Event 1",600,600);
  gpts[1]->SetMarkerStyle(7);
  gpts[1]->SetTitle(0);
  gpts[1]->Draw("AP");

  TCanvas *cev2;
  cev2=new TCanvas("cev2","Event 2",600,600);
  gpts[2]->SetMarkerStyle(7);
  gpts[2]->SetTitle(0);
  gpts[2]->Draw("AP");


  return 0;
}
Commit	Line	Data
29d63bb9	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