[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 "AliGeomManager.h"
#include "AliHeader.h"
#include "AliITS.h"
#include "AliITSDetTypeRec.h"
#include "AliITSgeomTGeo.h"
#include "AliITSRecPoint.h"
#include "AliRun.h"
#endif

Int_t ShowITSRecPoints(Int_t nevfordisp=0){
  ///////////////////////////////////////////////////////////////////////
  // 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");
  }
 
  // 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;
  }

  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.}; 
  AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();

  Int_t totmod=AliITSgeomTGeo::GetNModules();
  Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
  Int_t modmax=AliITSgeomTGeo::GetNModules()-1;

  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,5.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 *gptsXY=new TGraph(0);
  TGraph *gptsRZ=new TGraph(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 mod=modmin; mod<=modmax; 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==nevfordisp){
	    gptsXY->SetPoint(ipt,cluglo[0],cluglo[1]);
	    if(cluglo[1]>0) gptsRZ->SetPoint(ipt,cluglo[2],rad);
	    else gptsRZ->SetPoint(ipt,cluglo[2],-rad);
	    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;
  sprintf(ctit,"Event %d XY",nevfordisp);
  cev0=new TCanvas("cev0",ctit,600,600);
  if(gptsXY->GetN()>0){
    gptsXY->SetMarkerStyle(7);
    gptsXY->SetTitle(0);
    gptsXY->Draw("AP");
    gptsXY->GetXaxis()->SetTitle("Xglob");
    gptsXY->GetYaxis()->SetTitle("Yglob");
   }

  TCanvas *cev1;
  sprintf(ctit,"Event %d Zr",nevfordisp);
  cev1=new TCanvas("cev1",ctit,600,600);
  if(gptsRZ->GetN()>0){
    gptsRZ->SetMarkerStyle(7);
    gptsRZ->SetTitle(0);
    gptsRZ->Draw("AP");
    gptsRZ->GetXaxis()->SetTitle("Zglob");
    gptsRZ->GetYaxis()->SetTitle("Radius");
  }

  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>
29d63bb9	9	#include "AliGeomManager.h"
	10	#include "AliHeader.h"
	11	#include "AliITS.h"
	12	#include "AliITSDetTypeRec.h"
7285ac7a	13	#include "AliITSgeomTGeo.h"
29d63bb9	14	#include "AliITSRecPoint.h"
	15	#include "AliRun.h"
	16	#endif
	17
7285ac7a	18	Int_t ShowITSRecPoints(Int_t nevfordisp=0){
29d63bb9	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	}
29d63bb9	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	}
29d63bb9	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.};
29d63bb9	49	AliITSDetTypeRec* detTypeRec = new AliITSDetTypeRec();
29d63bb9	50
7285ac7a	51	Int_t totmod=AliITSgeomTGeo::GetNModules();
	52	Int_t modmin=AliITSgeomTGeo::GetModuleIndex(1,1,1);
	53	Int_t modmax=AliITSgeomTGeo::GetNModules()-1;
	54
29d63bb9	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
97e78ca8	58	TH1F* hlayer=new TH1F("hlayer","",6,-0.5,5.5);
29d63bb9	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
7285ac7a	109	TGraph *gptsXY=new TGraph(0);
	110	TGraph *gptsRZ=new TGraph(0);
	111
29d63bb9	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	Int_t nparticles = rl->GetHeader()->GetNtrack();
	125	cout<<"Event #"<<iev<<" #Particles="<<nparticles<<endl;
	126
	127
	128	Int_t ipt=0;
7285ac7a	129	for (Int_t mod=modmin; mod<=modmax; mod++){
	130	detTypeRec->ResetRecPoints();
	131	branch->GetEvent(mod);
	132	Int_t nrecp = ITSrec->GetEntries();
	133	if(nrecp>0){
	134	for(Int_t irec=0;irec<nrecp;irec++) {
	135	AliITSRecPoint recp = (AliITSRecPoint)ITSrec->At(irec);
	136	Int_t lay=recp->GetLayer();
	137	hlayer->Fill(lay);
	138	recp->GetGlobalXYZ(cluglo);
	139	Float_t rad=TMath::Sqrt(cluglo[0]cluglo[0]+cluglo[1]cluglo[1]);
	140	Float_t phi=TMath::ATan2(cluglo[1],cluglo[0]);
	141	if(iev==nevfordisp){
	142	gptsXY->SetPoint(ipt,cluglo[0],cluglo[1]);
	143	if(cluglo[1]>0) gptsRZ->SetPoint(ipt,cluglo[2],rad);
	144	else gptsRZ->SetPoint(ipt,cluglo[2],-rad);
	145	ipt++;
29d63bb9	146	}
7285ac7a	147	hmod[lay]->Fill(mod);
	148	hzl[lay]->Fill(recp->GetDetLocalZ());
	149	hxl[lay]->Fill(recp->GetDetLocalX());
	150	hzg[lay]->Fill(cluglo[2]);
	151	hyg[lay]->Fill(cluglo[1]);
	152	hxg[lay]->Fill(cluglo[0]);
	153	hr[lay]->Fill(rad);
	154	hphi[lay]->Fill(phi);
	155	hq[lay]->Fill(recp->GetQ());
29d63bb9	156	}
	157	}
	158	}
	159	}
7285ac7a	160
29d63bb9	161	TCanvas *c=new TCanvas[6];
	162	Char_t ctit[30];
	163	for(Int_t iLay=0;iLay<6;iLay++){
	164	sprintf(name,"can%d",iLay+1);
	165	sprintf(ctit,"Layer %d",iLay+1);
	166	c[iLay]=new TCanvas(name,ctit,1200,900);
	167	c[iLay]->Divide(3,3,0.001,0.001);
	168	c[iLay]->cd(1);
	169	hmod[iLay]->Draw();
	170	c[iLay]->cd(2);
	171	hxl[iLay]->Draw();
	172	c[iLay]->cd(3);
	173	hzl[iLay]->Draw();
	174	c[iLay]->cd(4);
	175	hxg[iLay]->Draw();
	176	c[iLay]->cd(5);
	177	hyg[iLay]->Draw();
	178	c[iLay]->cd(6);
	179	hzg[iLay]->Draw();
	180	c[iLay]->cd(7);
	181	hr[iLay]->Draw();
	182	c[iLay]->cd(8);
	183	hphi[iLay]->Draw();
	184	c[iLay]->cd(9);
	185	hq[iLay]->Draw();
	186	}
	187
	188	TCanvas *cev0;
7285ac7a	189	sprintf(ctit,"Event %d XY",nevfordisp);
	190	cev0=new TCanvas("cev0",ctit,600,600);
	191	if(gptsXY->GetN()>0){
	192	gptsXY->SetMarkerStyle(7);
	193	gptsXY->SetTitle(0);
	194	gptsXY->Draw("AP");
	195	gptsXY->GetXaxis()->SetTitle("Xglob");
	196	gptsXY->GetYaxis()->SetTitle("Yglob");
	197	}
29d63bb9	198
29d63bb9	199	TCanvas *cev1;
7285ac7a	200	sprintf(ctit,"Event %d Zr",nevfordisp);
	201	cev1=new TCanvas("cev1",ctit,600,600);
	202	if(gptsRZ->GetN()>0){
	203	gptsRZ->SetMarkerStyle(7);
	204	gptsRZ->SetTitle(0);
	205	gptsRZ->Draw("AP");
	206	gptsRZ->GetXaxis()->SetTitle("Zglob");
	207	gptsRZ->GetYaxis()->SetTitle("Radius");
	208	}
29d63bb9	209
	210	return 0;
	211	}
	212
	213