]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/SPDclusterTestDubna.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliGenSTRANGElib.cxx first commit.
[u/mrichter/AliRoot.git] / ITS / SPDclusterTestDubna.C
CommitLineData
409f8c84 1#include "iostream.h"
2
3void SPDclusterTestDubna (Int_t evNumber1=0,Int_t evNumber2=0)
4{
5/////////////////////////////////////////////////////////////////////////
6// This macro is a small example of a ROOT macro
7// illustrating how to read the output of GALICE
8// and do some analysis.
9//
10/////////////////////////////////////////////////////////////////////////
11
12// Dynamically link some shared libs
13
14 if (gClassTable->GetID("AliRun") < 0) {
15 gROOT->LoadMacro("loadlibs.C");
16 loadlibs();
17 } else {
18 delete gAlice;
19 gAlice=0;
20 }
21
22// Connect the Root Galice file containing Geometry, Kine and Hits
23
24 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
25 if (!file) file = new TFile("galice.root");
26 file->ls();
27
28// Get AliRun object from file or create it if not on file
29
30 if (!gAlice) {
31 gAlice = (AliRun*)file->Get("gAlice");
32 if (gAlice) printf("AliRun object found on file\n");
33 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
34 }
35
36 // ------------ Cluster and point analysis histogramms ------------
37
38 TH1F *Nxpix1 = new TH1F("Nxpix1","Cluster size in x(r*phi) direction for layer 1",20,0.,20.);
39 TH1F *Nxpix2 = new TH1F("Nxpix2","Cluster size in x(r*phi) direction for layer 2",20,0.,20.);
40 TH1F *Nzpix1 = new TH1F("Nzpix1","Cluster size in z direction for layer 1",15,0.,15.);
41 TH1F *Nzpix2 = new TH1F("Nzpix2","Cluster size in z direction for layer 2",15,0.,15.);
42 TH1F *Xpix1 = new TH1F("Xpix1","Local x coordinate (mm) for layer 1",20,-2.,18.);
43 TH1F *Xpix2 = new TH1F("Xpix2","Local x coordinate (mm) for layer 2",20,-2.,18.);
44 TH1F *Zpix1 = new TH1F("Zpix1","Local z coordinate (mm) for layer 1",90,-2.,88.);
45 TH1F *Zpix2 = new TH1F("Zpix2","Lolac z coordinate (mm) for layer 2",90,-2.,88.);
46
47 TH1F *Xres1 = new TH1F("Xres1","Xrec and Xgen difference (micr) for layers 1",100,-200.,200.);
48 TH1F *Xres2 = new TH1F("Xres2","Xrec and Xgen difference (micr) for layers 2",100,-200.,200.);
49 TH1F *Zres1 = new TH1F("Zres1","Zrec and Zgen difference (micr) for layers 1",100,-800.,800.);
50 TH1F *Zres2 = new TH1F("Zres2","Zrec and Zgen difference (micr) for layers 2",100,-800.,800.);
51
52 TH1F *Ptot1 = new TH1F("Ptot1","Total momentum (GeV/C) for layers 1",100,0.,5.);
53 TH1F *Pz1 = new TH1F("Pz1","Pz (GeV/C) for layers 1",100,-5.,5.);
54 TH1F *Theta1 = new TH1F("Theta1","Theta angle (rad) for layers 1",100,0.,4.);
55 TH1F *Y1 = new TH1F("Y1","Rapidity for layers 1",100,-4.,4.);
56 TH1F *Eta1 = new TH1F("Eta1","PseudoRapidity for layers 1",100,-4.,4.);
57 TH1F *Y1Den = new TH1F("Y1Den","Rapidity for layers 1",100,-0.5,0.5);
58 TH1F *Eta1Den = new TH1F("Eta1Den","PseudoRapidity for layers 1",100,-0.5,0.5);
59 TH1F *Y1DenA = new TH1F("Y1DenA","Rapidity for layers 1",100,-0.5,0.5);
60 TH1F *Eta1DenA = new TH1F("Eta1DenA","PseudoRapidity for layers 1",100,-0.5,0.5);
61 TH1F *Phi1 = new TH1F("Phi1","Phi angle (rad) for layers 1",100,0.,7.);
62 TH1F *Ptot2 = new TH1F("Ptot2","Total momentum (GeV/C) for layers 2",100,0.,5.);
63 TH1F *Pz2 = new TH1F("Pz2","Pz (GeV/C) for layers 2",100,-5.,5.);
64 TH1F *Theta2 = new TH1F("Theta2","Theta angle (rad) for layers 2",100,0.,4.);
65 TH1F *Y2 = new TH1F("Y2","Rapidity for layers 2",100,-4.,4.);
66 TH1F *Eta2 = new TH1F("Eta2","PseudoRapidity for layers 2",100,-4.,4.);
67 TH1F *Y2Den = new TH1F("Y2Den","Rapidity for layers 2",100,-0.5,0.5);
68 TH1F *Eta2Den = new TH1F("Eta2Den","PseudoRapidity for layers 2",100,-0.5,0.5);
69 TH1F *Y2DenA = new TH1F("Y2DenA","Rapidity for layers 2",100,-0.5,0.5);
70 TH1F *Eta2DenA = new TH1F("Eta2DenA","PseudoRapidity for layers 2",100,-0.5,0.5);
71 TH1F *Phi2 = new TH1F("Phi2","Phi angle (rad) for layers 2",100,0.,7.);
72
73 // -------------- Create ntuples --------------------
74 // ntuple structures:
75
76 struct {
77 Int_t lay;
78 Int_t nx;
79 Int_t nz;
80 Int_t hitprim;
81 Int_t partcode;
82 Int_t ntrover;
83 Float_t dx;
84 Float_t dz;
85 Float_t pmod;
86 } ntuple_st;
87
88 struct {
89 Int_t lay;
90 Int_t lad;
91 Int_t det;
92 Int_t nx;
93 Int_t nz;
94 Int_t ntrover;
95 Int_t noverlaps;
96 Int_t noverprim;
97 Float_t qcl;
98 Float_t dx;
99 Float_t dz;
100 Float_t x;
101 Float_t z;
102 } ntuple1_st;
103
104 struct {
105 // Int_t lay;
106 Int_t lay;
107 Int_t nx;
108 Int_t nz;
109 Float_t x;
110 Float_t z;
111 Float_t qcl;
112 } ntuple2_st;
113
114 ntuple = new TTree("ntuple","Demo ntuple");
115 ntuple->Branch("lay",&ntuple_st.lay,"lay/I");
116 ntuple->Branch("nx",&ntuple_st.nx,"nx/I");
117 ntuple->Branch("nz",&ntuple_st.nz,"nz/I");
118 ntuple->Branch("hitprim",&ntuple_st.hitprim,"hitprim/I");
119 ntuple->Branch("partcode",&ntuple_st.partcode,"partcode/I");
120 ntuple->Branch("ntrover",&ntuple_st.ntrover,"ntrover/I");
121 ntuple->Branch("dx",&ntuple_st.dx,"dx/F");
122 ntuple->Branch("dz",&ntuple_st.dz,"dz/F");
123 ntuple->Branch("pmod",&ntuple_st.pmod,"pmod/F");
124
125 ntuple1 = new TTree("ntuple1","Demo ntuple1");
126 ntuple1->Branch("lay",&ntuple1_st.lay,"lay/I");
127 ntuple1->Branch("lad",&ntuple1_st.lad,"lad/I");
128 ntuple1->Branch("det",&ntuple1_st.det,"det/I");
129 ntuple1->Branch("nx",&ntuple1_st.nx,"nx/I");
130 ntuple1->Branch("nz",&ntuple1_st.nz,"nz/I");
131 ntuple1->Branch("qcl",&ntuple1_st.qcl,"qcl/F");
132 ntuple1->Branch("ntrover",&ntuple1_st.ntrover,"ntrover/I");
133 ntuple1->Branch("noverlaps",&ntuple1_st.noverlaps,"noverlaps/I");
134 ntuple1->Branch("noverprim",&ntuple1_st.noverprim,"noverprim/I");
135 ntuple1->Branch("x",&ntuple1_st.x,"x/F");
136 ntuple1->Branch("z",&ntuple1_st.z,"z/F");
137 ntuple1->Branch("dx",&ntuple1_st.dx,"dx/F");
138 ntuple1->Branch("dz",&ntuple1_st.dz,"dz/F");
139
140
141 ntuple2 = new TTree("ntuple2","Demo ntuple2");
142 // ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
143 ntuple2->Branch("lay",&ntuple2_st.lay,"lay/I");
144 ntuple2->Branch("x",&ntuple2_st.x,"x/F");
145 ntuple2->Branch("z",&ntuple2_st.z,"z/F");
146 ntuple2->Branch("nx",&ntuple2_st.nx,"nx/I");
147 ntuple2->Branch("nz",&ntuple2_st.nz,"nz/I");
148 ntuple2->Branch("qcl",&ntuple2_st.qcl,"qcl/F");
149
150// ------------------------------------------------------------------------
151//
152// Loop over events
153//
154 for (int nev=0; nev<= evNumber2; nev++) {
155 Int_t nparticles = gAlice->GetEvent(nev);
156 cout << "nev " << nev <<endl;
157 cout << "nparticles " << nparticles <<endl;
158 if (nev < evNumber1) continue;
159 if (nparticles <= 0) return;
160
161 TTree *TH = gAlice->TreeH();
162 Int_t ntracks = TH->GetEntries();
163 cout<<"ntracks "<<ntracks<<endl;
164
165// Get pointers to Alice detectors and Digits containers
166 AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
167 TClonesArray *Particles = gAlice->Particles();
168
169 if (ITS) {
170 // fill modules with sorted by module hits
171 Int_t nmodules;
172 ITS->InitModules(-1,nmodules);
173 // ITS->FillModules(nev,-1,evNumber2,nmodules," "," ");
174 ITS->FillModules(nev,evNumber2,nmodules," "," ");
175 //get pointer to modules array
176 TObjArray *ITSmodules = ITS->GetModules();
177 AliITShit *itsHit;
178
179 // get the Tree for clusters
180 ITS->GetTreeC(nev);
181 TTree *TC=ITS->TreeC();
182 Int_t nent=TC->GetEntries();
183 printf("Found %d entries in the tree (must be one per module per event!)\n",nent);
184 Int_t lay, lad, det;
185 AliITSgeom *geom = ITS->GetITSgeom();
338e4f06 186
187 AliITSDetType *iDetType=ITS->DetType(0);
188 AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
189 printf("SPD dimensions %f %f %f \n",seg0->Dx(),seg0->Dz(),seg0->Dy());
190 printf("SPD npixels %d %d \n",seg0->Npz(),seg0->Npx());
191 printf("SPD pitches %d %d \n",seg0->Dpz(0),seg0->Dpx(0));
192 Float_t SPDlength = seg0->Dz();
193 Float_t SPDwidth = seg0->Dx();
194 Float_t SPDthickness = seg0->Dy();
195 Float_t xpitch = seg0->Dpx(0);
196 if(SPDlength > 80000) SPDthickness = 150;
197 Float_t ylim;
198 if(SPDthickness < 200) {
199 ylim = SPDthickness/2 - 4;
200 }else{
201 ylim = SPDthickness/2 - 10;
202 }
203
409f8c84 204 for (Int_t idettype=0;idettype<3;idettype++) {
205
206 TClonesArray *ITSclusters = ITS->ClustersAddress(idettype);
207 //printf ("ITSclusters %p \n",ITSclusters);
208
209 if (idettype != 0) continue;
210
211 Float_t occup1 = 0;
212 Float_t occup2 = 0;
213
214 // Module loop
215 for (Int_t mod=0; mod<nent; mod++) {
216 AliITSmodule *itsModule = (AliITSmodule*)ITSmodules->At(mod);
217 geom->GetModuleId(mod,lay,lad,det);
218
219 Int_t nhits = itsModule->GetNhits();
220 //if(nhits) printf("module nhits %d %d\n",mod,nhits);
221 if(!nhits) continue;
222
223 ITS->ResetClusters();
224 TC->GetEvent(mod);
225 Int_t nclust = ITSclusters->GetEntries();
226 if (!nclust) continue;
227
228 // cluster/hit loops
229 //cout<<"mod,lay,nclust,nhits ="<<mod<<","<<lay<<","<<nclust<<","<<nhits<<endl;
230 for (Int_t clu=0;clu<nclust;clu++) {
231 itsclu = (AliITSRawClusterSPD*)ITSclusters->UncheckedAt(clu);
232
233 Int_t noverlaps = 0;
234 Int_t noverprim = 0;
235
236 Int_t clustersizex = itsclu->NclX();
237 Int_t clustersizez = itsclu->NclZ();
409f8c84 238 Int_t xstart = itsclu->XStartf();
239 Int_t xstop = itsclu->XStopf();
338e4f06 240 Float_t fxstart = xstart*xpitch;
241 Float_t fxstop = (xstop+1)*xpitch;
409f8c84 242 Float_t zstart = itsclu->ZStart();
243 Float_t zstop = itsclu->ZStop();
244 Int_t zend = itsclu->Zend();
245 Int_t ntrover = itsclu->NTracks();
338e4f06 246 Float_t clusterx = 0;
409f8c84 247 Float_t clusterz = itsclu->Z();
338e4f06 248
249 Int_t tr0, tr1, tr2;
250 itsclu->GetTracks(tr0,tr1,tr2);
251
252 for(Int_t ii=0;ii<clustersizex;ii++) {
253 clusterx += (xstart+0.5+ii)*xpitch;
254 }
255 clusterx /= clustersizex;
256 clusterz /= clustersizez;
257
409f8c84 258 Float_t clusterQ = itsclu->Q();
259
260 if(lay == 1) occup1 += clusterQ;
261 if(lay == 2) occup2 += clusterQ;
262
263 ntuple2_st.lay = lay;
264 ntuple2_st.x = clusterx/1000.;
265 ntuple2_st.z = clusterz/1000.;
266 ntuple2_st.nx = clustersizex;
267 ntuple2_st.nz = clustersizez;
268 ntuple2_st.qcl = clusterQ;
269
270 ntuple2->Fill();
271
272 Int_t icl = 0;
273 Float_t dxprimlast = 10.e+6;
274 Float_t dzprimlast = 10.e+6;
275
338e4f06 276 Float_t xhit0 = 1e+6;
277 Float_t zhit0 = 1e+6;
409f8c84 278 for (Int_t hit=0;hit<nhits;hit++) {
279
280 // Find coordinate differences between the hit and cluster positions
281 // for the resolution determination.
282
283 itsHit = (AliITShit*)itsModule->GetHit(hit);
284
285 Int_t hitlayer = itsHit->GetLayer();
286 Int_t hitladder= itsHit->GetLadder();
287 Int_t hitdet= itsHit->GetDetector();
288 Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
289 Int_t track = itsHit->GetTrack();
290 Int_t dray = 0;
291 Int_t hitstat = itsHit->GetTrackStatus();
409f8c84 292 Float_t zhit = 10000*itsHit->GetZL();
293 Float_t xhit = 10000*itsHit->GetXL();
338e4f06 294 Float_t yhit = 10000*itsHit->GetYL();
295 Int_t parent = itsHit->GetParticle()->GetFirstMother();
296 Int_t partcode = itsHit->GetParticle()->GetPdgCode();
297
298 Int_t hitprim = 0;
299
300 if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
301 // particles
409f8c84 302
303 Float_t pxsimL = itsHit->GetPXL(); // the momenta at GEANT points
304 Float_t pysimL = itsHit->GetPYL();
305 Float_t pzsimL = itsHit->GetPZL();
306 Float_t psimL = TMath::Sqrt(pxsimL*pxsimL+pysimL*pysimL+pzsimL*pzsimL);
307
308 // Check boundaries
309 if(zhit > SPDlength/2) {
310 //cout<<"!!! z outside ="<<zhit<<endl;
311 zhit = SPDlength/2 - 10;
312 }
313 if(zhit < 0 && zhit < -SPDlength/2) {
314 //cout<<"!!! z outside ="<<zhit<<endl;
315 zhit = -SPDlength/2 + 10;
316 }
317 if(xhit > SPDwidth/2) {
318 //cout<<"!!! x outside ="<<xhit<<endl;
319 xhit = SPDwidth/2 - 10;
320 }
321 if(xhit < 0 && xhit < -SPDwidth/2) {
322 //cout<<"!!! x outside ="<<xhit<<endl;
323 xhit = -SPDwidth/2 + 10;
324 }
325
326 zhit += SPDlength/2;
327 xhit += SPDwidth/2;
409f8c84 328
338e4f06 329
330
331 if(hitlayer == 1 && hitstat == 66 && yhit > ylim) {
332 //if(hitstat == 66 && hitprim ==1) {
409f8c84 333 xhit0 = xhit;
334 zhit0 = zhit;
335 }
338e4f06 336
337
338 if(hitlayer == 2 && hitstat == 66 && yhit < -ylim) {
409f8c84 339 xhit0 = xhit;
340 zhit0 = zhit;
338e4f06 341 }
342
409f8c84 343 if(hitstat != 68) continue; // Take only the hit if the last
344 // track point went out from
345 // the detector.
346
338e4f06 347 if(xhit0 > 1e+5 || zhit0 > 1e+5) continue;
409f8c84 348
349 Float_t xmed = (xhit + xhit0)/2;
350 Float_t zmed = (zhit + zhit0)/2;
351
352 Float_t xdif = xmed - clusterx;
353 Float_t zdif = zmed - clusterz;
338e4f06 354 xhit0 = 1e+6;
355 zhit0 = 1e+6;
356
357 // Consider the hits inside of cluster region only b.b.
409f8c84 358
338e4f06 359 if((xmed >= fxstart && xmed <= fxstop) && (zmed >= zstart && zmed <= zstop)) {
409f8c84 360
338e4f06 361 // Consider the hits only with the track number equaled to one
362 // of the cluster
363 //if((track == tr0) || (track == tr1) || (track == tr2)) {
409f8c84 364
409f8c84 365 icl = 1;
366
367 // part = (TParticle *)particles.UncheckedAt(track);
368 // Int_t partcode = part->GetPdgCode();
369 // Int_t primery = gAlice->GetPrimary(track);
370
409f8c84 371
372// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
373// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
374
375 Float_t pmod = itsHit->GetParticle()->P(); // total momentum at the
376 // vertex
377 Float_t energy = itsHit->GetParticle()->Energy(); // energy at the
378 // vertex
379 Float_t mass = itsHit->GetParticle()->GetMass(); // particle mass
380
381 Float_t pz = itsHit->GetParticle()->Pz(); // z momentum componetnt
382 // at the vertex
383 Float_t px = itsHit->GetParticle()->Px(); // z momentum componetnt
384 // at the vertex
385 Float_t py = itsHit->GetParticle()->Py(); // z momentum componetnt
386 // at the vertex
387 Float_t phi = itsHit->GetParticle()->Phi(); // Phi angle at the
388 // vertex
389 Float_t theta = itsHit->GetParticle()->Theta(); // Theta angle at the
390 // vertex
391 //Float_t eta = itsHit->GetParticle()->Eta(); // Pseudo rapidity at the
338e4f06 392 // vertex
393 Float_t y;
409f8c84 394 if((energy-pz) > 0) {
338e4f06 395 y = 0.5*TMath::Log((energy+pz)/(energy-pz));
409f8c84 396 }else{
397 cout<<" Warning: energy < pz ="<<energy<<","<<pz<<endl;
398 y = 10;
399 }
400 Float_t eta = -TMath::Log(TMath::Tan(theta/2));
401 pmod *= 1.0e+3;
402
403
404 Float_t pxsim = itsHit->GetPXG(); // the momenta at this GEANT point
405 Float_t pysim = itsHit->GetPYG();
406 Float_t pzsim = itsHit->GetPZG();
407 Float_t psim = TMath::Sqrt(pxsim*pxsim+pysim*pysim+pzsim*pzsim);
408
338e4f06 409 //Int_t hitprim = 0;
409f8c84 410
411 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
412 // at p < 6 MeV/c
413
414 if(dray == 0) noverlaps = noverlaps + 1; // overlapps for all hits but
415 // not for delta ray which
416 // also went out from the
417 // detector and returned
418 // again
419
338e4f06 420 //if(parent < 0) hitprim = hitprim + 1; // hitprim=1 for the primery
409f8c84 421 // particles
422
423 if(hitprim > 0) noverprim = noverprim + 1;
424
425 if(hitprim > 0) {
426 dxprimlast = xdif;
427 dzprimlast = zdif;
428 }
429
430 // fill ntuple
431
432 ntuple_st.lay = hitlayer;
433 ntuple_st.nx = clustersizex;
434 ntuple_st.nz = clustersizez;
435 ntuple_st.hitprim = hitprim;
436 ntuple_st.partcode = partcode;
437 ntuple_st.ntrover = ntrover;
438 ntuple_st.dx = xdif;
439 ntuple_st.dz = zdif;
440 ntuple_st.pmod = pmod;
441
442 ntuple->Fill();
443
444 if(hitlayer == 1) {
445 Y1DenA->Fill(y);
446 Eta1DenA->Fill(eta);
447 }
448 if(hitlayer == 2) {
449 Y2DenA->Fill(y);
450 Eta2DenA->Fill(eta);
451 }
452
453
454
455 if(hitprim > 0) { // for primary particles
456
457 if(hitlayer == 1) {
458 Xres1->Fill(xdif);
459 Zres1->Fill(zdif);
460 Ptot1->Fill(pmod/1000.);
461 Pz1->Fill(pz);
462 Theta1->Fill(theta);
463 Y1->Fill(y);
464 Eta1->Fill(eta);
465 Y1Den->Fill(y);
466 Eta1Den->Fill(eta);
467 Phi1->Fill(phi);
468 }
469 if(hitlayer == 2) {
470 Xres2->Fill(xdif);
471 Zres2->Fill(zdif);
472 Ptot2->Fill(pmod/1000.);
473 Pz2->Fill(pz);
474 Theta2->Fill(theta);
475 Y2->Fill(y);
476 Eta2->Fill(eta);
477 Y2Den->Fill(y);
478 Eta2Den->Fill(eta);
479 Phi2->Fill(phi);
480 }
481 } // primery particles
482
483 } // end of cluster region
338e4f06 484
409f8c84 485 } // end of hit loop
486
487 if(icl == 1) {
488
489 // fill ntuple1
490
409f8c84 491 if(noverlaps == 0) noverlaps = 1; // cluster contains one or more
492 // delta rays only
493
494 ntuple1_st.lay = lay;
495 ntuple1_st.lad = lad;
496 ntuple1_st.det = det;
497 ntuple1_st.x = clusterx*1000.;
498 ntuple1_st.z = clusterz*1000.;
499 ntuple1_st.nx = clustersizex;
500 ntuple1_st.nz = clustersizez;
501 ntuple1_st.qcl = clusterQ;
502 ntuple1_st.ntrover = ntrover;
503 ntuple1_st.noverlaps = noverlaps;
504 ntuple1_st.noverprim = noverprim;
505 ntuple1_st.dx = dxprimlast;
506 ntuple1_st.dz = dzprimlast;
507
508 ntuple1->Fill();
509
510 } // icl = 1
511 } // cluster loop
512 } // module loop
513
514 cout<<" Occupancy for layer-1 ="<<occup1<<endl;
515 cout<<" Occupancy for layer-2 ="<<occup2<<endl;
516 // The real occupancy values are:
517 // (for full ALICE event at the full SPD acceptence)
518 // occup1 /= 3932160;
519 // occup2 /= 7864320;
520
521 } // idettype loop
522 } // end if ITS
523} // event loop
524
525
526 // Write and Draw Histogramms and ntuples
527
528
529
530 TFile fhistos("SPD_his_dubna.root","RECREATE");
531
532 ntuple->Write();
533 ntuple1->Write();
534 ntuple2->Write();
535
536 Nxpix1->Write();
537 Nzpix1->Write();
538 Nxpix2->Write();
539 Nzpix2->Write();
540
541 Xpix1->Write();
542 Zpix1->Write();
543 Xpix2->Write();
544 Zpix2->Write();
545
546 Xres1->Write();
547 Zres1->Write();
548 Xres2->Write();
549 Zres2->Write();
550
551 Ptot1->Write();
552 Pz1->Write();
553 Theta1->Write();
554 Y1->Write();
555 Eta1->Write();
556 Y1Den->Write();
557 Eta1Den->Write();
558 Y1DenA->Write();
559 Eta1DenA->Write();
560 Phi1->Write();
561
562 Ptot2->Write();
563 Pz2->Write();
564 Theta2->Write();
565 Y2->Write();
566 Eta2->Write();
567 Y2Den->Write();
568 Eta2Den->Write();
569 Y2DenA->Write();
570 Eta2DenA->Write();
571 Phi2->Write();
572
573 fhistos.Close();
574 cout<<"!!! Histogramms and ntuples were written"<<endl;
575
576 TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
577 c1->Divide(2,2);
578 c1->cd(1);
579 gPad->SetFillColor(33);
580 Xres1->SetFillColor(42);
581 Xres1->Draw();
582 c1->cd(2);
583 gPad->SetFillColor(33);
584 Zres1->SetFillColor(46);
585 Zres1->Draw();
586 c1->cd(3);
587 gPad->SetFillColor(33);
588 Xres2->SetFillColor(42);
589 Xres2->Draw();
590 c1->cd(4);
591 gPad->SetFillColor(33);
592 Zres2->SetFillColor(46);
593 Zres2->Draw();
594
595 cout<<"END test for clusters and hits "<<endl;
596
597// file->Close();
598}
599
600
601