[u/mrichter/AliRoot.git] / TRD / AliTRDclusterErrors.C

#ifndef __CINT__
  #include <iostream.h>

  #include "AliTRDtracker.h"
  #include "AliTRDcluster.h"
  #include "AliTRDhit.h"
  #include "AliTRDv1.h"
  #include "AliTRDgeometry.h"
  #include "AliTRDparameter.h"

  #include "alles.h"
  #include "TFile.h"
  #include "TStopwatch.h"

#endif    

void AliTRDclusterErrors() {

  Int_t No_of_tracks_to_analyze = 40;

  const Int_t tbpp = 15;
  const Int_t nPlanes = 6;
  const Int_t ntb = tbpp * nPlanes; 

  TH1F *hy = new TH1F("delta R*phi","Cluster displacement in R*phi",200,-1.,1.); 
  TH1F *hyp = new TH1F("delta R*phi pos","delta R*phi, positive",200,-1.,1.); 
  TH1F *hym = new TH1F("delta R*phi neg","delta R*phi, negative",200,-1.,1.); 

  TH1F *hyn = new TH1F("Norm., d(R*phi)","Norm. cluster displacement in R*phi",400,-8.,8.); 
  TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.); 
  TH2F *hy2 = new TH2F("Amp vs delta R*phi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.); 
  TH2F *herr = new TH2F("sigmaY vs delta R*phi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1); 
  TH2F *hy3 = new TH2F("Position within pad vs delta R*phi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5); 
  TH2F *hy4 = new TH2F("local tb vs delta R*phi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5); 

  hy->SetXTitle("Displacement, cm"); 
  hyn->SetXTitle("Displacement, SigmaY"); 
  hy2->SetXTitle("Displacement, cm"); 
  hy2->SetYTitle("Amplitude"); 
  hz->SetXTitle("Displacement, cm"); 
  hy3->SetXTitle("Displacement, cm"); 
  hy3->SetYTitle("Position, cm"); 
  hy4->SetXTitle("Displacement, cm"); 
  hy4->SetYTitle("local time bin"); 

  /*
  // Dynamically link some shared libs
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
    cout << "Loaded shared libraries" << endl;
  } 
  */      

  // Load clusters
  Char_t *clusterfile = "AliTRDclusters.root";
  Int_t   nEvent  = 0;

  TObjArray carray(2000);
  TObjArray *ClustersArray = &carray;
  TFile *geofile =TFile::Open("AliTRDclusters.root");   

  AliTRDparameter *par = (AliTRDparameter*) geofile->Get("TRDparameter");
  AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry");   

  AliTRDtracker *Tracker = new AliTRDtracker(geofile);
  Tracker->SetEventNumber(nEvent);
  Tracker->ReadClusters(ClustersArray,clusterfile);
  Int_t nClusters = carray.GetEntriesFast();

  printf("Total number of clusters %d \n", nClusters);

  Char_t *alifile = "galice.root";
  // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
  TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
  if (!gafl) {
    cout << "Open the ALIROOT-file " << alifile << endl;
    gafl = new TFile(alifile);
  }
  else {
    cout << alifile << " is already open" << endl;
  }

  // Get AliRun object from file or create it if not on file
  gAlice = (AliRun*) gafl->Get("gAlice");
  if (gAlice)
    cout << "AliRun object found on file" << endl;
  else
    gAlice = new AliRun("gAlice","Alice test program");

  AliTRDv1       *fTRD           = (AliTRDv1*) gAlice->GetDetector("TRD");

  // Import the Trees for the event nEvent in the file
  Int_t nparticles = gAlice->GetEvent(nEvent);

  TObjArray *particles=gAlice->Particles();

  TTree *hitTree = gAlice->TreeH();
  Int_t  nBytes = 0;

  // Get the number of entries in the hit tree
  // (Number of primary particles creating a hit somewhere)
  Int_t nTrack = (Int_t) hitTree->GetEntries();
  cout << " Found " << nTrack << " primary particles with hits" << endl;   
  No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack);


  // Loop through particles and fill histoes

  Float_t hitY[ntb];
  Float_t hitZ[ntb];
  Float_t hitO[ntb];
  Float_t clusterY[ntb];
  Float_t clusterZ[ntb];
  Float_t clusterQ[ntb];
  Float_t clusterSigmaY[ntb];
  Float_t pos[3];
  Float_t rot[3];
  Float_t global[3];
  Int_t det = 0;
  Int_t track_index[3];

  printf("\n");

  for (Int_t ii=0; ii<No_of_tracks_to_analyze; ii++) {
    printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze); 

    for(Int_t plane = 0; plane < nPlanes; plane++) {
      for(Int_t ltb = 14; ltb > -1; ltb--) {

	if(ii >= nTrack) continue;

	TParticle *p = gAlice->Particle(ii);
	if (p->GetFirstMother()>=0) continue;
	TParticlePDG *pdg = p->GetPDG();
	Float_t charge=pdg->Charge(); 
	if(TMath::Abs(charge) < 0.5) continue;

	Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
	Double_t x = Tracker->GetX(0,plane,ltb);

	// loop through clusters

	Bool_t cluster_found = kFALSE;
	Int_t nw = 0;

	for (Int_t i = 0; i < nClusters; i++) {

	  AliTRDcluster *cl = (AliTRDcluster *) carray.UncheckedAt(i);

	  nw = cl->GetDetector(); 

	  nw = fGeo->GetPlane(nw);

	  if(nw != plane) continue; 
	  for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j); 
	  if((track_index[0] != ii) &&
	     (track_index[1] != ii) &&
	     (track_index[2] != ii)) continue;

	  nw = cl->GetLocalTimeBin(); 
	  if(nw != ltb) continue;

	  clusterY[gtb] = cl->GetY();
	  clusterZ[gtb] = cl->GetZ();
	  clusterQ[gtb] = cl->GetQ();

	  clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2());
	  cluster_found = kTRUE;
	  break;	  
	}

	if(!cluster_found) continue;

	gAlice->ResetHits();

	//	nBytes += hitTree->GetEvent(nPrimaries - ii - 1);
	nBytes += hitTree->GetEvent(nTrack - ii - 1);

	// Loop through the TRD hits
	Bool_t found_hit = kFALSE;


	for(AliTRDhit *hit = (AliTRDhit *) fTRD->FirstHit(-1); 
	    hit; 
	    hit = (AliTRDhit *) fTRD->NextHit()) {
	  nw = hit->Track();
	  if(nw != ii) continue;
	  det   = hit->GetDetector();
	  nw = fGeo->GetPlane(det);
	  if(nw != plane) continue;           


	  pos[0]=hit->X(); 
	  pos[1]=hit->Y();
	  pos[2]=hit->Z();
	  fGeo->Rotate(det,pos,rot);

	  if(TMath::Abs(rot[0]-x) > 0.01) continue;
	  hitY[gtb] = rot[1];
	  hitZ[gtb] = rot[2];

	  Float_t col0        = par->GetCol0(plane);
	  Float_t colPadSize  = par->GetColPadSize(plane);         
	  Float_t colH = (Int_t ((rot[1] -  col0)/colPadSize)) * colPadSize;  
	  hitO[gtb] = (rot[1] -  col0) - colH;
	  found_hit = kTRUE;
	  break;
	}

	if(!found_hit) continue;

	/*	
	printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n",
	       gtb, x, rot[0], rot[1], clusterY[gtb]);
	printf("\n                            Zhit - Zcl = %f - %f = %f\n",
	       rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]);


	printf("found hit within dx = %f - %f \n",rot[0],x);
	printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]);
	printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]);
	printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]);
	*/


	hy->Fill(hitY[gtb]-clusterY[gtb]);
	if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]);
	else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]);

	if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0)) 
	  hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]);
	hz->Fill(hitZ[gtb]-clusterZ[gtb]);
	hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]);
	hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]);  
	hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp));
	herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]);
      }
    }
  }
  
  gStyle->SetOptStat(1);
  gStyle->SetOptFit(1); 

  TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840);
  c->SetFillColor(10);
  c->Divide(2,2);
  c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw();
  c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw();
  c->cd(3); hyn->Draw();
  c->cd(4); hy4->Draw();

  TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
  c1->SetFillColor(10);
  c1->Divide(2,2);
  c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus");
  c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus");
  c1->cd(3); hy3->Draw();
  c1->cd(4); hy2->Draw();
   
}
Commit	Line	Data
59dfcadd	1	#ifndef __CINT__
	2	#include <iostream.h>
	3
	4	#include "AliTRDtracker.h"
	5	#include "AliTRDcluster.h"
	6	#include "AliTRDhit.h"
	7	#include "AliTRDv1.h"
	8	#include "AliTRDgeometry.h"
	9	#include "AliTRDparameter.h"
	10
	11	#include "alles.h"
	12	#include "TFile.h"
	13	#include "TStopwatch.h"
	14
	15	#endif
	16
	17	void AliTRDclusterErrors() {
	18
	19	Int_t No_of_tracks_to_analyze = 40;
	20
	21	const Int_t tbpp = 15;
	22	const Int_t nPlanes = 6;
	23	const Int_t ntb = tbpp * nPlanes;
	24
	25	TH1F hy = new TH1F("delta Rphi","Cluster displacement in R*phi",200,-1.,1.);
	26	TH1F hyp = new TH1F("delta Rphi pos","delta R*phi, positive",200,-1.,1.);
	27	TH1F hym = new TH1F("delta Rphi neg","delta R*phi, negative",200,-1.,1.);
	28
	29	TH1F hyn = new TH1F("Norm., d(Rphi)","Norm. cluster displacement in R*phi",400,-8.,8.);
	30	TH1F *hz = new TH1F("delta Z","Cluster displacement in Z",300,-10.,50.);
	31	TH2F hy2 = new TH2F("Amp vs delta Rphi","Amplitude versus delta R*phi",200,-5.,5.,200,0.,600.);
	32	TH2F herr = new TH2F("sigmaY vs delta Rphi","sigmaY vs delta R*phi",200,-1,1,200,0.,0.1);
	33	TH2F hy3 = new TH2F("Position within pad vs delta Rphi","Position within pad vs delta R*phi",200,-1.,1.,200,-0.5,1.5);
	34	TH2F hy4 = new TH2F("local tb vs delta Rphi","local tb vs delta R*phi",200,-1.,1.,20,-2.5,17.5);
	35
	36	hy->SetXTitle("Displacement, cm");
	37	hyn->SetXTitle("Displacement, SigmaY");
	38	hy2->SetXTitle("Displacement, cm");
	39	hy2->SetYTitle("Amplitude");
	40	hz->SetXTitle("Displacement, cm");
	41	hy3->SetXTitle("Displacement, cm");
	42	hy3->SetYTitle("Position, cm");
	43	hy4->SetXTitle("Displacement, cm");
	44	hy4->SetYTitle("local time bin");
	45
	46	/*
	47	// Dynamically link some shared libs
	48	if (gClassTable->GetID("AliRun") < 0) {
	49	gROOT->LoadMacro("loadlibs.C");
	50	loadlibs();
	51	cout << "Loaded shared libraries" << endl;
	52	}
	53	*/
	54
	55	// Load clusters
	56	Char_t *clusterfile = "AliTRDclusters.root";
	57	Int_t nEvent = 0;
	58
	59	TObjArray carray(2000);
	60	TObjArray *ClustersArray = &carray;
	61	TFile *geofile =TFile::Open("AliTRDclusters.root");
	62
	63	AliTRDparameter par = (AliTRDparameter) geofile->Get("TRDparameter");
	64	AliTRDgeometry fGeo = (AliTRDgeometry) geofile->Get("TRDgeometry");
65
66	AliTRDtracker *Tracker = new AliTRDtracker(geofile);
67	Tracker->SetEventNumber(nEvent);
68	Tracker->ReadClusters(ClustersArray,clusterfile);
69	Int_t nClusters = carray.GetEntriesFast();
70
71	printf("Total number of clusters %d \n", nClusters);
72
73	Char_t *alifile = "galice.root";
74	// Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
75	TFile gafl = (TFile) gROOT->GetListOfFiles()->FindObject(alifile);
76	if (!gafl) {
77	cout << "Open the ALIROOT-file " << alifile << endl;
78	gafl = new TFile(alifile);
79	}
80	else {
81	cout << alifile << " is already open" << endl;
82	}
83
84	// Get AliRun object from file or create it if not on file
85	gAlice = (AliRun*) gafl->Get("gAlice");
86	if (gAlice)
87	cout << "AliRun object found on file" << endl;
88	else
89	gAlice = new AliRun("gAlice","Alice test program");
90
91	AliTRDv1 fTRD = (AliTRDv1) gAlice->GetDetector("TRD");
92
93	// Import the Trees for the event nEvent in the file
94	Int_t nparticles = gAlice->GetEvent(nEvent);
95
96	TObjArray *particles=gAlice->Particles();
97
98	TTree *hitTree = gAlice->TreeH();
99	Int_t nBytes = 0;
100
101	// Get the number of entries in the hit tree
102	// (Number of primary particles creating a hit somewhere)
103	Int_t nTrack = (Int_t) hitTree->GetEntries();
104	cout << " Found " << nTrack << " primary particles with hits" << endl;
105	No_of_tracks_to_analyze = TMath::Min(No_of_tracks_to_analyze, nTrack);
106
107
108	// Loop through particles and fill histoes
109
110	Float_t hitY[ntb];
111	Float_t hitZ[ntb];
112	Float_t hitO[ntb];
113	Float_t clusterY[ntb];
114	Float_t clusterZ[ntb];
115	Float_t clusterQ[ntb];
116	Float_t clusterSigmaY[ntb];
117	Float_t pos[3];
118	Float_t rot[3];
119	Float_t global[3];
120	Int_t det = 0;
121	Int_t track_index[3];
122
123	printf("\n");
124
125	for (Int_t ii=0; ii<No_of_tracks_to_analyze; ii++) {
126	printf("track %d out of %d \n", ii+1 , No_of_tracks_to_analyze);
127
128	for(Int_t plane = 0; plane < nPlanes; plane++) {
129	for(Int_t ltb = 14; ltb > -1; ltb--) {
130
131	if(ii >= nTrack) continue;
132
133	TParticle *p = gAlice->Particle(ii);
134	if (p->GetFirstMother()>=0) continue;
135	TParticlePDG *pdg = p->GetPDG();
136	Float_t charge=pdg->Charge();
137	if(TMath::Abs(charge) < 0.5) continue;
138
139	Int_t gtb = Tracker->GetGlobalTimeBin(0,plane,ltb);
140	Double_t x = Tracker->GetX(0,plane,ltb);
141
142	// loop through clusters
143
144	Bool_t cluster_found = kFALSE;
145	Int_t nw = 0;
146
147	for (Int_t i = 0; i < nClusters; i++) {
148
149	AliTRDcluster cl = (AliTRDcluster ) carray.UncheckedAt(i);
150
151	nw = cl->GetDetector();
152
153	nw = fGeo->GetPlane(nw);
154
155	if(nw != plane) continue;
156	for(Int_t j=0; j<3; j++) track_index[j] = cl->GetLabel(j);
157	if((track_index[0] != ii) &&
158	(track_index[1] != ii) &&
159	(track_index[2] != ii)) continue;
160
161	nw = cl->GetLocalTimeBin();
162	if(nw != ltb) continue;
163
164	clusterY[gtb] = cl->GetY();
165	clusterZ[gtb] = cl->GetZ();
166	clusterQ[gtb] = cl->GetQ();
167
168	clusterSigmaY[gtb] = TMath::Sqrt(cl->GetSigmaY2());
169	cluster_found = kTRUE;
170	break;
171	}
172
173	if(!cluster_found) continue;
174
175	gAlice->ResetHits();
176
177	// nBytes += hitTree->GetEvent(nPrimaries - ii - 1);
178	nBytes += hitTree->GetEvent(nTrack - ii - 1);
179
180	// Loop through the TRD hits
181	Bool_t found_hit = kFALSE;
182
183
184	for(AliTRDhit hit = (AliTRDhit ) fTRD->FirstHit(-1);
185	hit;
186	hit = (AliTRDhit *) fTRD->NextHit()) {
187	nw = hit->Track();
188	if(nw != ii) continue;
189	det = hit->GetDetector();
190	nw = fGeo->GetPlane(det);
191	if(nw != plane) continue;
192
193
194	pos[0]=hit->X();
195	pos[1]=hit->Y();
196	pos[2]=hit->Z();
197	fGeo->Rotate(det,pos,rot);
198
199	if(TMath::Abs(rot[0]-x) > 0.01) continue;
200	hitY[gtb] = rot[1];
201	hitZ[gtb] = rot[2];
202
203	Float_t col0 = par->GetCol0(plane);
204	Float_t colPadSize = par->GetColPadSize(plane);
205	Float_t colH = (Int_t ((rot[1] - col0)/colPadSize)) * colPadSize;
206	hitO[gtb] = (rot[1] - col0) - colH;
207	found_hit = kTRUE;
208	break;
209	}
210
211	if(!found_hit) continue;
212
213	/*
214	printf("gtb: %d, x: %f, rot[0]: %f, Yhit: %f, Ycl: %f\n",
215	gtb, x, rot[0], rot[1], clusterY[gtb]);
216	printf("\n Zhit - Zcl = %f - %f = %f\n",
217	rot[2], clusterZ[gtb], rot[2] - clusterZ[gtb]);
218
219
220
221	printf("found hit within dx = %f - %f \n",rot[0],x);
222	printf("pos: %f, %f, %f \n",pos[0],pos[1],pos[2]);
223	printf("rot: %f, %f, %f \n",rot[0],rot[1],rot[2]);
224	printf("cluster: %d, %f, %f \n",gtb,clusterY[gtb],clusterZ[gtb]);
225	*/
226
227
228	hy->Fill(hitY[gtb]-clusterY[gtb]);
229	if(charge > 0) hyp->Fill(hitY[gtb]-clusterY[gtb]);
230	else if(charge < 0) hym->Fill(hitY[gtb]-clusterY[gtb]);
231
232	if((clusterQ[gtb]>10)&&(clusterSigmaY[gtb]>0))
233	hyn->Fill((hitY[gtb]-clusterY[gtb])/clusterSigmaY[gtb]);
234	hz->Fill(hitZ[gtb]-clusterZ[gtb]);
235	hy2->Fill(hitY[gtb]-clusterY[gtb],clusterQ[gtb]);
236	hy3->Fill(hitY[gtb]-clusterY[gtb],hitO[gtb]);
237	hy4->Fill(hitY[gtb]-clusterY[gtb],(Float_t)(tbpp - 1 - gtb%tbpp));
238	herr->Fill(hitY[gtb]-clusterY[gtb],clusterSigmaY[gtb]);
239	}
240	}
241	}
242
243	gStyle->SetOptStat(1);
244	gStyle->SetOptFit(1);
245
246	TCanvas* c = new TCanvas("c", "c", 110, 110, 810, 840);
247	c->SetFillColor(10);
248	c->Divide(2,2);
249	c->cd(1); hy->SetLineWidth(2); hy->SetFillColor(29); hy->Draw();
250	c->cd(2); hz->SetLineWidth(2); hz->SetFillColor(29); hz->Draw();
251	c->cd(3); hyn->Draw();
252	c->cd(4); hy4->Draw();
253
254	TCanvas* c1 = new TCanvas("c1", "c1", 210, 210, 910, 940);
255	c1->SetFillColor(10);
256	c1->Divide(2,2);
257	c1->cd(1); hyp->SetLineWidth(2); hyp->SetFillColor(29); hyp->Fit("gaus");
258	c1->cd(2); hym->SetLineWidth(2); hym->SetFillColor(29); hym->Fit("gaus");
259	c1->cd(3); hy3->Draw();
260	c1->cd(4); hy2->Draw();
261
262	}
263
264
265
266
267
268