]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/ITShitOccupancy.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
General updates
[u/mrichter/AliRoot.git] / ITS / ITShitOccupancy.C
CommitLineData
58005f18 1/*
2// Some time it is nice to compile the code to find all of the typos.
3#include <stdio.h>
4#include <math.h>
5#include "TROOT.h"
6#include "TSystem.h"
7#include "TClassTable.h"
8#include "TFile.h"
9#include "AliITS.h"
10#include "TTree.h"
11#include "TClonesArray.h"
12#include "TArray.h"
13#include "TCanvas.h"
14#include "AliRun.h"
15#include "AliITSgeom.h"
16#include "TParticlet.h"
17
18extern TSystem *gSystem;
19extern TROOT *gROOT;
20extern TClassTable *gClassTable;
21extern AliRun *gAlice;
22*/
23
24Double_t dEtadZ(Double_t z,Double_t r){
25 Float_t a;
26
27 a = TMath::Sqrt(z*z+r*r);
28 if(a!=0.0) a = 1./a;
29 return(a);
30}
31
32Double_t PdEtadZ(Double_t *z,Double_t *r){
33 return(dEtadZ(*z,*r));
34}
35
36Double_t ZfromEtaR(Double_t eta,Double_t r){
37 Double_t z;
38
39 z = 2.0*TMath::ATan(TMath::Exp(-eta));
40 z = r/TMath::Tan(z);
41 return(z);
42}
43
44Double_t EtafromZr(Double_t z,Double_t r){
45 Double_t a,eta;
46
47 a = 0.5*TMath::ATan2(r,z); // Tan2(y,x) = y/x
48 eta = -TMath::Log(TMath::Abs(TMath::Tan(a)));
49 if(a<0.0) eta = -eta;
50 return (eta);
51}
52
53Double_t PEtafromZr(Double_t *z,Double_t *r){
54 return (EtafromZr(*z,*r));
55}
56
57Double_t Weight(Double_t *a,Double_t *b){
58 Double_t eta,z,Etap,rt,dEta,dEtap,weight,REtap;
59
60 eta = a[0];
61 Etap = a[1];
62 rt = b[0];
63 dEta = b[1];
64 dEtap = b[2];
65 REtap = b[3];
66 z = ZfromEtaR(eta,rt);
67 weight = TMath::Abs(dEtadZ(z,rt));
68 weight = weight*REtap/(rt*2.0*TMath::Pi()*dEta*dEtap);
69 return weight;
70}
71
72void ITShitOccupancy (const char *filename="galice.root",Int_t evNumber=0){
73/////////////////////////////////////////////////////////////////////////
74//
75// Written by Roberto Barbera
76// Modified by Bjorn S. Nilsen May 19 1999
77//
78// This macro is a small example of a ROOT macro
79// illustrating how to read the output of aliroot
80// and fill some histograms. The Occupancy values produce are only
81// approximate and are known to have problems. No warranty of any kind
82// may be give or assumed.
83//
84// Root > .L ITShitOccupancy.C//this loads the macro in memory
85// Root > ITShitOccupancy(); //by default process first event
86// Root > ITShitOccupancy("galiceA.root"); //process file galiceA.root
87// Root > ITShitOccupancy("galiceA.root",2);//process file galiceA.root
88// third event
89//
90// Root > .x ITShitOccupancy.C(); //by default process first event
91// Root > .x ITShitOccupancy.C("galiceA.root"); //process file galiceA.root
92// Root > .x ITShitOccupancy.C("galiceA.root",2);//process file
93// galiceA.root third event
94// The test figures for this macro are kept in
95// $ALICE_ROOT/ITS/figures/ITShitOccupancy.figures as a series of post-
96// script files and one ITShitOccupancy.root file. This is because ther
97// are so many figures produced.
98//
99////////////////////////////////////////////////////////////////////////
100//
101// The code computes the Occupancy of HITS in the different ITS detectors
102// as a function of z. This does not represent the true occupancy expected
103// in the ITS since there will be more than one hit per charge cluster in
104// the different detector of the ITS. To do a proper occupancy calculation
105// a full simulation of all of the detectors in the ITS will be needed. In
106// addition a proper background needs to be used.
107//
108// This code has been modified to be as compatible with the standard (May 26
109// 1999) version of the official released ITS code as possible. Functions
110// introduce in the unofficial release of the ITS code (March 199) have
111// been commented out and replaced by the variables defined in the ITShits
112// class. As soon as these functions are officially part of the aliroot
113// release. These functions should be used instead to maximise future
114// portability.
115//
116// As with all ITS software produced under the deadline of the ITS TDR, there
117// is no real guarantee that this code will function, produce meaning full
118// results, solve in any way any possible problem or question, or be supported
119// in any kind. Questions or comments can be addressed to
120// Nilsen@mps.ohio-state.edu (for as long as I am there) or
121// Roberto.Barbera@ct.infn.it. A response may never come.
122//
123// This exclamation should be copy written and thereby not used my any other
124// software writer, especially those writing code for ALICE. We expect
125// everyone else's code to work perfectly and do the job as advertised, either
126// explicitly, by word of mouth, or by any other mystical means.
127//
128////////////////////////////////////////////////////////////////////////
129// Dynamically link some shared libs
130 if (gClassTable->GetID("AliRun") < 0) {
131 gROOT->LoadMacro("loadlibs.C");
132 loadlibs();
133 } // end if gClassTable...
134//
135// Connect the Root Galice file containing Geometry, Kine and Hits
136 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
137 if (!file) file = new TFile(filename);
138 printf("Reading from root file %s\n",filename);
139//
140// Get AliRun object from file or create it if not on file
141 if (!gAlice) {
142 gAlice = (AliRun*)file->Get("gAlice");
143 if (gAlice) printf("AliRun object found on file\n");
144 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
145 } // end if !gAlice
146//
147// Get pointers to Alice detectors and Hits containers
148 Int_t nparticles = gAlice->GetEvent(evNumber);
149 if (nparticles <= 0) return;
150 AliITS *ITS = gAlice->GetDetector("ITS");
151 if(!ITS) {printf("ITS detector not found. Exiting\n"); return;}
152 TClonesArray *Particles = gAlice->Particles();
153 TClonesArray *ITShits = ITS->Hits();
154 TTree *TH = gAlice->TreeH();
155 Int_t ntracks = TH->GetEntries();
156 AliITSgeom *gm = ITS->GetITSgeom();
157 Int_t Nlayers = gm->GetNlayers();
158//
159 printf("%d primary tracks and their secondary tracks read in\n",ntracks);
160//
161// Import the Kine and Hits Trees for the event evNumber in the file
162 Float_t x,y,z,mass,e,r,rt;
163 Float_t px,py,pz,pt;
164 Float_t destep;
165 Float_t theta,phi;
166 Float_t eta,etap;
167 Float_t dzeta,weight;
168 Int_t idpart;
169//
170 Float_t ztot[] = {33.52,// total length of layer 1
171 33.52,// total length of layer 2
172 46.10,// total length of layer 3
173 60.77,// total length of layer 4
174 90.42,// total length of layer 5
175 101.95};// total length of layer 6
176 Float_t raprx[] = {4.0,7.0,14.9,23.8,39.1,43.6}; // aproximate layer radii
177//
178 Float_t totdet[]={256*256*4*20,//total number of channels layer 1 = 5242880
179 256*256*4*40,//total number of channels layer 2 = 10485760
180 2*256*256*6*14,//total number of channels layer 3 = 11010048
181 2*256*256*8*22,//total number of channels layer 4 = 23068672
182 768*23*34,//total number of channels layer 5 = 1201152
183 768*26*38};//total number of channels layer 6 = 1517568
184//
185 Float_t clusize,clusize_z,clusize_rphi;
186//
187 Float_t threshold[]={8.0e-6,8.0e-6, // noise= 200 e- ; thr = 2000 e- ;
188 3.0e-6,3.0e-6, // noise= 250 e- ; thr= 3 sigma noise
189 6.0e-6,6.0e-6};// noise= 500 e- ; thr= 3 sigma noise
190//
191 Int_t nhitlayer[6];
192 Int_t ntrklayer[6];
193 Int_t layer,ladder,det;
194 Int_t nbytes = 0;
195 Int_t i,j,jj,hit,ipart,ipartold,iparent;
196 Int_t nhits;
197 Int_t sector,plane;
198 TParticle *particle,*pparticle;
199 AliITShit *itsHit;
200
201// Create histograms
202
203 Float_t etamin = -1.;
204 Float_t etamax = 1.;
205 Float_t ptmin = 0.;
206 Float_t ptmax = 3.;
207 Int_t nbin = 40;
208 Int_t nbinz = 200;
209 Float_t deta = (etamax-etamin)/nbin;
210 Float_t EtapMax=2.436246,EtapMin=-2.436246; // from Config.C file
211 TVector etabin(nbin+1);
212 TVector zbin(nbin+1);
213
214// Open a root file to save histograms in
215 char *Hfilename = "ITShitOccupancy.root";
216 TFile *Hfile = (TFile *)gROOT->GetListOfFiles()->FindObject(Hfilename);
217 if(Hfile) Hfile->Close();
218 Hfile = new TFile(Hfilename,"RECREATE","Histograms for ITS Occupancy");
219 printf("Writting histograms to root file %s\n",Hfilename);
220//
221 TH2F *hzetap[6];
222 TH1D *hzetap_z[6],*hzetap_etap[6];
223 TH2F *hzetapden[6];
224 TH1D *hzetapden_z[6],*hzetapden_etap[6];
225 TH2F *hEtaEtapden[6];
226 TH1D *hEtaEtapden_Eta[6],*hEtaEtapden_Etap[6];
227 TH2F *hEtaEtapdenWeight[6];
228 TH1D *hEtaEtapdenWeight_Eta[6];
229 TF2 *fweight;
230 TF1 *fdEtadZ,*fEta;
231 TH1F *hdEtadZ[6],*hEta[6];
232
233 { // "Book" histograms
234 char histid[7],htit[40];
235 for(i=0;i<Nlayers;i++){
236 sprintf(histid,"hzetap%1.1d",i+1);
237 sprintf(htit,"Hits for Layer %1.1d",i+1);
238 hzetap[i] = new TH2F(histid,htit,nbinz,-0.6*ztot[i],0.6*ztot[i],
239 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
240 hzetap[i]->SetXTitle("Z (cm)");
241 hzetap[i]->SetYTitle("Beam particle pseudorapidity");
242 hzetap[i]->Sumw2();
243 //
244 sprintf(histid,"hzetapden%1.1d",i+1);
245 sprintf(htit,"Density of Hits for Layer %1.1d",i+1);
246 hzetapden[i] = new TH2F(histid,htit,nbinz,-0.6*ztot[i],0.6*ztot[i],
247 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
248 hzetapden[i]->Sumw2();
249 hzetapden[i]->SetXTitle("Z (cm)");
250 hzetapden[i]->SetYTitle("Beam particle pseudorapidity");
251 nhitlayer[i] = 0.0;
252 ntrklayer[i] = 0.0;
253 //
254 sprintf(histid,"hEtaEtapden%1.1d",i+1);
255 sprintf(htit,"Hits for Layer %1.1d by Eta",i+1);
256 hEtaEtapden[i] = new TH2F(histid,htit,nbinz,
257 -EtafromZr(0.6*ztot[i],raprx[i]),
258 EtafromZr(0.6*ztot[i],raprx[i]),
259 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
260 hEtaEtapden[i]->SetXTitle("Track Pseudorapidity");
261 hEtaEtapden[i]->SetYTitle("Beam particle pseudorapidity");
262 hEtaEtapden[i]->Sumw2();
263 } // end for i
264
265 sprintf(histid,"fweight");
266 fweight = new TF2(histid,Weight,-EtafromZr(0.6*ztot[0],raprx[0]),
267 EtafromZr(0.6*ztot[0],raprx[0]),
268 floor(EtapMin)-1.0,ceil(EtapMax)+1.0,4);
269 Double_t params[4];
270 sprintf(histid,"fdEtadZ");
271 fdEtadZ = new TF1(histid,PdEtadZ,floor(EtapMin)-1.0,
272 ceil(EtapMax)+1.0,1);
273 sprintf(histid,"fEta");
274 fEta = new TF1(histid,PEtafromZr,floor(EtapMin)-1.0,
275 ceil(EtapMax)+1.0,1);
276 for(i=0;i<Nlayers;i++){
277 params[0] = raprx[i];
278 params[1] = hEtaEtapden[i]->GetXaxis()->GetBinWidth(1);
279 params[2] = hEtaEtapden[i]->GetYaxis()->GetBinWidth(1);
280 params[3] = EtapMax-EtapMin;
281 fweight->SetParameters(params);
282 sprintf(histid,"hEtaEtapdenWeight%1.1d",i+1);
283 sprintf(htit,"Weight for Layer %1.1d by Eta",i+1);
284 hEtaEtapdenWeight[i] = new TH2F(histid,htit,nbinz,
285 -EtafromZr(0.6*ztot[i],raprx[i]),
286 EtafromZr(0.6*ztot[i],raprx[i]),
287 20,floor(EtapMin)-1.0,ceil(EtapMax)+1.0);
288 hEtaEtapdenWeight[i]->SetXTitle("Track Pseudorapidity");
289 hEtaEtapdenWeight[i]->SetYTitle("Beam particle pseudorapidity");
290 hEtaEtapdenWeight[i]->Sumw2();
291 hEtaEtapdenWeight[i]->Eval(fweight);
292//
293 fdEtadZ->SetParameters(params);
294 sprintf(histid,"hdEtadZ%1.1d",i+1);
295 sprintf(htit,"dEtadZ(z) for Layer %1.1d",i+1);
296 hdEtadZ[i] = new TH1F(histid,htit,nbinz,
297 -0.6*ztot[i],0.6*ztot[i]);
298 hdEtadZ[i]->SetXTitle("Z (cm)");
299 hdEtadZ[i]->Sumw2();
300 hdEtadZ[i]->Eval(fdEtadZ);
301//
302 fEta->SetParameters(params);
303 sprintf(histid,"hEta%1.1d",i+1);
304 sprintf(htit,"Eta(z) for Layer %1.1d",i+1);
305 hEta[i] = new TH1F(histid,htit,nbinz,
306 -0.6*ztot[i],0.6*ztot[i]);
307 hEta[i]->SetXTitle("Z (cm)");
308 hEta[i]->Sumw2();
309 hEta[i]->Eval(fEta);
310 } // end for i
311
312 } // end "Book" histograms
313// Start loop on tracks in the hits container
314
315 Double_t xa[2],par[4];
316
317 for (Int_t track=0; track<ntracks;track++) {
318 gAlice->ResetHits();
319 nbytes += TH->GetEvent(track);
320 nhits = ITShits->GetEntriesFast();
321 for (hit=0;hit<nhits;hit++) {
322 itsHit = (AliITShit*)ITShits->UncheckedAt(hit);
323 destep = itsHit->GetIonization();
324 // With this new version, to be able to do proper detector
325 // simulations, the requirment that a "hit" leave some
326 // energy deposited has been removed.
327 if(destep<=0.0) continue; // skip hits without energy loss.
328 ipart = itsHit->fTrack;
329 particle = (TParticle*)Particles->UncheckedAt(ipart);
330 iparent = particle->GetFirstMother();
331 pparticle= (TParticle*)Particles->UncheckedAt(ntracks-track);
332 etap = TMath::ATan2(pparticle->Pt(),pparticle->Pz());
333 etap = -TMath::Log(TMath::Tan(0.5*etap));
334 itsHit->GetPositionG(x,y,z);
335 x -= pparticle->Vx(); // correct for primary vertex position
336 y -= pparticle->Vy(); // correct for primary vertex position
337 z -= pparticle->Vz(); // correct for primary vertex position
338 itsHit->GetMomentumG(px,py,pz);
339 itsHit->GetDetectorID(layer,ladder,det);
340 r = sqrt(x*x+y*y+z*z);
341 rt = sqrt(x*x+y*y);
342 theta = TMath::ACos(z/r)*180./TMath::Pi();
343 phi = TMath::ATan2(y,x); if(phi<0.0) phi += TMath::Pi();
344 eta = EtafromZr(z,rt);
345 pt = TMath::Sqrt(px*px+py*py);
346 if(ipart!=ipartold)ntrklayer[layer-1] += 1;
347 ipartold = ipart;
348 nhitlayer[layer-1] += 1;
349
350 switch (layer) {
351 // Calculate average cluster sizes
352 case 1:
353 clusize_rphi = 1.4;
354 if(TMath::Abs(eta)>=0.0 && TMath::Abs(eta)< 0.2) clusize_z=1.2;
355 if(TMath::Abs(eta)>=0.2 && TMath::Abs(eta)< 0.4) clusize_z=1.2;
356 if(TMath::Abs(eta)>=0.4 && TMath::Abs(eta)< 0.6) clusize_z=1.3;
357 if(TMath::Abs(eta)>=0.6 && TMath::Abs(eta)< 0.8) clusize_z=1.4;
358 if(TMath::Abs(eta)>=0.8 && TMath::Abs(eta)< 1.0) clusize_z=1.5;
359 clusize = clusize_z*clusize_rphi;
360 case 2:
361 clusize_rphi = 1.4;
362 if(TMath::Abs(eta)>=0.0 && TMath::Abs(eta)< 0.2) clusize_z=1.2;
363 if(TMath::Abs(eta)>=0.2 && TMath::Abs(eta)< 0.4) clusize_z=1.2;
364 if(TMath::Abs(eta)>=0.4 && TMath::Abs(eta)< 0.6) clusize_z=1.3;
365 if(TMath::Abs(eta)>=0.6 && TMath::Abs(eta)< 0.8) clusize_z=1.4;
366 if(TMath::Abs(eta)>=0.8 && TMath::Abs(eta)< 1.0) clusize_z=1.5;
367 clusize = clusize_z*clusize_rphi;
368 case 3:
369 clusize_rphi = 8.0;
370 clusize_z = 3.0;
371 clusize = clusize_z*clusize_rphi;
372 case 4:
373 clusize_rphi = 8.0;
374 clusize_z = 3.0;
375 clusize = clusize_z*clusize_rphi;
376 case 5:
377 clusize_z = 1.0;
378 if(pt>=0.0 && pt< 0.2) clusize_rphi=2.0;
379 if(pt>=0.2 && pt< 0.4) clusize_rphi=1.9;
380 if(pt>=0.4 && pt< 0.6) clusize_rphi=1.6;
381 if(pt>=0.6 && pt< 0.8) clusize_rphi=1.6;
382 if(pt>=0.8 && pt< 1.0) clusize_rphi=1.7;
383 if(pt>=1.0 && pt< 1.5) clusize_rphi=1.7;
384 if(pt>=1.5 && pt< 2.5) clusize_rphi=1.7;
385 if(pt>=2.5 ) clusize_rphi=1.7;
386 clusize = clusize_z*clusize_rphi;
387 case 6:
388 clusize_z = 1.0;
389 if(pt>=0.0 && pt< 0.2) clusize_rphi=2.0;
390 if(pt>=0.2 && pt< 0.4) clusize_rphi=1.9;
391 if(pt>=0.4 && pt< 0.6) clusize_rphi=1.6;
392 if(pt>=0.6 && pt< 0.8) clusize_rphi=1.6;
393 if(pt>=0.8 && pt< 1.0) clusize_rphi=1.7;
394 if(pt>=1.0 && pt< 1.5) clusize_rphi=1.7;
395 if(pt>=1.5 && pt< 2.5) clusize_rphi=1.7;
396 if(pt>=2.5 ) clusize_rphi=1.7;
397 clusize = clusize_z*clusize_rphi;
398 } // end switch layer
399// weightE = clusize*ztot[layer-1]/totdet[layer-1];
400 xa[0] = eta; xa[1] = etap;
401 par[0] = rt;
402 par[1] = hEtaEtapden[layer-1]->GetXaxis()->GetBinWidth(1);
403 par[2] = hEtaEtapden[layer-1]->GetYaxis()->GetBinWidth(1);
404 par[3] = EtapMax-EtapMin;
405 weight = Weight(xa,par);
406 hEtaEtapden[layer-1]->Fill(eta,etap,weight);
407 weight = 1.;
408 hzetap[layer-1]->Fill(z,etap,weight);
409 weight = par[3]/(rt*2.0*TMath::Pi()*
410 hzetapden[layer-1]->GetXaxis()->GetBinWidth(1)*
411 hzetapden[layer-1]->GetYaxis()->GetBinWidth(1));
412 hzetapden[layer-1]->Fill(z,etap,weight);
413 } // end for hit
414 } // end for track
415
416 for(i=0;i<Nlayers;i++)printf("No. of hits in layer %d = %d\n",i+1,nhitlayer[i]);
417 for(i=0;i<Nlayers;i++)printf("No. of tracks in layer %d = %d\n",i+1,ntrklayer[i]);
418
419//Create canvases, set the view ranges and show histograms
420
421 TCanvas *canvas = new TCanvas("canvas","ITShitOccupancy",1);
422
423 Bool_t printit = kTRUE;
424 char psfilename[40];
425 for(i=0;i<Nlayers;i++){
426 hzetap[i]->Draw("colz");
427 sprintf(psfilename,"ITShitOccupancy_%1.1d_z_etap.ps",i+1);
428 if(printit) canvas->Print(psfilename);
429 sprintf(psfilename,"hzetap_z%1.1d",i+1);
430 hzetap_z[i] = hzetap[i]->ProjectionX(psfilename,0,
431 hzetap[i]->GetNbinsY()+1,"E");
432 hzetap_z[i]->SetXTitle("Z (cm)");
433 hzetap_z[i]->SetYTitle("Number of Hits per cm");
434 sprintf(psfilename,"hzetap_etap%1.1d",i+1);
435 hzetap_etap[i] = hzetap[i]->ProjectionY(psfilename,0,
436 hzetap[i]->GetNbinsX()+1,"E");
437 hzetap_etap[i]->SetXTitle("Beam particle pseudorapidity");
438 hzetap_etap[i]->SetYTitle("Number of Hits");
439 hzetap_z[i]->Draw();
440 sprintf(psfilename,"ITShitOccupancy_%1.1d_z.ps",i+1);
441 if(printit) canvas->Print(psfilename);
442 hzetap_etap[i]->Draw();
443 sprintf(psfilename,"ITShitOccupancy_%1.1d_etap.ps",i+1);
444 if(printit) canvas->Print(psfilename);
445//
446 hzetapden[i]->Draw("colz");
447 sprintf(psfilename,"ITSHitDen_%1.1d_z_etap.ps",i+1);
448 if(printit) canvas->Print(psfilename);
449 sprintf(psfilename,"hzetapden_z%1.1d",i+1);
450 hzetapden_z[i] = hzetapden[i]->ProjectionX(psfilename,0,
451 hzetapden[i]->GetNbinsY()+1,"E");
452 hzetapden_z[i]->SetXTitle("Z (cm)");
453 hzetapden_z[i]->SetYTitle("Number of Hits (cm^-2)");
454 sprintf(psfilename,"hzetapden_etap%1.1d",i+1);
455 hzetapden_etap[i] = hzetapden[i]->ProjectionY(psfilename,0,
456 hzetapden[i]->GetNbinsX()+1,"E");
457 hzetapden_etap[i]->SetXTitle("Beam particle pseudorapidity");
458 hzetapden_etap[i]->SetYTitle("Number of Hits (cm^-2)");
459 hzetapden_z[i]->Draw();
460 sprintf(psfilename,"ITSHitDen_%1.1d_z.ps",i+1);
461 if(printit) canvas->Print(psfilename);
462 hzetapden_etap[i]->Draw();
463 sprintf(psfilename,"ITSHitDen_%1.1d_etap.ps",i+1);
464 if(printit) canvas->Print(psfilename);
465//
466 hEtaEtapden[i]->Draw("colz");
467 sprintf(psfilename,"ITShitDensity_%1.1d_Eta_Etap.ps",i+1);
468 if(printit) canvas->Print(psfilename);
469 sprintf(psfilename,"hEtaEtapden_Eta%1.1d",i+1);
470 hEtaEtapden_Eta[i] = hEtaEtapden[i]->ProjectionX(psfilename,0,
471 hEtaEtapden[i]->GetNbinsY()+1,"E");
472 hEtaEtapden_Eta[i]->SetXTitle("Pseudorapidity");
473 hEtaEtapden_Eta[i]->SetYTitle("Number of Hits (cm^-2)");
474 sprintf(psfilename,"hEtaEtapden_Etap%1.1d",i+1);
475 hEtaEtapden_Etap[i] = hEtaEtapden[i]->ProjectionY(psfilename,0,
476 hEtaEtapden[i]->GetNbinsX()+1,"E");
477 hEtaEtapden_Etap[i]->SetXTitle("Beam particle pseudorapidity");
478 hEtaEtapden_Etap[i]->SetYTitle("Number of Hits (cm^-2)");
479 hEtaEtapden_Eta[i]->Draw();
480 sprintf(psfilename,"ITShitDensity_%1.1d_Eta.ps",i+1);
481 if(printit) canvas->Print(psfilename);
482 hEtaEtapden_Etap[i]->Draw();
483 sprintf(psfilename,"ITShitDensity_%1.1d_Etap.ps",i+1);
484 if(printit) canvas->Print(psfilename);
485//
486 hEtaEtapdenWeight[i]->Draw("colz");
487 sprintf(psfilename,"ITShitDensityWeight_%1.1d_Eta_Etap.ps",i+1);
488 if(printit) canvas->Print(psfilename);
489 sprintf(psfilename,"hEtaEtapdenWeight_Eta%1.1d",i+1);
490 hEtaEtapdenWeight_Eta[i] =
491 hEtaEtapdenWeight[i]->ProjectionX(psfilename,
492 0,hEtaEtapdenWeight[i]->GetNbinsY()+1,"E");
493 hEtaEtapdenWeight_Eta[i]->SetXTitle("Pseudorapidity");
494 hEtaEtapdenWeight_Eta[i]->SetYTitle("Weight (cm^-2)");
495 hEtaEtapdenWeight_Eta[i]->Draw();
496 sprintf(psfilename,"ITShitDensityWeight_%1.1d_Eta.ps",i+1);
497 if(printit) canvas->Print(psfilename);
498 } // end for i
499//
500 Hfile->Write();
501//
502 return;
503}