]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - RICH/AliRICH.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Print() added
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17 $Log$
18 Revision 1.56 2001/11/02 15:37:25 hristov
19 Digitizer class created. Code cleaning and bug fixes (J.Chudoba)
20
21 Revision 1.55 2001/10/23 13:03:35 hristov
22 The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)
23
24 Revision 1.54 2001/09/07 08:38:10 hristov
25 Pointers initialised to 0 in the default constructors
26
27 Revision 1.53 2001/08/30 09:51:23 hristov
28 The operator[] is replaced by At() or AddAt() in case of TObjArray.
29
30 Revision 1.52 2001/05/16 14:57:20 alibrary
31 New files for folders and Stack
32
33 Revision 1.51 2001/05/14 10:18:55 hristov
34 Default arguments declared once
35
36 Revision 1.50 2001/05/10 14:44:16 jbarbosa
37 Corrected some overlaps (thanks I. Hrivnacovna).
38
39 Revision 1.49 2001/05/10 12:23:49 jbarbosa
40 Repositioned the RICH modules.
41 Eliminated magic numbers.
42 Incorporated diagnostics (from macros).
43
44 Revision 1.48 2001/03/15 10:35:00 jbarbosa
45 Corrected bug in MakeBranch (was using a different version of STEER)
46
47 Revision 1.47 2001/03/14 18:13:56 jbarbosa
48 Several changes to adapt to new IO.
49 Removed digitising function, using AliRICHMerger::Digitise from now on.
50
51 Revision 1.46 2001/03/12 17:46:33 hristov
52 Changes needed on Sun with CC 5.0
53
54 Revision 1.45 2001/02/27 22:11:46 jbarbosa
55 Testing TreeS, removing of output.
56
57 Revision 1.44 2001/02/27 15:19:12 jbarbosa
58 Transition to SDigits.
59
60 Revision 1.43 2001/02/23 17:19:06 jbarbosa
61 Corrected photocathode definition in BuildGeometry().
62
63 Revision 1.42 2001/02/13 20:07:23 jbarbosa
64 Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum
65 when entering the freon. Corrected calls to particle stack.
66
67 Revision 1.41 2001/01/26 20:00:20 hristov
68 Major upgrade of AliRoot code
69
70 Revision 1.40 2001/01/24 20:58:03 jbarbosa
71 Enhanced BuildGeometry. Now the photocathodes are drawn.
72
73 Revision 1.39 2001/01/22 21:40:24 jbarbosa
74 Removing magic numbers
75
76 Revision 1.37 2000/12/20 14:07:25 jbarbosa
77 Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
78
79 Revision 1.36 2000/12/18 17:45:54 jbarbosa
80 Cleaned up PadHits object.
81
82 Revision 1.35 2000/12/15 16:49:40 jbarbosa
83 Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
84
85 Revision 1.34 2000/11/10 18:12:12 jbarbosa
86 Bug fix for AliRICHCerenkov (thanks to P. Hristov)
87
88 Revision 1.33 2000/11/02 10:09:01 jbarbosa
89 Minor bug correction (some pointers were not initialised in the default constructor)
90
91 Revision 1.32 2000/11/01 15:32:55 jbarbosa
92 Updated to handle both reconstruction algorithms.
93
94 Revision 1.31 2000/10/26 20:18:33 jbarbosa
95 Supports for methane and freon vessels
96
97 Revision 1.30 2000/10/24 13:19:12 jbarbosa
98 Geometry updates.
99
100 Revision 1.29 2000/10/19 19:39:25 jbarbosa
101 Some more changes to geometry. Further correction of digitisation "per part. type"
102
103 Revision 1.28 2000/10/17 20:50:57 jbarbosa
104 Inversed digtise by particle type (now, only the selected particle type is not digitsed).
105 Corrected several geometry minor bugs.
106 Added new parameter (opaque quartz thickness).
107
108 Revision 1.27 2000/10/11 10:33:55 jbarbosa
109 Corrected bug introduced by earlier revisions (CerenkovData array cannot be reset to zero on wach call of StepManager)
110
111 Revision 1.26 2000/10/03 21:44:08 morsch
112 Use AliSegmentation and AliHit abstract base classes.
113
114 Revision 1.25 2000/10/02 21:28:12 fca
115 Removal of useless dependecies via forward declarations
116
117 Revision 1.24 2000/10/02 15:43:17 jbarbosa
118 Fixed forward declarations.
119 Fixed honeycomb density.
120 Fixed cerenkov storing.
121 New electronics.
122
123 Revision 1.23 2000/09/13 10:42:14 hristov
124 Minor corrections for HP, DEC and Sun; strings.h included
125
126 Revision 1.22 2000/09/12 18:11:13 fca
127 zero hits area before using
128
129 Revision 1.21 2000/07/21 10:21:07 morsch
130 fNrawch = 0; and fNrechits = 0; in the default constructor.
131
132 Revision 1.20 2000/07/10 15:28:39 fca
133 Correction of the inheritance scheme
134
135 Revision 1.19 2000/06/30 16:29:51 dibari
136 Added kDebugLevel variable to control output size on demand
137
138 Revision 1.18 2000/06/12 15:15:46 jbarbosa
139 Cleaned up version.
140
141 Revision 1.17 2000/06/09 14:58:37 jbarbosa
142 New digitisation per particle type
143
144 Revision 1.16 2000/04/19 12:55:43 morsch
145 Newly structured and updated version (JB, AM)
146
147*/
148
149
150////////////////////////////////////////////////
151// Manager and hits classes for set:RICH //
152////////////////////////////////////////////////
153
154#include <TBRIK.h>
155#include <TTUBE.h>
156#include <TNode.h>
157#include <TRandom.h>
158#include <TObject.h>
159#include <TVector.h>
160#include <TObjArray.h>
161#include <TArrayF.h>
162#include <TFile.h>
163#include <TParticle.h>
164#include <TGeometry.h>
165#include <TTree.h>
166#include <TH1.h>
167#include <TH2.h>
168#include <TCanvas.h>
169//#include <TPad.h>
170#include <TF1.h>
171
172#include <iostream.h>
173#include <strings.h>
174
175#include "AliRICH.h"
176#include "AliSegmentation.h"
177#include "AliRICHSegmentationV0.h"
178#include "AliRICHHit.h"
179#include "AliRICHCerenkov.h"
180#include "AliRICHSDigit.h"
181#include "AliRICHDigit.h"
182#include "AliRICHTransientDigit.h"
183#include "AliRICHRawCluster.h"
184#include "AliRICHRecHit1D.h"
185#include "AliRICHRecHit3D.h"
186#include "AliRICHHitMapA1.h"
187#include "AliRICHClusterFinder.h"
188#include "AliRICHMerger.h"
189#include "AliRun.h"
190#include "AliMC.h"
191#include "AliMagF.h"
192#include "AliConst.h"
193#include "AliPDG.h"
194#include "AliPoints.h"
195#include "AliCallf77.h"
196
197
198// Static variables for the pad-hit iterator routines
199static Int_t sMaxIterPad=0;
200static Int_t sCurIterPad=0;
201
202ClassImp(AliRICH)
203
204//___________________________________________
205AliRICH::AliRICH()
206{
207// Default constructor for RICH manager class
208
209 fIshunt = 0;
210 fHits = 0;
211 fSDigits = 0;
212 fNSDigits = 0;
213 fNcerenkovs = 0;
214 fDchambers = 0;
215 fRecHits1D = 0;
216 fRecHits3D = 0;
217 fRawClusters = 0;
218 fChambers = 0;
219 fCerenkovs = 0;
220 for (Int_t i=0; i<7; i++)
221 {
222 fNdch[i] = 0;
223 fNrawch[i] = 0;
224 fNrechits1D[i] = 0;
225 fNrechits3D[i] = 0;
226 }
227
228 fFileName = 0;
229 fMerger = 0;
230}
231
232//___________________________________________
233AliRICH::AliRICH(const char *name, const char *title)
234 : AliDetector(name,title)
235{
236//Begin_Html
237/*
238 <img src="gif/alirich.gif">
239*/
240//End_Html
241
242 fHits = new TClonesArray("AliRICHHit",1000 );
243 gAlice->AddHitList(fHits);
244 fSDigits = new TClonesArray("AliRICHSDigit",100000);
245 fCerenkovs = new TClonesArray("AliRICHCerenkov",1000);
246 gAlice->AddHitList(fCerenkovs);
247 //gAlice->AddHitList(fHits);
248 fNSDigits = 0;
249 fNcerenkovs = 0;
250 fIshunt = 0;
251
252 //fNdch = new Int_t[kNCH];
253
254 fDchambers = new TObjArray(kNCH);
255
256 fRecHits1D = new TObjArray(kNCH);
257 fRecHits3D = new TObjArray(kNCH);
258
259 Int_t i;
260
261 for (i=0; i<kNCH ;i++) {
262 //PH (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
263 fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i);
264 fNdch[i]=0;
265 }
266
267 //fNrawch = new Int_t[kNCH];
268
269 fRawClusters = new TObjArray(kNCH);
270 //printf("Created fRwClusters with adress:%p",fRawClusters);
271
272 for (i=0; i<kNCH ;i++) {
273 //PH (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
274 fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i);
275 fNrawch[i]=0;
276 }
277
278 //fNrechits = new Int_t[kNCH];
279
280 for (i=0; i<kNCH ;i++) {
281 //PH (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
282 fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
283 }
284 for (i=0; i<kNCH ;i++) {
285 //PH (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
286 fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
287 }
288 //printf("Created fRecHits with adress:%p",fRecHits);
289
290
291 SetMarkerColor(kRed);
292
293 /*fChambers = new TObjArray(kNCH);
294 for (i=0; i<kNCH; i++)
295 (*fChambers)[i] = new AliRICHChamber();*/
296
297 fFileName = 0;
298 fMerger = 0;
299}
300
301AliRICH::AliRICH(const AliRICH& RICH)
302{
303// Copy Constructor
304}
305
306
307//___________________________________________
308AliRICH::~AliRICH()
309{
310
311// Destructor of RICH manager class
312
313 fIshunt = 0;
314 delete fHits;
315 delete fSDigits;
316 delete fCerenkovs;
317
318 //PH Delete TObjArrays
319 if (fChambers) {
320 fChambers->Delete();
321 delete fChambers;
322 }
323 if (fDchambers) {
324 fDchambers->Delete();
325 delete fDchambers;
326 }
327 if (fRawClusters) {
328 fRawClusters->Delete();
329 delete fRawClusters;
330 }
331 if (fRecHits1D) {
332 fRecHits1D->Delete();
333 delete fRecHits1D;
334 }
335 if (fRecHits3D) {
336 fRecHits3D->Delete();
337 delete fRecHits3D;
338 }
339
340}
341
342
343//_____________________________________________________________________________
344Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
345{
346//
347// Calls the charge disintegration method of the current chamber and adds
348// the simulated cluster to the root treee
349//
350 Int_t clhits[5];
351 Float_t newclust[4][500];
352 Int_t nnew;
353
354//
355// Integrated pulse height on chamber
356
357 clhits[0]=fNhits+1;
358
359 //PH ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
360 ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
361 Int_t ic=0;
362
363//
364// Add new clusters
365 for (Int_t i=0; i<nnew; i++) {
366 if (Int_t(newclust[0][i]) > 0) {
367 ic++;
368// Cluster Charge
369 clhits[1] = Int_t(newclust[0][i]);
370// Pad: ix
371 clhits[2] = Int_t(newclust[1][i]);
372// Pad: iy
373 clhits[3] = Int_t(newclust[2][i]);
374// Pad: chamber sector
375 clhits[4] = Int_t(newclust[3][i]);
376
377 //printf(" %d %d %d %d %d\n", clhits[0], clhits[1], clhits[2], clhits[3], clhits[4]);
378
379 AddSDigit(clhits);
380 }
381 }
382
383 if (gAlice->TreeS())
384 {
385 gAlice->TreeS()->Fill();
386 gAlice->TreeS()->Write(0,TObject::kOverwrite);
387 //printf("Filled SDigits...\n");
388 }
389
390return nnew;
391}
392//___________________________________________
393void AliRICH::Hits2SDigits()
394{
395
396// Dummy: sdigits are created during transport.
397// Called from alirun.
398
399 int nparticles = gAlice->GetNtrack();
400 cout << "Particles (RICH):" <<nparticles<<endl;
401 if (nparticles > 0) printf("SDigits were already generated.\n");
402
403}
404
405//___________________________________________
406void AliRICH::SDigits2Digits(Int_t nev, Int_t flag)
407{
408
409//
410// Generate digits.
411// Called from macro. Multiple events, more functionality.
412
413 AliRICHChamber* iChamber;
414
415 printf("Generating tresholds...\n");
416
417 for(Int_t i=0;i<7;i++) {
418 iChamber = &(Chamber(i));
419 iChamber->GenerateTresholds();
420 }
421
422 int nparticles = gAlice->GetNtrack();
423 cout << "Particles (RICH):" <<nparticles<<endl;
424 if (nparticles <= 0) return;
425 if (!fMerger) {
426 fMerger = new AliRICHMerger();
427 }
428 fMerger->Init();
429 fMerger->Digitise(nev,flag);
430}
431//___________________________________________
432void AliRICH::SDigits2Digits()
433{
434 SDigits2Digits(0,0);
435}
436//___________________________________________
437void AliRICH::Digits2Reco()
438{
439
440// Generate clusters
441// Called from alirun, single event only.
442
443 int nparticles = gAlice->GetNtrack();
444 cout << "Particles (RICH):" <<nparticles<<endl;
445 if (nparticles > 0) FindClusters(0,0);
446
447}
448
449//___________________________________________
450void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
451{
452
453//
454// Adds a hit to the Hits list
455
456 TClonesArray &lhits = *fHits;
457 new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
458}
459//_____________________________________________________________________________
460void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
461{
462
463//
464// Adds a RICH cerenkov hit to the Cerenkov Hits list
465//
466
467 TClonesArray &lcerenkovs = *fCerenkovs;
468 new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
469 //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
470}
471//___________________________________________
472void AliRICH::AddSDigit(Int_t *clhits)
473{
474
475//
476// Add a RICH pad hit to the list
477//
478
479 //printf("fsdigits:%p, data: %d\n",fSDigits,clhits[2]);
480 TClonesArray &lSDigits = *fSDigits;
481 new(lSDigits[fNSDigits++]) AliRICHSDigit(clhits);
482}
483//_____________________________________________________________________________
484void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
485{
486
487 //
488 // Add a RICH digit to the list
489 //
490
491 //printf("fdigits:%p, data: %d\n",((TClonesArray*)(*fDchambers)[id]),digits[0]);
492 //PH TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
493 TClonesArray &ldigits = *((TClonesArray*)fDchambers->At(id));
494 new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
495}
496
497//_____________________________________________________________________________
498void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
499{
500 //
501 // Add a RICH digit to the list
502 //
503
504 //PH TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
505 TClonesArray &lrawcl = *((TClonesArray*)fRawClusters->At(id));
506 new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
507}
508
509//_____________________________________________________________________________
510void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
511{
512
513 //
514 // Add a RICH reconstructed hit to the list
515 //
516
517 //PH TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
518 TClonesArray &lrec1D = *((TClonesArray*)fRecHits1D->At(id));
519 new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
520}
521
522//_____________________________________________________________________________
523void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
524{
525
526 //
527 // Add a RICH reconstructed hit to the list
528 //
529
530 //PH TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
531 TClonesArray &lrec3D = *((TClonesArray*)fRecHits3D->At(id));
532 new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
533}
534
535//___________________________________________
536void AliRICH::BuildGeometry()
537
538{
539
540 //
541 // Builds a TNode geometry for event display
542 //
543 TNode *node, *subnode, *top;
544
545 const int kColorRICH = kRed;
546 //
547 top=gAlice->GetGeometry()->GetNode("alice");
548
549 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
550 AliRICHSegmentationV0* segmentation;
551 AliRICHChamber* iChamber;
552 AliRICHGeometry* geometry;
553
554 iChamber = &(pRICH->Chamber(0));
555 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
556 geometry=iChamber->GetGeometryModel();
557
558 new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
559
560 Float_t padplane_width = segmentation->GetPadPlaneWidth();
561 Float_t padplane_length = segmentation->GetPadPlaneLength();
562
563 //printf("\n\n\n\n\n In BuildGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f\n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy());
564
565 new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2);
566
567 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2);
568 //printf("\n\n\n\n\n Padplane w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength());
569
570 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
571 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
572 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
573 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
574 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
575 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
576 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
577
578 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
579
580 new TRotMatrix("rot993","rot993",90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
581 new TRotMatrix("rot994","rot994",90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
582 new TRotMatrix("rot995","rot995",90., 0. , 90. , 90. , 0. , 0. );
583 new TRotMatrix("rot996","rot996",90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
584 new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
585 new TRotMatrix("rot998","rot998",90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
586 new TRotMatrix("rot999","rot999",90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
587
588 Float_t pos1[3]={0. , offset*cosphi , offset*sinphi};
589 Float_t pos2[3]={offset*sintheta , offset*costheta , 0. };
590 Float_t pos3[3]={0. , offset , 0.};
591 Float_t pos4[3]={-offset*sintheta , offset*costheta , 0.};
592 Float_t pos5[3]={offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
593 Float_t pos6[3]={0. , offset*cosphi , -offset*sinphi};
594 Float_t pos7[3]={ -offset*sinphi , offset*costheta*cosphi, -offset*sinphi};
595
596
597 top->cd();
598 //Float_t pos1[3]={0,471.8999,165.2599};
599 //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
600 //new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
601 node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
602 node->SetLineColor(kColorRICH);
603 node->cd();
604 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
605 subnode->SetLineColor(kGreen);
606 fNodes->Add(subnode);
607 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
608 subnode->SetLineColor(kGreen);
609 fNodes->Add(subnode);
610 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
611 subnode->SetLineColor(kGreen);
612 fNodes->Add(subnode);
613 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
614 subnode->SetLineColor(kGreen);
615 fNodes->Add(subnode);
616 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
617 subnode->SetLineColor(kGreen);
618 fNodes->Add(subnode);
619 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
620 subnode->SetLineColor(kGreen);
621 fNodes->Add(subnode);
622 fNodes->Add(node);
623
624
625 top->cd();
626 //Float_t pos2[3]={171,470,0};
627 //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
628 //new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
629 node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
630 node->SetLineColor(kColorRICH);
631 node->cd();
632 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
633 subnode->SetLineColor(kGreen);
634 fNodes->Add(subnode);
635 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
636 subnode->SetLineColor(kGreen);
637 fNodes->Add(subnode);
638 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
639 subnode->SetLineColor(kGreen);
640 fNodes->Add(subnode);
641 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
642 subnode->SetLineColor(kGreen);
643 fNodes->Add(subnode);
644 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
645 subnode->SetLineColor(kGreen);
646 fNodes->Add(subnode);
647 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
648 subnode->SetLineColor(kGreen);
649 fNodes->Add(subnode);
650 fNodes->Add(node);
651
652
653 top->cd();
654 //Float_t pos3[3]={0,500,0};
655 //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
656 //new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
657 node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
658 node->SetLineColor(kColorRICH);
659 node->cd();
660 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
661 subnode->SetLineColor(kGreen);
662 fNodes->Add(subnode);
663 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
664 subnode->SetLineColor(kGreen);
665 fNodes->Add(subnode);
666 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
667 subnode->SetLineColor(kGreen);
668 fNodes->Add(subnode);
669 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
670 subnode->SetLineColor(kGreen);
671 fNodes->Add(subnode);
672 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
673 subnode->SetLineColor(kGreen);
674 fNodes->Add(subnode);
675 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
676 subnode->SetLineColor(kGreen);
677 fNodes->Add(subnode);
678 fNodes->Add(node);
679
680 top->cd();
681 //Float_t pos4[3]={-171,470,0};
682 //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
683 //new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
684 node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
685 node->SetLineColor(kColorRICH);
686 node->cd();
687 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
688 subnode->SetLineColor(kGreen);
689 fNodes->Add(subnode);
690 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
691 subnode->SetLineColor(kGreen);
692 fNodes->Add(subnode);
693 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
694 subnode->SetLineColor(kGreen);
695 fNodes->Add(subnode);
696 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
697 subnode->SetLineColor(kGreen);
698 fNodes->Add(subnode);
699 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
700 subnode->SetLineColor(kGreen);
701 fNodes->Add(subnode);
702 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
703 subnode->SetLineColor(kGreen);
704 fNodes->Add(subnode);
705 fNodes->Add(node);
706
707
708 top->cd();
709 //Float_t pos5[3]={161.3999,443.3999,-165.3};
710 //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
711 //new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
712 node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
713 node->SetLineColor(kColorRICH);
714 node->cd();
715 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
716 subnode->SetLineColor(kGreen);
717 fNodes->Add(subnode);
718 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
719 subnode->SetLineColor(kGreen);
720 fNodes->Add(subnode);
721 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
722 subnode->SetLineColor(kGreen);
723 fNodes->Add(subnode);
724 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
725 subnode->SetLineColor(kGreen);
726 fNodes->Add(subnode);
727 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
728 subnode->SetLineColor(kGreen);
729 fNodes->Add(subnode);
730 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
731 subnode->SetLineColor(kGreen);
732 fNodes->Add(subnode);
733 fNodes->Add(node);
734
735
736 top->cd();
737 //Float_t pos6[3]={0., 471.9, -165.3,};
738 //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
739 //new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
740 node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
741 node->SetLineColor(kColorRICH);
742 fNodes->Add(node);node->cd();
743 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
744 subnode->SetLineColor(kGreen);
745 fNodes->Add(subnode);
746 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
747 subnode->SetLineColor(kGreen);
748 fNodes->Add(subnode);
749 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
750 subnode->SetLineColor(kGreen);
751 fNodes->Add(subnode);
752 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
753 subnode->SetLineColor(kGreen);
754 fNodes->Add(subnode);
755 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
756 subnode->SetLineColor(kGreen);
757 fNodes->Add(subnode);
758 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
759 subnode->SetLineColor(kGreen);
760 fNodes->Add(subnode);
761
762
763 top->cd();
764 //Float_t pos7[3]={-161.399,443.3999,-165.3};
765 //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
766 //new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
767 node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
768 node->SetLineColor(kColorRICH);
769 node->cd();
770 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
771 subnode->SetLineColor(kGreen);
772 fNodes->Add(subnode);
773 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
774 subnode->SetLineColor(kGreen);
775 fNodes->Add(subnode);
776 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
777 subnode->SetLineColor(kGreen);
778 fNodes->Add(subnode);
779 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
780 subnode->SetLineColor(kGreen);
781 fNodes->Add(subnode);
782 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
783 subnode->SetLineColor(kGreen);
784 fNodes->Add(subnode);
785 subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
786 subnode->SetLineColor(kGreen);
787 fNodes->Add(subnode);
788 fNodes->Add(node);
789
790}
791
792//___________________________________________
793void AliRICH::CreateGeometry()
794{
795 //
796 // Create the geometry for RICH version 1
797 //
798 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
799 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
800 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
801 //
802 //Begin_Html
803 /*
804 <img src="picts/AliRICHv1.gif">
805 */
806 //End_Html
807 //Begin_Html
808 /*
809 <img src="picts/AliRICHv1Tree.gif">
810 */
811 //End_Html
812
813 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
814 AliRICHSegmentationV0* segmentation;
815 AliRICHGeometry* geometry;
816 AliRICHChamber* iChamber;
817
818 iChamber = &(pRICH->Chamber(0));
819 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
820 geometry=iChamber->GetGeometryModel();
821
822 Float_t distance;
823 distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
824 geometry->SetRadiatorToPads(distance);
825
826 //Opaque quartz thickness
827 Float_t oqua_thickness = .5;
828 //CsI dimensions
829
830 //Float_t csi_length = 160*.8 + 2.6;
831 //Float_t csi_width = 144*.84 + 2*2.6;
832
833 Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
834 Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
835
836 //printf("\n\n\n\n\n In CreateGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f deadzone: %f \n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy(),segmentation->DeadZone());
837
838 Int_t *idtmed = fIdtmed->GetArray()-999;
839
840 Int_t i;
841 Float_t zs;
842 Int_t idrotm[1099];
843 Float_t par[3];
844
845 // --- Define the RICH detector
846 // External aluminium box
847 par[0] = 68.8;
848 par[1] = 13; //Original Settings
849 par[2] = 70.86;
850 /*par[0] = 73.15;
851 par[1] = 11.5;
852 par[2] = 71.1;*/
853 gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
854
855 // Air
856 par[0] = 66.3;
857 par[1] = 13; //Original Settings
858 par[2] = 68.35;
859 /*par[0] = 66.55;
860 par[1] = 11.5;
861 par[2] = 64.8;*/
862 gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
863
864 // Air 2 (cutting the lower part of the box)
865
866 par[0] = 1.25;
867 par[1] = 3; //Original Settings
868 par[2] = 70.86;
869 gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
870
871 // Air 3 (cutting the lower part of the box)
872
873 par[0] = 66.3;
874 par[1] = 3; //Original Settings
875 par[2] = 1.2505;
876 gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
877
878 // Honeycomb
879 par[0] = 66.3;
880 par[1] = .188; //Original Settings
881 par[2] = 68.35;
882 /*par[0] = 66.55;
883 par[1] = .188;
884 par[2] = 63.1;*/
885 gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
886
887 // Aluminium sheet
888 par[0] = 66.3;
889 par[1] = .025; //Original Settings
890 par[2] = 68.35;
891 /*par[0] = 66.5;
892 par[1] = .025;
893 par[2] = 63.1;*/
894 gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
895
896 // Quartz
897 par[0] = geometry->GetQuartzWidth()/2;
898 par[1] = geometry->GetQuartzThickness()/2;
899 par[2] = geometry->GetQuartzLength()/2;
900 /*par[0] = 63.1;
901 par[1] = .25; //Original Settings
902 par[2] = 65.5;*/
903 /*par[0] = geometry->GetQuartzWidth()/2;
904 par[1] = geometry->GetQuartzThickness()/2;
905 par[2] = geometry->GetQuartzLength()/2;*/
906 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
907 gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
908
909 // Spacers (cylinders)
910 par[0] = 0.;
911 par[1] = .5;
912 par[2] = geometry->GetFreonThickness()/2;
913 gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
914
915 // Feet (freon slabs supports)
916
917 par[0] = .7;
918 par[1] = .3;
919 par[2] = 1.9;
920 gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
921
922 // Opaque quartz
923 par[0] = geometry->GetQuartzWidth()/2;
924 par[1] = .2;
925 par[2] = geometry->GetQuartzLength()/2;
926 /*par[0] = 61.95;
927 par[1] = .2; //Original Settings
928 par[2] = 66.5;*/
929 /*par[0] = 66.5;
930 par[1] = .2;
931 par[2] = 61.95;*/
932 gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
933
934 // Frame of opaque quartz
935 par[0] = geometry->GetOuterFreonWidth()/2;
936 //+ oqua_thickness;
937 par[1] = geometry->GetFreonThickness()/2;
938 par[2] = geometry->GetOuterFreonLength()/2;
939 //+ oqua_thickness;
940 /*par[0] = 20.65;
941 par[1] = .5; //Original Settings
942 par[2] = 66.5;*/
943 /*par[0] = 66.5;
944 par[1] = .5;
945 par[2] = 20.65;*/
946 gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
947
948 par[0] = geometry->GetInnerFreonWidth()/2;
949 par[1] = geometry->GetFreonThickness()/2;
950 par[2] = geometry->GetInnerFreonLength()/2;
951 gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
952
953 // Little bar of opaque quartz
954 //par[0] = .275;
955 //par[1] = geometry->GetQuartzThickness()/2;
956 //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
957 //par[2] = geometry->GetInnerFreonLength()/2;
958 //+ oqua_thickness;
959 /*par[0] = .275;
960 par[1] = .25; //Original Settings
961 par[2] = 63.1;*/
962 /*par[0] = 63.1;
963 par[1] = .25;
964 par[2] = .275;*/
965 //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
966
967 // Freon
968 par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
969 par[1] = geometry->GetFreonThickness()/2;
970 par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
971 /*par[0] = 20.15;
972 par[1] = .5; //Original Settings
973 par[2] = 65.5;*/
974 /*par[0] = 65.5;
975 par[1] = .5;
976 par[2] = 20.15;*/
977 gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
978
979 par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
980 par[1] = geometry->GetFreonThickness()/2;
981 par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
982 gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
983
984 // Methane
985 //par[0] = 64.8;
986 par[0] = csi_width/2;
987 par[1] = geometry->GetGapThickness()/2;
988 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
989 //par[2] = 64.8;
990 par[2] = csi_length/2;
991 gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
992
993 // Methane gap
994 //par[0] = 64.8;
995 par[0] = csi_width/2;
996 par[1] = geometry->GetProximityGapThickness()/2;
997 //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
998 //par[2] = 64.8;
999 par[2] = csi_length/2;
1000 gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1001
1002 // CsI photocathode
1003 //par[0] = 64.8;
1004 par[0] = csi_width/2;
1005 par[1] = .25;
1006 //par[2] = 64.8;
1007 par[2] = csi_length/2;
1008 gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1009
1010 // Anode grid
1011 par[0] = 0.;
1012 par[1] = .001;
1013 par[2] = 20.;
1014 gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1015
1016 // Wire supports
1017 // Bar of metal
1018
1019 par[0] = csi_width/2;
1020 par[1] = 1.05;
1021 par[2] = 1.05;
1022 gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
1023
1024 // Ceramic pick up (base)
1025
1026 par[0] = csi_width/2;
1027 par[1] = .25;
1028 par[2] = 1.05;
1029 gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
1030
1031 // Ceramic pick up (head)
1032
1033 par[0] = csi_width/2;
1034 par[1] = .1;
1035 par[2] = .1;
1036 gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
1037
1038 // Aluminium supports for methane and CsI
1039 // Short bar
1040
1041 par[0] = csi_width/2;
1042 par[1] = geometry->GetGapThickness()/2 + .25;
1043 par[2] = (68.35 - csi_length/2)/2;
1044 gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
1045
1046 // Long bar
1047
1048 par[0] = (66.3 - csi_width/2)/2;
1049 par[1] = geometry->GetGapThickness()/2 + .25;
1050 par[2] = csi_length/2 + 68.35 - csi_length/2;
1051 gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
1052
1053 // Aluminium supports for freon
1054 // Short bar
1055
1056 par[0] = geometry->GetQuartzWidth()/2;
1057 par[1] = .3;
1058 par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
1059 gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
1060
1061 // Long bar
1062
1063 par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
1064 par[1] = .3;
1065 par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
1066 gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
1067
1068 // PCB backplane
1069
1070 par[0] = csi_width/2;
1071 par[1] = .25;
1072 par[2] = csi_length/4 -.5025;
1073 gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
1074
1075
1076 // Backplane supports
1077
1078 // Aluminium slab
1079
1080 par[0] = 33.15;
1081 par[1] = 2;
1082 par[2] = 21.65;
1083 gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
1084
1085 // Big hole
1086
1087 par[0] = 9.05;
1088 par[1] = 2;
1089 par[2] = 4.4625;
1090 gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
1091
1092 // Small hole
1093
1094 par[0] = 5.7;
1095 par[1] = 2;
1096 par[2] = 4.4625;
1097 gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
1098
1099 // Place holes inside backplane support
1100
1101 gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
1102 gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
1103 gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
1104 gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
1105 gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1106 gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1107 gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1108 gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1109 gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
1110 gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
1111 gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1112 gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1113 gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1114 gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1115 gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1116 gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1117
1118
1119
1120 // --- Places the detectors defined with GSVOLU
1121 // Place material inside RICH
1122 gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
1123 gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1124 gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1125 gMC->Gspos("AIR3", 1, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
1126 gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 68.35 + 1.25, 0, "ONLY");
1127
1128
1129 gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
1130 gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
1131 gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
1132 gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1133 gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
1134 gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1135 gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1136 gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1137 gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1138 gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1139 gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1140 gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
1141
1142 // Supports placing
1143
1144 // Methane supports
1145 gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1146 gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1147 gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
1148 gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
1149
1150 //Freon supports
1151
1152 Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
1153
1154 gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1155 gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1156 gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1157 gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1158
1159 AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1160
1161 //Placing of the spacers inside the freon slabs
1162
1163 Int_t nspacers = 30;
1164 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers);
1165
1166 //printf("Nspacers: %d", nspacers);
1167
1168 for (i = 0; i < nspacers/3; i++) {
1169 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1170 gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1171 }
1172
1173 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1174 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1175 gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
1176 }
1177
1178 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1179 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1180 gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1181 }
1182
1183 for (i = 0; i < nspacers/3; i++) {
1184 zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1185 gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1186 }
1187
1188 for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1189 zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1190 gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY"); //Original settings
1191 }
1192
1193 for (i = (nspacers*2)/3; i < nspacers; ++i) {
1194 zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1195 gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1196 }
1197
1198
1199 gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
1200 gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
1201 gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
1202// printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2);
1203 gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings
1204 gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (-31.3)
1205 //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (-21.65)
1206 //gMC->Gspos("BARR", 2, "QUAR", geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY"); //Original settings (21.65)
1207 gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
1208 gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
1209 gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1210 gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
1211 printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25);
1212
1213 // Wire support placing
1214
1215 gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
1216 gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
1217 gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
1218
1219 // Backplane placing
1220
1221 gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
1222 gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
1223 gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1224 gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1225 gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1226 gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1227
1228 // PCB placing
1229
1230 gMC->Gspos("PCB ", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY");
1231 gMC->Gspos("PCB ", 2, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
1232
1233
1234
1235 //printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2);
1236
1237 // Place RICH inside ALICE apparatus
1238
1239 /* old values
1240
1241 AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1242 AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1243 AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1244 AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1245 AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1246 AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1247 AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1248
1249 gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26, idrotm[1000], "ONLY");
1250 gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0., idrotm[1001], "ONLY");
1251 gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0., idrotm[1002], "ONLY");
1252 gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0., idrotm[1003], "ONLY");
1253 gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3, idrotm[1004], "ONLY");
1254 gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3, idrotm[1005], "ONLY");
1255 gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");*/
1256
1257 // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm)
1258
1259 Float_t offset = 490 + 1.276 - geometry->GetGapThickness()/2; //distance from center of mother volume to methane
1260 Float_t deltaphi = 19.5; //phi angle between center of chambers - z direction
1261 Float_t deltatheta = 20; //theta angle between center of chambers - x direction
1262 Float_t cosphi = TMath::Cos(deltaphi*TMath::Pi()/180);
1263 Float_t sinphi = TMath::Sin(deltaphi*TMath::Pi()/180);
1264 Float_t costheta = TMath::Cos(deltatheta*TMath::Pi()/180);
1265 Float_t sintheta = TMath::Sin(deltatheta*TMath::Pi()/180);
1266
1267 //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
1268
1269 AliMatrix(idrotm[1000], 90., 0. , 90. - deltaphi, 90. , deltaphi, -90. );
1270 AliMatrix(idrotm[1001], 90., -deltatheta , 90. , 90.- deltatheta , 0. , 0. );
1271 AliMatrix(idrotm[1002], 90., 0. , 90. , 90. , 0. , 0. );
1272 AliMatrix(idrotm[1003], 90., deltatheta , 90. , 90 + deltatheta , 0. , 0. );
1273 AliMatrix(idrotm[1004], 90., 360. - deltatheta, 108.2 , 90.- deltatheta ,18.2 , 90 - deltatheta);
1274 AliMatrix(idrotm[1005], 90., 0. , 90 + deltaphi , 90. , deltaphi, 90. );
1275 AliMatrix(idrotm[1006], 90., deltatheta , 108.2 , 90.+ deltatheta ,18.2 , 90 + deltatheta);
1276
1277 gMC->Gspos("RICH", 1, "ALIC", 0. , offset*cosphi , offset*sinphi ,idrotm[1000], "ONLY");
1278 gMC->Gspos("RICH", 2, "ALIC", (offset)*sintheta , offset*costheta , 0. ,idrotm[1001], "ONLY");
1279 gMC->Gspos("RICH", 3, "ALIC", 0. , offset , 0. ,idrotm[1002], "ONLY");
1280 gMC->Gspos("RICH", 4, "ALIC", -(offset)*sintheta, offset*costheta , 0. ,idrotm[1003], "ONLY");
1281 gMC->Gspos("RICH", 5, "ALIC", (offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1004], "ONLY");
1282 gMC->Gspos("RICH", 6, "ALIC", 0. , offset*cosphi , -offset*sinphi,idrotm[1005], "ONLY");
1283 gMC->Gspos("RICH", 7, "ALIC", -(offset)*sinphi , offset*costheta*cosphi, -offset*sinphi,idrotm[1006], "ONLY");
1284
1285}
1286
1287
1288//___________________________________________
1289void AliRICH::CreateMaterials()
1290{
1291 //
1292 // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1293 // ORIGIN : NICK VAN EIJNDHOVEN
1294 // Modified by: N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1295 // R.A. Fini (INFN - BARI, Rosanna.Fini@ba.infn.it)
1296 // R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1297 //
1298 Int_t isxfld = gAlice->Field()->Integ();
1299 Float_t sxmgmx = gAlice->Field()->Max();
1300 Int_t i;
1301
1302 /************************************Antonnelo's Values (14-vectors)*****************************************/
1303 /*
1304 Float_t ppckov[14] = { 5.63e-9,5.77e-9,5.9e-9,6.05e-9,6.2e-9,6.36e-9,6.52e-9,
1305 6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1306 Float_t rIndexQuarz[14] = { 1.528309,1.533333,
1307 1.538243,1.544223,1.550568,1.55777,
1308 1.565463,1.574765,1.584831,1.597027,
1309 1.611858,1.6277,1.6472,1.6724 };
1310 Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1311 Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1312 Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1313 Float_t abscoFreon[14] = { 179.0987,179.0987,
1314 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
1315 //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1316 // 1e-5,1e-5,1e-5,1e-5,1e-5 };
1317 Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
1318 14.177,9.282,4.0925,1.149,.3627,.10857 };
1319 Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1320 1e-5,1e-5,1e-5,1e-5,1e-5 };
1321 Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1322 1e-4,1e-4,1e-4,1e-4 };
1323 Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1324 1e6,1e6,1e6 };
1325 Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1326 1e-4,1e-4,1e-4,1e-4 };
1327 Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1328 Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
1329 .17425,.1785,.1836,.1904,.1938,.221 };
1330 Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1331 */
1332
1333
1334 /**********************************End of Antonnelo's Values**********************************/
1335
1336 /**********************************Values from rich_media.f (31-vectors)**********************************/
1337
1338
1339 //Photons energy intervals
1340 Float_t ppckov[26];
1341 for (i=0;i<26;i++)
1342 {
1343 ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
1344 //printf ("Energy intervals: %e\n",ppckov[i]);
1345 }
1346
1347
1348 //Refraction index for quarz
1349 Float_t rIndexQuarz[26];
1350 Float_t e1= 10.666;
1351 Float_t e2= 18.125;
1352 Float_t f1= 46.411;
1353 Float_t f2= 228.71;
1354 for (i=0;i<26;i++)
1355 {
1356 Float_t ene=ppckov[i]*1e9;
1357 Float_t a=f1/(e1*e1 - ene*ene);
1358 Float_t b=f2/(e2*e2 - ene*ene);
1359 rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
1360 //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
1361 }
1362
1363 //Refraction index for opaque quarz, methane and grid
1364 Float_t rIndexOpaqueQuarz[26];
1365 Float_t rIndexMethane[26];
1366 Float_t rIndexGrid[26];
1367 for (i=0;i<26;i++)
1368 {
1369 rIndexOpaqueQuarz[i]=1;
1370 rIndexMethane[i]=1.000444;
1371 rIndexGrid[i]=1;
1372 //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
1373 }
1374
1375 //Absorption index for freon
1376 Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987, 179.0987,
1377 179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
1378 7.069031, 4.461292, 2.028366, 1.293013, .577267, .40746, .334964, 0., 0., 0.};
1379
1380 //Absorption index for quarz
1381 /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1382 .906,.907,.907,.907};
1383 Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
1384 215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1385 Float_t abscoQuarz[31];
1386 for (Int_t i=0;i<31;i++)
1387 {
1388 Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1389 if (Xlam <= 160) abscoQuarz[i] = 0;
1390 if (Xlam > 250) abscoQuarz[i] = 1;
1391 else
1392 {
1393 for (Int_t j=0;j<21;j++)
1394 {
1395 //printf ("Passed\n");
1396 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1397 {
1398 Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1399 Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1400 abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1401 }
1402 }
1403 }
1404 printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1405 }*/
1406
1407 /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1408 39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1409 17.94531, 11.88753, 5.99128, 3.83503, 2.36661, 1.53155, 1.30582, 1.08574, .8779708,
1410 .675275, 0., 0., 0.};
1411
1412 for (Int_t i=0;i<31;i++)
1413 {
1414 abscoQuarz[i] = abscoQuarz[i]/10;
1415 }*/
1416
1417 Float_t abscoQuarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
1418 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1419 .192, .1497, .10857};
1420
1421 //Absorption index for methane
1422 Float_t abscoMethane[26];
1423 for (i=0;i<26;i++)
1424 {
1425 abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1426 //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1427 }
1428
1429 //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1430 Float_t abscoOpaqueQuarz[26];
1431 Float_t abscoCsI[26];
1432 Float_t abscoGrid[26];
1433 Float_t efficAll[26];
1434 Float_t efficGrid[26];
1435 for (i=0;i<26;i++)
1436 {
1437 abscoOpaqueQuarz[i]=1e-5;
1438 abscoCsI[i]=1e-4;
1439 abscoGrid[i]=1e-4;
1440 efficAll[i]=1;
1441 efficGrid[i]=1;
1442 //printf ("All must be 1: %e, %e, %e, %e, %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1443 }
1444
1445 //Efficiency for csi
1446
1447 Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1448 0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1449 0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1450 0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1451
1452
1453
1454 //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1455 //UNPOLARIZED PHOTONS
1456
1457 for (i=0;i<26;i++)
1458 {
1459 efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1460 //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1461 }
1462
1463 /*******************************************End of rich_media.f***************************************/
1464
1465
1466
1467
1468
1469
1470 Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1471 zmet[2], zqua[2];
1472 Int_t nlmatfre;
1473 Float_t densquao;
1474 Int_t nlmatmet, nlmatqua;
1475 Float_t wmatquao[2], rIndexFreon[26];
1476 Float_t aquao[2], epsil, stmin, zquao[2];
1477 Int_t nlmatquao;
1478 Float_t radlal, densal, tmaxfd, deemax, stemax;
1479 Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1480
1481 Int_t *idtmed = fIdtmed->GetArray()-999;
1482
1483 // --- Photon energy (GeV)
1484 // --- Refraction indexes
1485 for (i = 0; i < 26; ++i) {
1486 rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1487 //rIndexFreon[i] = 1;
1488 //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1489 }
1490
1491 // --- Detection efficiencies (quantum efficiency for CsI)
1492 // --- Define parameters for honeycomb.
1493 // Used carbon of equivalent rad. lenght
1494
1495 ahon = 12.01;
1496 zhon = 6.;
1497 denshon = 0.1;
1498 radlhon = 18.8;
1499
1500 // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1501
1502 aqua[0] = 28.09;
1503 aqua[1] = 16.;
1504 zqua[0] = 14.;
1505 zqua[1] = 8.;
1506 densqua = 2.64;
1507 nlmatqua = -2;
1508 wmatqua[0] = 1.;
1509 wmatqua[1] = 2.;
1510
1511 // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1512
1513 aquao[0] = 28.09;
1514 aquao[1] = 16.;
1515 zquao[0] = 14.;
1516 zquao[1] = 8.;
1517 densquao = 2.64;
1518 nlmatquao = -2;
1519 wmatquao[0] = 1.;
1520 wmatquao[1] = 2.;
1521
1522 // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1523
1524 afre[0] = 12.;
1525 afre[1] = 19.;
1526 zfre[0] = 6.;
1527 zfre[1] = 9.;
1528 densfre = 1.7;
1529 nlmatfre = -2;
1530 wmatfre[0] = 6.;
1531 wmatfre[1] = 14.;
1532
1533 // --- Parameters to include in GSMIXT, relative to methane (CH4)
1534
1535 amet[0] = 12.01;
1536 amet[1] = 1.;
1537 zmet[0] = 6.;
1538 zmet[1] = 1.;
1539 densmet = 7.17e-4;
1540 nlmatmet = -2;
1541 wmatmet[0] = 1.;
1542 wmatmet[1] = 4.;
1543
1544 // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1545
1546 agri = 63.54;
1547 zgri = 29.;
1548 densgri = 8.96;
1549 radlgri = 1.43;
1550
1551 // --- Parameters to include in GSMATE related to aluminium sheet
1552
1553 aal = 26.98;
1554 zal = 13.;
1555 densal = 2.7;
1556 radlal = 8.9;
1557
1558 // --- Glass parameters
1559
1560 Float_t aglass[5]={12.01, 28.09, 16., 10.8, 23.};
1561 Float_t zglass[5]={ 6., 14., 8., 5., 11.};
1562 Float_t wglass[5]={ 0.5, 0.105, 0.355, 0.03, 0.01};
1563 Float_t dglass=1.74;
1564
1565
1566 AliMaterial(1, "Air $", 14.61, 7.3, .001205, 30420., 67500);
1567 AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1568 AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1569 AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1570 AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1571 AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1572 AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1573 AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1574 AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1575 AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1576 AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1577 AliMaterial(31, "COPPER$", 63.54, 29., 8.96, 1.4, 0.);
1578
1579 tmaxfd = -10.;
1580 stemax = -.1;
1581 deemax = -.2;
1582 epsil = .001;
1583 stmin = -.001;
1584
1585 AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1586 AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1587 AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1588 AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1589 AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1590 AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1591 AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1592 AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1593 AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1594 AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1595 AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1596 AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1597
1598
1599 gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1600 gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1601 gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1602 gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1603 gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1604 gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1605 gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1606 gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1607 gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1608 gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1609 gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1610}
1611
1612//___________________________________________
1613
1614Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1615{
1616
1617 //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1618
1619 Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
1620 6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
1621 7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1622
1623
1624 Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1625 2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1626 2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1627 1.72,1.53};
1628
1629 Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1630 0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1631 0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1632 1.714,1.498};
1633 Float_t xe=ene;
1634 Int_t j=Int_t(xe*10)-49;
1635 Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1636 Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1637
1638 //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1639 //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1640
1641 Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1642 Float_t tanin=sinin/pdoti;
1643
1644 Float_t c1=cn*cn-ck*ck-sinin*sinin;
1645 Float_t c2=4*cn*cn*ck*ck;
1646 Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1647 Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1648
1649 Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1650 Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1651
1652
1653 //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1654 //B.J. STAGG APPLIED OPTICS, 30(1991),4113
1655
1656 Float_t sigraf=18.;
1657 Float_t lamb=1240/ene;
1658 Float_t fresn;
1659
1660 Float_t rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1661
1662 if(pola)
1663 {
1664 Float_t pdotr=0.8; //DEGREE OF POLARIZATION : 1->P , -1->S
1665 fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1666 }
1667 else
1668 fresn=0.5*(rp+rs);
1669
1670 fresn = fresn*rO;
1671 return(fresn);
1672}
1673
1674//__________________________________________
1675Float_t AliRICH::AbsoCH4(Float_t x)
1676{
1677
1678 //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1679 Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.}; //MB X 10^22
1680 //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1681 Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1682 const Float_t kLosch=2.686763E19; // LOSCHMIDT NUMBER IN CM-3
1683 const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1684 Float_t pn=kPressure/760.;
1685 Float_t tn=kTemperature/273.16;
1686
1687
1688// ------- METHANE CROSS SECTION -----------------
1689// ASTROPH. J. 214, L47 (1978)
1690
1691 Float_t sm=0;
1692 if (x<7.75)
1693 sm=.06e-22;
1694
1695 if(x>=7.75 && x<=8.1)
1696 {
1697 Float_t c0=-1.655279e-1;
1698 Float_t c1=6.307392e-2;
1699 Float_t c2=-8.011441e-3;
1700 Float_t c3=3.392126e-4;
1701 sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1702 }
1703
1704 if (x> 8.1)
1705 {
1706 Int_t j=0;
1707 while (x<=em[j] && x>=em[j+1])
1708 {
1709 j++;
1710 Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1711 sm=(sch4[j]+a*(x-em[j]))*1e-22;
1712 }
1713 }
1714
1715 Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1716 Float_t abslm=1./sm/dm;
1717
1718// ------- ISOBUTHANE CROSS SECTION --------------
1719// i-C4H10 (ai) abs. length from curves in
1720// Lu-McDonald paper for BARI RICH workshop .
1721// -----------------------------------------------------------
1722
1723 Float_t ai;
1724 Float_t absli;
1725 if (kIgas2 != 0)
1726 {
1727 if (x<7.25)
1728 ai=100000000.;
1729
1730 if(x>=7.25 && x<7.375)
1731 ai=24.3;
1732
1733 if(x>=7.375)
1734 ai=.0000000001;
1735
1736 Float_t si = 1./(ai*kLosch*273.16/293.); // ISOB. CRO.SEC.IN CM2
1737 Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1738 absli =1./si/di;
1739 }
1740 else
1741 absli=1.e18;
1742// ---------------------------------------------------------
1743//
1744// transmission of O2
1745//
1746// y= path in cm, x=energy in eV
1747// so= cross section for UV absorption in cm2
1748// do= O2 molecular density in cm-3
1749// ---------------------------------------------------------
1750
1751 Float_t abslo;
1752 Float_t so=0;
1753 if(x>=6.0)
1754 {
1755 if(x>=6.0 && x<6.5)
1756 {
1757 so=3.392709e-13 * TMath::Exp(2.864104 *x);
1758 so=so*1e-18;
1759 }
1760
1761 if(x>=6.5 && x<7.0)
1762 {
1763 so=2.910039e-34 * TMath::Exp(10.3337*x);
1764 so=so*1e-18;
1765 }
1766
1767
1768 if (x>=7.0)
1769 {
1770 Float_t a0=-73770.76;
1771 Float_t a1=46190.69;
1772 Float_t a2=-11475.44;
1773 Float_t a3=1412.611;
1774 Float_t a4=-86.07027;
1775 Float_t a5=2.074234;
1776 so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1777 so=so*1e-18;
1778 }
1779
1780 Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1781 abslo=1./so/dox;
1782 }
1783 else
1784 abslo=1.e18;
1785// ---------------------------------------------------------
1786//
1787// transmission of H2O
1788//
1789// y= path in cm, x=energy in eV
1790// sw= cross section for UV absorption in cm2
1791// dw= H2O molecular density in cm-3
1792// ---------------------------------------------------------
1793
1794 Float_t abslw;
1795
1796 Float_t b0=29231.65;
1797 Float_t b1=-15807.74;
1798 Float_t b2=3192.926;
1799 Float_t b3=-285.4809;
1800 Float_t b4=9.533944;
1801
1802 if(x>6.75)
1803 {
1804 Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1805 sw=sw*1e-18;
1806 Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1807 abslw=1./sw/dw;
1808 }
1809 else
1810 abslw=1.e18;
1811
1812// ---------------------------------------------------------
1813
1814 Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1815 return (alength);
1816}
1817
1818
1819
1820//___________________________________________
1821Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1822{
1823
1824// Default value
1825
1826 return 9999;
1827}
1828
1829//___________________________________________
1830void AliRICH::MakeBranch(Option_t* option, const char *file)
1831{
1832 // Create Tree branches for the RICH.
1833
1834 const Int_t kBufferSize = 4000;
1835 char branchname[20];
1836
1837 AliDetector::MakeBranch(option,file);
1838
1839 const char *cH = strstr(option,"H");
1840 const char *cD = strstr(option,"D");
1841 const char *cR = strstr(option,"R");
1842 const char *cS = strstr(option,"S");
1843
1844
1845 if (cH) {
1846 sprintf(branchname,"%sCerenkov",GetName());
1847 if (fCerenkovs && gAlice->TreeH()) {
1848 //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1849 MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1850 //branch->SetAutoDelete(kFALSE);
1851 }
1852 sprintf(branchname,"%sSDigits",GetName());
1853 if (fSDigits && gAlice->TreeH()) {
1854 //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1855 MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1856 //branch->SetAutoDelete(kFALSE);
1857 //printf("Making branch %sSDigits in TreeH\n",GetName());
1858 }
1859 }
1860
1861 if (cS) {
1862 sprintf(branchname,"%sSDigits",GetName());
1863 if (fSDigits && gAlice->TreeS()) {
1864 //TBranch* branch = MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1865 MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1866 //branch->SetAutoDelete(kFALSE);
1867 //printf("Making branch %sSDigits in TreeS\n",GetName());
1868 }
1869 }
1870
1871 if (cD) {
1872 //
1873 // one branch for digits per chamber
1874 //
1875 Int_t i;
1876
1877 for (i=0; i<kNCH ;i++) {
1878 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1879 if (fDchambers && gAlice->TreeD()) {
1880 //TBranch* branch = MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1881 MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1882 //branch->SetAutoDelete(kFALSE);
1883 //printf("Making Branch %sDigits%d\n",GetName(),i+1);
1884 }
1885 }
1886 }
1887
1888 if (cR) {
1889 //
1890 // one branch for raw clusters per chamber
1891 //
1892
1893 //printf("Called MakeBranch for TreeR\n");
1894
1895 Int_t i;
1896
1897 for (i=0; i<kNCH ;i++) {
1898 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1899 if (fRawClusters && gAlice->TreeR()) {
1900 //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1901 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1902 //branch->SetAutoDelete(kFALSE);
1903 }
1904 }
1905 //
1906 // one branch for rec hits per chamber
1907 //
1908 for (i=0; i<kNCH ;i++) {
1909 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1910 if (fRecHits1D && gAlice->TreeR()) {
1911 //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1912 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1913 //branch->SetAutoDelete(kFALSE);
1914 }
1915 }
1916 for (i=0; i<kNCH ;i++) {
1917 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1918 if (fRecHits3D && gAlice->TreeR()) {
1919 MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
1920 //branch->SetAutoDelete(kFALSE);
1921 }
1922 }
1923 }
1924}
1925
1926//___________________________________________
1927void AliRICH::SetTreeAddress()
1928{
1929 // Set branch address for the Hits and Digits Tree.
1930 char branchname[20];
1931 Int_t i;
1932
1933 AliDetector::SetTreeAddress();
1934
1935 TBranch *branch;
1936 TTree *treeH = gAlice->TreeH();
1937 TTree *treeD = gAlice->TreeD();
1938 TTree *treeR = gAlice->TreeR();
1939 TTree *treeS = gAlice->TreeS();
1940
1941 if (treeH) {
1942 if (fCerenkovs) {
1943 branch = treeH->GetBranch("RICHCerenkov");
1944 if (branch) branch->SetAddress(&fCerenkovs);
1945 }
1946 if (fSDigits) {
1947 branch = treeH->GetBranch("RICHSDigits");
1948 if (branch)
1949 {
1950 branch->SetAddress(&fSDigits);
1951 //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits);
1952 }
1953 }
1954 }
1955
1956 if (treeS) {
1957 if (fSDigits) {
1958 branch = treeS->GetBranch("RICHSDigits");
1959 if (branch)
1960 {
1961 branch->SetAddress(&fSDigits);
1962 //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits);
1963 }
1964 }
1965 }
1966
1967
1968 if (treeD) {
1969 for (int i=0; i<kNCH; i++) {
1970 sprintf(branchname,"%sDigits%d",GetName(),i+1);
1971 if (fDchambers) {
1972 branch = treeD->GetBranch(branchname);
1973 if (branch) branch->SetAddress(&((*fDchambers)[i]));
1974 }
1975 }
1976 }
1977 if (treeR) {
1978 for (i=0; i<kNCH; i++) {
1979 sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1980 if (fRawClusters) {
1981 branch = treeR->GetBranch(branchname);
1982 if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1983 }
1984 }
1985
1986 for (i=0; i<kNCH; i++) {
1987 sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1988 if (fRecHits1D) {
1989 branch = treeR->GetBranch(branchname);
1990 if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
1991 }
1992 }
1993
1994 for (i=0; i<kNCH; i++) {
1995 sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1996 if (fRecHits3D) {
1997 branch = treeR->GetBranch(branchname);
1998 if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
1999 }
2000 }
2001
2002 }
2003}
2004//___________________________________________
2005void AliRICH::ResetHits()
2006{
2007 // Reset number of clusters and the cluster array for this detector
2008 AliDetector::ResetHits();
2009 fNSDigits = 0;
2010 fNcerenkovs = 0;
2011 if (fSDigits) fSDigits->Clear();
2012 if (fCerenkovs) fCerenkovs->Clear();
2013}
2014
2015
2016//____________________________________________
2017void AliRICH::ResetDigits()
2018{
2019 //
2020 // Reset number of digits and the digits array for this detector
2021 //
2022 for ( int i=0;i<kNCH;i++ ) {
2023 //PH if (fDchambers && (*fDchambers)[i]) (*fDchambers)[i]->Clear();
2024 if (fDchambers && fDchambers->At(i)) fDchambers->At(i)->Clear();
2025 if (fNdch) fNdch[i]=0;
2026 }
2027}
2028
2029//____________________________________________
2030void AliRICH::ResetRawClusters()
2031{
2032 //
2033 // Reset number of raw clusters and the raw clust array for this detector
2034 //
2035 for ( int i=0;i<kNCH;i++ ) {
2036 //PH if ((*fRawClusters)[i]) ((TClonesArray*)(*fRawClusters)[i])->Clear();
2037 if (fRawClusters->At(i)) ((TClonesArray*)fRawClusters->At(i))->Clear();
2038 if (fNrawch) fNrawch[i]=0;
2039 }
2040}
2041
2042//____________________________________________
2043void AliRICH::ResetRecHits1D()
2044{
2045 //
2046 // Reset number of raw clusters and the raw clust array for this detector
2047 //
2048
2049 for ( int i=0;i<kNCH;i++ ) {
2050 //PH if ((*fRecHits1D)[i]) ((TClonesArray*)(*fRecHits1D)[i])->Clear();
2051 if (fRecHits1D->At(i)) ((TClonesArray*)fRecHits1D->At(i))->Clear();
2052 if (fNrechits1D) fNrechits1D[i]=0;
2053 }
2054}
2055
2056//____________________________________________
2057void AliRICH::ResetRecHits3D()
2058{
2059 //
2060 // Reset number of raw clusters and the raw clust array for this detector
2061 //
2062
2063 for ( int i=0;i<kNCH;i++ ) {
2064 //PH if ((*fRecHits3D)[i]) ((TClonesArray*)(*fRecHits3D)[i])->Clear();
2065 if (fRecHits3D->At(i)) ((TClonesArray*)fRecHits3D->At(i))->Clear();
2066 if (fNrechits3D) fNrechits3D[i]=0;
2067 }
2068}
2069
2070
2071//___________________________________________
2072void AliRICH::StepManager()
2073{
2074
2075// Full Step Manager
2076
2077 Int_t copy, id;
2078 static Int_t idvol;
2079 static Int_t vol[2];
2080 Int_t ipart;
2081 static Float_t hits[22];
2082 static Float_t ckovData[19];
2083 TLorentzVector position;
2084 TLorentzVector momentum;
2085 Float_t pos[3];
2086 Float_t mom[4];
2087 Float_t localPos[3];
2088 Float_t localMom[4];
2089 Float_t localTheta,localPhi;
2090 Float_t theta,phi;
2091 Float_t destep, step;
2092 Float_t ranf[2];
2093 Int_t nPads;
2094 Float_t coscerenkov;
2095 static Float_t eloss, xhit, yhit, tlength;
2096 const Float_t kBig=1.e10;
2097
2098 TClonesArray &lhits = *fHits;
2099 TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
2100
2101 //if (current->Energy()>1)
2102 //{
2103
2104 // Only gas gap inside chamber
2105 // Tag chambers and record hits when track enters
2106
2107 idvol=-1;
2108 id=gMC->CurrentVolID(copy);
2109 Float_t cherenkovLoss=0;
2110 //gAlice->KeepTrack(gAlice->CurrentTrack());
2111
2112 gMC->TrackPosition(position);
2113 pos[0]=position(0);
2114 pos[1]=position(1);
2115 pos[2]=position(2);
2116 //bzero((char *)ckovData,sizeof(ckovData)*19);
2117 ckovData[1] = pos[0]; // X-position for hit
2118 ckovData[2] = pos[1]; // Y-position for hit
2119 ckovData[3] = pos[2]; // Z-position for hit
2120 ckovData[6] = 0; // dummy track length
2121 //ckovData[11] = gAlice->CurrentTrack();
2122
2123 //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
2124
2125 //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
2126
2127 /********************Store production parameters for Cerenkov photons************************/
2128//is it a Cerenkov photon?
2129 if (gMC->TrackPid() == 50000050) {
2130
2131 //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
2132 //{
2133 Float_t ckovEnergy = current->Energy();
2134 //energy interval for tracking
2135 if (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
2136 //if (ckovEnergy > 0)
2137 {
2138 if (gMC->IsTrackEntering()){ //is track entering?
2139 //printf("Track entered (1)\n");
2140 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2141 { //is it in freo?
2142 if (gMC->IsNewTrack()){ //is it the first step?
2143 //printf("I'm in!\n");
2144 Int_t mother = current->GetFirstMother();
2145
2146 //printf("Second Mother:%d\n",current->GetSecondMother());
2147
2148 ckovData[10] = mother;
2149 ckovData[11] = gAlice->CurrentTrack();
2150 ckovData[12] = 1; //Media where photon was produced 1->Freon, 2->Quarz
2151 //printf("Produced in FREO\n");
2152 fCkovNumber++;
2153 fFreonProd=1;
2154 //printf("Index: %d\n",fCkovNumber);
2155 } //first step question
2156 } //freo question
2157
2158 if (gMC->IsNewTrack()){ //is it first step?
2159 if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy)) //is it in quarz?
2160 {
2161 ckovData[12] = 2;
2162 //printf("Produced in QUAR\n");
2163 } //quarz question
2164 } //first step question
2165
2166 //printf("Before %d\n",fFreonProd);
2167 } //track entering question
2168
2169 if (ckovData[12] == 1) //was it produced in Freon?
2170 //if (fFreonProd == 1)
2171 {
2172 if (gMC->IsTrackEntering()){ //is track entering?
2173 //printf("Track entered (2)\n");
2174 //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
2175 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
2176 if (gMC->VolId("META")==gMC->CurrentVolID(copy)) //is it in gap?
2177 {
2178 //printf("Got in META\n");
2179 gMC->TrackMomentum(momentum);
2180 mom[0]=momentum(0);
2181 mom[1]=momentum(1);
2182 mom[2]=momentum(2);
2183 mom[3]=momentum(3);
2184 // Z-position for hit
2185
2186
2187 /**************** Photons lost in second grid have to be calculated by hand************/
2188
2189 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2190 Float_t t = (1. - .025 / cophi) * (1. - .05 / cophi);
2191 gMC->Rndm(ranf, 1);
2192 //printf("grid calculation:%f\n",t);
2193 if (ranf[0] > t) {
2194 gMC->StopTrack();
2195 ckovData[13] = 5;
2196 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2197 //printf("Added One (1)!\n");
2198 //printf("Lost one in grid\n");
2199 }
2200 /**********************************************************************************/
2201 } //gap
2202
2203 //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
2204 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2205 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy)) //is it in csi?
2206 {
2207 //printf("Got in CSI\n");
2208 gMC->TrackMomentum(momentum);
2209 mom[0]=momentum(0);
2210 mom[1]=momentum(1);
2211 mom[2]=momentum(2);
2212 mom[3]=momentum(3);
2213
2214 /********* Photons lost by Fresnel reflection have to be calculated by hand********/
2215 /***********************Cerenkov phtons (always polarised)*************************/
2216
2217 Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2218 Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
2219 gMC->Rndm(ranf, 1);
2220 if (ranf[0] < t) {
2221 gMC->StopTrack();
2222 ckovData[13] = 6;
2223 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2224 //printf("Added One (2)!\n");
2225 //printf("Lost by Fresnel\n");
2226 }
2227 /**********************************************************************************/
2228 }
2229 } //track entering?
2230
2231
2232 /********************Evaluation of losses************************/
2233 /******************still in the old fashion**********************/
2234
2235 TArrayI procs;
2236 Int_t i1 = gMC->StepProcesses(procs); //number of physics mechanisms acting on the particle
2237 for (Int_t i = 0; i < i1; ++i) {
2238 // Reflection loss
2239 if (procs[i] == kPLightReflection) { //was it reflected
2240 ckovData[13]=10;
2241 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2242 ckovData[13]=1;
2243 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2244 ckovData[13]=2;
2245 //gMC->StopTrack();
2246 //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2247 } //reflection question
2248
2249 // Absorption loss
2250 else if (procs[i] == kPLightAbsorption) { //was it absorbed?
2251 //printf("Got in absorption\n");
2252 ckovData[13]=20;
2253 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2254 ckovData[13]=11;
2255 if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2256 ckovData[13]=12;
2257 if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
2258 ckovData[13]=13;
2259 if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
2260 ckovData[13]=13;
2261
2262 if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
2263 ckovData[13]=15;
2264
2265 // CsI inefficiency
2266 if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
2267 ckovData[13]=16;
2268 }
2269 gMC->StopTrack();
2270 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2271 //printf("Added One (3)!\n");
2272 //printf("Added cerenkov %d\n",fCkovNumber);
2273 } //absorption question
2274
2275
2276 // Photon goes out of tracking scope
2277 else if (procs[i] == kPStop) { //is it below energy treshold?
2278 ckovData[13]=21;
2279 gMC->StopTrack();
2280 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2281 //printf("Added One (4)!\n");
2282 } // energy treshold question
2283 } //number of mechanisms cycle
2284 /**********************End of evaluation************************/
2285 } //freon production question
2286 } //energy interval question
2287 //}//inside the proximity gap question
2288 } //cerenkov photon question
2289
2290 /**************************************End of Production Parameters Storing*********************/
2291
2292
2293 /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
2294
2295 if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
2296 //printf("Cerenkov\n");
2297
2298 //if (gMC->TrackPid() == 50000051)
2299 //printf("Tracking a feedback\n");
2300
2301 if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
2302 {
2303 //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
2304 //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2305 //printf("Got in CSI\n");
2306 //printf("Tracking a %d\n",gMC->TrackPid());
2307 if (gMC->Edep() > 0.){
2308 gMC->TrackPosition(position);
2309 gMC->TrackMomentum(momentum);
2310 pos[0]=position(0);
2311 pos[1]=position(1);
2312 pos[2]=position(2);
2313 mom[0]=momentum(0);
2314 mom[1]=momentum(1);
2315 mom[2]=momentum(2);
2316 mom[3]=momentum(3);
2317 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2318 Double_t rt = TMath::Sqrt(tc);
2319 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2320 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2321 gMC->Gmtod(pos,localPos,1);
2322 gMC->Gmtod(mom,localMom,2);
2323
2324 gMC->CurrentVolOffID(2,copy);
2325 vol[0]=copy;
2326 idvol=vol[0]-1;
2327
2328 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2329 //->Sector(localPos[0], localPos[2]);
2330 //printf("Sector:%d\n",sector);
2331
2332 /*if (gMC->TrackPid() == 50000051){
2333 fFeedbacks++;
2334 printf("Feedbacks:%d\n",fFeedbacks);
2335 }*/
2336
2337 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2338 ((AliRICHChamber*)fChambers->At(idvol))
2339 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2340 if(idvol<kNCH) {
2341 ckovData[0] = gMC->TrackPid(); // particle type
2342 ckovData[1] = pos[0]; // X-position for hit
2343 ckovData[2] = pos[1]; // Y-position for hit
2344 ckovData[3] = pos[2]; // Z-position for hit
2345 ckovData[4] = theta; // theta angle of incidence
2346 ckovData[5] = phi; // phi angle of incidence
2347 ckovData[8] = (Float_t) fNSDigits; // first sdigit
2348 ckovData[9] = -1; // last pad hit
2349 ckovData[13] = 4; // photon was detected
2350 ckovData[14] = mom[0];
2351 ckovData[15] = mom[1];
2352 ckovData[16] = mom[2];
2353
2354 destep = gMC->Edep();
2355 gMC->SetMaxStep(kBig);
2356 cherenkovLoss += destep;
2357 ckovData[7]=cherenkovLoss;
2358
2359 nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
2360
2361 if (fNSDigits > (Int_t)ckovData[8]) {
2362 ckovData[8]= ckovData[8]+1;
2363 ckovData[9]= (Float_t) fNSDigits;
2364 }
2365
2366 //printf("Cerenkov loss: %f\n", cherenkovLoss);
2367
2368 ckovData[17] = nPads;
2369 //printf("nPads:%d",nPads);
2370
2371 //TClonesArray *Hits = RICH->Hits();
2372 AliRICHHit *mipHit = (AliRICHHit*) (fHits->UncheckedAt(0));
2373 if (mipHit)
2374 {
2375 mom[0] = current->Px();
2376 mom[1] = current->Py();
2377 mom[2] = current->Pz();
2378 Float_t mipPx = mipHit->MomX();
2379 Float_t mipPy = mipHit->MomY();
2380 Float_t mipPz = mipHit->MomZ();
2381
2382 Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
2383 Float_t rt = TMath::Sqrt(r);
2384 Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
2385 Float_t mipRt = TMath::Sqrt(mipR);
2386 if ((rt*mipRt) > 0)
2387 {
2388 coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
2389 }
2390 else
2391 {
2392 coscerenkov = 0;
2393 }
2394 Float_t cherenkov = TMath::ACos(coscerenkov);
2395 ckovData[18]=cherenkov;
2396 }
2397 //if (sector != -1)
2398 //{
2399 AddHit(gAlice->CurrentTrack(),vol,ckovData);
2400 AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2401 //printf("Added One (5)!\n");
2402 //}
2403 }
2404 }
2405 }
2406 }
2407
2408 /***********************************************End of photon hits*********************************************/
2409
2410
2411 /**********************************************Charged particles treatment*************************************/
2412
2413 else if (gMC->TrackCharge())
2414 //else if (1 == 1)
2415 {
2416//If MIP
2417 /*if (gMC->IsTrackEntering())
2418 {
2419 hits[13]=20;//is track entering?
2420 }*/
2421 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2422 {
2423 gMC->TrackMomentum(momentum);
2424 mom[0]=momentum(0);
2425 mom[1]=momentum(1);
2426 mom[2]=momentum(2);
2427 mom[3]=momentum(3);
2428 hits [19] = mom[0];
2429 hits [20] = mom[1];
2430 hits [21] = mom[2];
2431 fFreonProd=1;
2432 }
2433
2434 if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2435// Get current particle id (ipart), track position (pos) and momentum (mom)
2436
2437 gMC->CurrentVolOffID(3,copy);
2438 vol[0]=copy;
2439 idvol=vol[0]-1;
2440
2441 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2442 //->Sector(localPos[0], localPos[2]);
2443 //printf("Sector:%d\n",sector);
2444
2445 gMC->TrackPosition(position);
2446 gMC->TrackMomentum(momentum);
2447 pos[0]=position(0);
2448 pos[1]=position(1);
2449 pos[2]=position(2);
2450 mom[0]=momentum(0);
2451 mom[1]=momentum(1);
2452 mom[2]=momentum(2);
2453 mom[3]=momentum(3);
2454 gMC->Gmtod(pos,localPos,1);
2455 gMC->Gmtod(mom,localMom,2);
2456
2457 ipart = gMC->TrackPid();
2458 //
2459 // momentum loss and steplength in last step
2460 destep = gMC->Edep();
2461 step = gMC->TrackStep();
2462
2463 //
2464 // record hits when track enters ...
2465 if( gMC->IsTrackEntering()) {
2466// gMC->SetMaxStep(fMaxStepGas);
2467 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2468 Double_t rt = TMath::Sqrt(tc);
2469 theta = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2470 phi = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2471
2472
2473 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2474 Double_t localRt = TMath::Sqrt(localTc);
2475 localTheta = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2476 localPhi = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2477
2478 hits[0] = Float_t(ipart); // particle type
2479 hits[1] = localPos[0]; // X-position for hit
2480 hits[2] = localPos[1]; // Y-position for hit
2481 hits[3] = localPos[2]; // Z-position for hit
2482 hits[4] = localTheta; // theta angle of incidence
2483 hits[5] = localPhi; // phi angle of incidence
2484 hits[8] = (Float_t) fNSDigits; // first sdigit
2485 hits[9] = -1; // last pad hit
2486 hits[13] = fFreonProd; // did id hit the freon?
2487 hits[14] = mom[0];
2488 hits[15] = mom[1];
2489 hits[16] = mom[2];
2490 hits[18] = 0; // dummy cerenkov angle
2491
2492 tlength = 0;
2493 eloss = 0;
2494 fFreonProd = 0;
2495
2496 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2497
2498
2499 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2500 xhit = localPos[0];
2501 yhit = localPos[2];
2502 // Only if not trigger chamber
2503 if(idvol<kNCH) {
2504 //
2505 // Initialize hit position (cursor) in the segmentation model
2506 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2507 ((AliRICHChamber*)fChambers->At(idvol))
2508 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2509 }
2510 }
2511
2512 //
2513 // Calculate the charge induced on a pad (disintegration) in case
2514 //
2515 // Mip left chamber ...
2516 if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2517 gMC->SetMaxStep(kBig);
2518 eloss += destep;
2519 tlength += step;
2520
2521
2522 // Only if not trigger chamber
2523 if(idvol<kNCH) {
2524 if (eloss > 0)
2525 {
2526 if(gMC->TrackPid() == kNeutron)
2527 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2528 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2529 hits[17] = nPads;
2530 //printf("nPads:%d",nPads);
2531 }
2532 }
2533
2534 hits[6]=tlength;
2535 hits[7]=eloss;
2536 if (fNSDigits > (Int_t)hits[8]) {
2537 hits[8]= hits[8]+1;
2538 hits[9]= (Float_t) fNSDigits;
2539 }
2540
2541 //if(sector !=-1)
2542 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2543 eloss = 0;
2544 //
2545 // Check additional signal generation conditions
2546 // defined by the segmentation
2547 // model (boundary crossing conditions)
2548 } else if
2549 //PH (((AliRICHChamber*) (*fChambers)[idvol])
2550 (((AliRICHChamber*)fChambers->At(idvol))
2551 ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2552 {
2553 //PH ((AliRICHChamber*) (*fChambers)[idvol])
2554 ((AliRICHChamber*)fChambers->At(idvol))
2555 ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2556 if (eloss > 0)
2557 {
2558 if(gMC->TrackPid() == kNeutron)
2559 printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2560 nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2561 hits[17] = nPads;
2562 //printf("Npads:%d",NPads);
2563 }
2564 xhit = localPos[0];
2565 yhit = localPos[2];
2566 eloss = destep;
2567 tlength += step ;
2568 //
2569 // nothing special happened, add up energy loss
2570 } else {
2571 eloss += destep;
2572 tlength += step ;
2573 }
2574 }
2575 }
2576 /*************************************************End of MIP treatment**************************************/
2577 //}
2578}
2579
2580void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
2581{
2582
2583//
2584// Loop on chambers and on cathode planes
2585//
2586 for (Int_t icat=1;icat<2;icat++) {
2587 gAlice->ResetDigits();
2588 gAlice->TreeD()->GetEvent(0);
2589 for (Int_t ich=0;ich<kNCH;ich++) {
2590 //PH AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
2591 AliRICHChamber* iChamber=(AliRICHChamber*)fChambers->At(ich);
2592 TClonesArray *pRICHdigits = this->DigitsAddress(ich);
2593 if (pRICHdigits == 0)
2594 continue;
2595 //
2596 // Get ready the current chamber stuff
2597 //
2598 AliRICHResponse* response = iChamber->GetResponseModel();
2599 AliSegmentation* seg = iChamber->GetSegmentationModel();
2600 AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2601 if (seg) {
2602 rec->SetSegmentation(seg);
2603 rec->SetResponse(response);
2604 rec->SetDigits(pRICHdigits);
2605 rec->SetChamber(ich);
2606 if (nev==0) rec->CalibrateCOG();
2607 rec->FindRawClusters();
2608 }
2609 TClonesArray *fRch;
2610 fRch=RawClustAddress(ich);
2611 fRch->Sort();
2612 } // for ich
2613
2614 gAlice->TreeR()->Fill();
2615 TClonesArray *fRch;
2616 for (int i=0;i<kNCH;i++) {
2617 fRch=RawClustAddress(i);
2618 int nraw=fRch->GetEntriesFast();
2619 printf ("Chamber %d, raw clusters %d\n",i,nraw);
2620 }
2621
2622 ResetRawClusters();
2623
2624 } // for icat
2625
2626 char hname[30];
2627 sprintf(hname,"TreeR%d",nev);
2628 gAlice->TreeR()->Write(hname,kOverwrite,0);
2629 gAlice->TreeR()->Reset();
2630
2631 //gObjectTable->Print();
2632}
2633
2634AliRICHSDigit* AliRICH::FirstPad(AliRICHHit* hit,TClonesArray *clusters )
2635{
2636//
2637 // Initialise the pad iterator
2638 // Return the address of the first sdigit for hit
2639 TClonesArray *theClusters = clusters;
2640 Int_t nclust = theClusters->GetEntriesFast();
2641 if (nclust && hit->PHlast() > 0) {
2642 sMaxIterPad=Int_t(hit->PHlast());
2643 sCurIterPad=Int_t(hit->PHfirst());
2644 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2645 } else {
2646 return 0;
2647 }
2648
2649}
2650
2651AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters)
2652{
2653
2654 // Iterates over pads
2655
2656 sCurIterPad++;
2657 if (sCurIterPad <= sMaxIterPad) {
2658 return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2659 } else {
2660 return 0;
2661 }
2662}
2663
2664AliRICH& AliRICH::operator=(const AliRICH& rhs)
2665{
2666// Assignment operator
2667 return *this;
2668
2669}
2670
2671void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
2672{
2673
2674 Int_t NpadX = 162; // number of pads on X
2675 Int_t NpadY = 162; // number of pads on Y
2676
2677 Int_t Pad[162][162];
2678 for (Int_t i=0;i<NpadX;i++) {
2679 for (Int_t j=0;j<NpadY;j++) {
2680 Pad[i][j]=0;
2681 }
2682 }
2683
2684 // Create some histograms
2685
2686 TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
2687 TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
2688 TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
2689 TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
2690 TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
2691 TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
2692 TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
2693 TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
2694 TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
2695 TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
2696 TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
2697 TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
2698 TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
2699 TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
2700 TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
2701 TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
2702 TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2703 TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2704 TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
2705 TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2706 TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2707 TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
2708 TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
2709 TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
2710 TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
2711 //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
2712 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
2713 TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
2714 TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
2715 TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
2716 TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
2717 TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
2718
2719
2720
2721
2722// Start loop over events
2723
2724 Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
2725 Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
2726 Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
2727 TRandom* random=0;
2728
2729 for (int nev=0; nev<= evNumber2; nev++) {
2730 Int_t nparticles = gAlice->GetEvent(nev);
2731
2732
2733 printf ("Event number : %d\n",nev);
2734 printf ("Number of particles: %d\n",nparticles);
2735 if (nev < evNumber1) continue;
2736 if (nparticles <= 0) return;
2737
2738// Get pointers to RICH detector and Hits containers
2739
2740 AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
2741
2742 TTree *treeH = gAlice->TreeH();
2743 Int_t ntracks =(Int_t) treeH->GetEntries();
2744
2745// Start loop on tracks in the hits containers
2746
2747 for (Int_t track=0; track<ntracks;track++) {
2748 printf ("Processing Track: %d\n",track);
2749 gAlice->ResetHits();
2750 treeH->GetEvent(track);
2751
2752 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2753 mHit;
2754 mHit=(AliRICHHit*)pRICH->NextHit())
2755 {
2756 //Int_t nch = mHit->fChamber; // chamber number
2757 //Float_t x = mHit->X(); // x-pos of hit
2758 //Float_t y = mHit->Z(); // y-pos
2759 //Float_t z = mHit->Y();
2760 //Float_t phi = mHit->Phi(); //Phi angle of incidence
2761 Float_t theta = mHit->Theta(); //Theta angle of incidence
2762 Float_t px = mHit->MomX();
2763 Float_t py = mHit->MomY();
2764 Int_t index = mHit->Track();
2765 Int_t particle = (Int_t)(mHit->Particle());
2766 Float_t R;
2767 Float_t PTfinal;
2768 Float_t PTvertex;
2769
2770 TParticle *current = gAlice->Particle(index);
2771
2772 //Float_t energy=current->Energy();
2773
2774 R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
2775 PTfinal=TMath::Sqrt(px*px + py*py);
2776 PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
2777
2778
2779
2780 if (TMath::Abs(particle) < 10000000)
2781 {
2782 hitsTheta->Fill(theta,(float) 1);
2783 if (R<5)
2784 {
2785 if (PTvertex>.5 && PTvertex<=1)
2786 {
2787 hitsTheta500MeV->Fill(theta,(float) 1);
2788 }
2789 if (PTvertex>1 && PTvertex<=2)
2790 {
2791 hitsTheta1GeV->Fill(theta,(float) 1);
2792 }
2793 if (PTvertex>2 && PTvertex<=3)
2794 {
2795 hitsTheta2GeV->Fill(theta,(float) 1);
2796 }
2797 if (PTvertex>3)
2798 {
2799 hitsTheta3GeV->Fill(theta,(float) 1);
2800 }
2801 }
2802
2803 }
2804
2805 //if (nch == 3)
2806 //{
2807
2808 //printf("Particle type: %d\n",current->GetPdgCode());
2809 if (TMath::Abs(particle) < 50000051)
2810 {
2811 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
2812 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
2813 {
2814 //gMC->Rndm(&random, 1);
2815 if (random->Rndm() < .1)
2816 production->Fill(current->Vz(),R,(float) 1);
2817 if (TMath::Abs(particle) == 50000050)
2818 //if (TMath::Abs(particle) > 50000000)
2819 {
2820 photons +=1;
2821 if (R<5)
2822 {
2823 primaryphotons +=1;
2824 if (current->Energy()>0.001)
2825 highprimaryphotons +=1;
2826 }
2827 }
2828 if (TMath::Abs(particle) == 2112)
2829 {
2830 neutron +=1;
2831 if (current->Energy()>0.0001)
2832 highneutrons +=1;
2833 }
2834 }
2835 if (TMath::Abs(particle) < 50000000)
2836 {
2837 production->Fill(current->Vz(),R,(float) 1);
2838 //printf("Adding %d at %f\n",particle,R);
2839 }
2840 //mip->Fill(x,y,(float) 1);
2841 }
2842
2843 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2844 {
2845 if (R<5)
2846 {
2847 pionptspectravertex->Fill(PTvertex,(float) 1);
2848 pionptspectrafinal->Fill(PTfinal,(float) 1);
2849 }
2850 }
2851
2852 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2853 || TMath::Abs(particle)==311)
2854 {
2855 if (R<5)
2856 {
2857 kaonptspectravertex->Fill(PTvertex,(float) 1);
2858 kaonptspectrafinal->Fill(PTfinal,(float) 1);
2859 }
2860 }
2861
2862
2863 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2864 {
2865 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2866 //printf ("fParticle: %d, PDG code:%d\n",particle,current->GetPdgCode());
2867 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2868 pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2869 if (R>250 && R<450)
2870 {
2871 pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2872 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2873 }
2874 //printf("Pion mass: %e\n",current->GetCalcMass());
2875 pion +=1;
2876 if (TMath::Abs(particle)==211)
2877 {
2878 chargedpions +=1;
2879 if (R<5)
2880 {
2881 primarypions +=1;
2882 if (current->Energy()>1)
2883 highprimarypions +=1;
2884 }
2885 }
2886 }
2887 if (TMath::Abs(particle)==2212)
2888 {
2889 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2890 //ptspectra->Fill(Pt,(float) 1);
2891 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2892 protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2893 if (R>250 && R<450)
2894 protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2895 //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
2896 proton +=1;
2897 }
2898 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2899 || TMath::Abs(particle)==311)
2900 {
2901 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2902 //ptspectra->Fill(Pt,(float) 1);
2903 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2904 kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2905 if (R>250 && R<450)
2906 kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2907 //printf("Kaon mass: %e\n",current->GetCalcMass());
2908 kaon +=1;
2909 if (TMath::Abs(particle)==321)
2910 {
2911 chargedkaons +=1;
2912 if (R<5)
2913 {
2914 primarykaons +=1;
2915 if (current->Energy()>1)
2916 highprimarykaons +=1;
2917 }
2918 }
2919 }
2920 if (TMath::Abs(particle)==11)
2921 {
2922 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2923 //ptspectra->Fill(Pt,(float) 1);
2924 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2925 electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2926 if (R>250 && R<450)
2927 electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2928 //printf("Electron mass: %e\n",current->GetCalcMass());
2929 if (particle == 11)
2930 electron +=1;
2931 if (particle == -11)
2932 positron +=1;
2933 }
2934 if (TMath::Abs(particle)==13)
2935 {
2936 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2937 //ptspectra->Fill(Pt,(float) 1);
2938 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2939 muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2940 if (R>250 && R<450)
2941 muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2942 //printf("Muon mass: %e\n",current->GetCalcMass());
2943 muon +=1;
2944 }
2945 if (TMath::Abs(particle)==2112)
2946 {
2947 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2948 //ptspectra->Fill(Pt,(float) 1);
2949 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2950 neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2951 if (R>250 && R<450)
2952 {
2953 neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2954 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2955 }
2956 //printf("Neutron mass: %e\n",current->GetCalcMass());
2957 neutron +=1;
2958 }
2959 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
2960 {
2961 if (current->Energy()-current->GetCalcMass()>1)
2962 {
2963 chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2964 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2965 chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2966 if (R>250 && R<450)
2967 chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2968 }
2969 }
2970 //printf("Hits:%d\n",hit);
2971 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
2972 // Fill the histograms
2973 //Nh1+=nhits;
2974 //h->Fill(x,y,(float) 1);
2975 //}
2976 //}
2977 }
2978
2979 }
2980
2981 }
2982 // }
2983
2984 //Create canvases, set the view range, show histograms
2985
2986 TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
2987 c2->Divide(2,2);
2988 //c2->SetFillColor(42);
2989
2990 c2->cd(1);
2991 hitsTheta500MeV->SetFillColor(5);
2992 hitsTheta500MeV->Draw();
2993 c2->cd(2);
2994 hitsTheta1GeV->SetFillColor(5);
2995 hitsTheta1GeV->Draw();
2996 c2->cd(3);
2997 hitsTheta2GeV->SetFillColor(5);
2998 hitsTheta2GeV->Draw();
2999 c2->cd(4);
3000 hitsTheta3GeV->SetFillColor(5);
3001 hitsTheta3GeV->Draw();
3002
3003
3004
3005 TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
3006 c15->cd();
3007 production->SetFillColor(42);
3008 production->SetXTitle("z (m)");
3009 production->SetYTitle("R (m)");
3010 production->Draw();
3011
3012 TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
3013 c10->Divide(2,2);
3014 c10->cd(1);
3015 pionptspectravertex->SetFillColor(5);
3016 pionptspectravertex->SetXTitle("Pt (GeV)");
3017 pionptspectravertex->Draw();
3018 c10->cd(2);
3019 pionptspectrafinal->SetFillColor(5);
3020 pionptspectrafinal->SetXTitle("Pt (GeV)");
3021 pionptspectrafinal->Draw();
3022 c10->cd(3);
3023 kaonptspectravertex->SetFillColor(5);
3024 kaonptspectravertex->SetXTitle("Pt (GeV)");
3025 kaonptspectravertex->Draw();
3026 c10->cd(4);
3027 kaonptspectrafinal->SetFillColor(5);
3028 kaonptspectrafinal->SetXTitle("Pt (GeV)");
3029 kaonptspectrafinal->Draw();
3030
3031
3032 TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
3033 c16->Divide(2,1);
3034
3035 c16->cd(1);
3036 //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
3037 electronspectra1->SetFillColor(5);
3038 electronspectra1->SetXTitle("log(GeV)");
3039 electronspectra2->SetFillColor(46);
3040 electronspectra2->SetXTitle("log(GeV)");
3041 electronspectra3->SetFillColor(10);
3042 electronspectra3->SetXTitle("log(GeV)");
3043 //c13->SetLogx();
3044 electronspectra1->Draw();
3045 electronspectra2->Draw("same");
3046 electronspectra3->Draw("same");
3047
3048 c16->cd(2);
3049 //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
3050 muonspectra1->SetFillColor(5);
3051 muonspectra1->SetXTitle("log(GeV)");
3052 muonspectra2->SetFillColor(46);
3053 muonspectra2->SetXTitle("log(GeV)");
3054 muonspectra3->SetFillColor(10);
3055 muonspectra3->SetXTitle("log(GeV)");
3056 //c14->SetLogx();
3057 muonspectra1->Draw();
3058 muonspectra2->Draw("same");
3059 muonspectra3->Draw("same");
3060
3061 //c16->cd(3);
3062 //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
3063 //neutronspectra1->SetFillColor(42);
3064 //neutronspectra1->SetXTitle("log(GeV)");
3065 //neutronspectra2->SetFillColor(46);
3066 //neutronspectra2->SetXTitle("log(GeV)");
3067 //neutronspectra3->SetFillColor(10);
3068 //neutronspectra3->SetXTitle("log(GeV)");
3069 //c16->SetLogx();
3070 //neutronspectra1->Draw();
3071 //neutronspectra2->Draw("same");
3072 //neutronspectra3->Draw("same");
3073
3074 TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
3075 //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
3076 c9->Divide(2,2);
3077
3078 c9->cd(1);
3079 pionspectra1->SetFillColor(5);
3080 pionspectra1->SetXTitle("log(GeV)");
3081 pionspectra2->SetFillColor(46);
3082 pionspectra2->SetXTitle("log(GeV)");
3083 pionspectra3->SetFillColor(10);
3084 pionspectra3->SetXTitle("log(GeV)");
3085 //c9->SetLogx();
3086 pionspectra1->Draw();
3087 pionspectra2->Draw("same");
3088 pionspectra3->Draw("same");
3089
3090 c9->cd(2);
3091 //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
3092 protonspectra1->SetFillColor(5);
3093 protonspectra1->SetXTitle("log(GeV)");
3094 protonspectra2->SetFillColor(46);
3095 protonspectra2->SetXTitle("log(GeV)");
3096 protonspectra3->SetFillColor(10);
3097 protonspectra3->SetXTitle("log(GeV)");
3098 //c10->SetLogx();
3099 protonspectra1->Draw();
3100 protonspectra2->Draw("same");
3101 protonspectra3->Draw("same");
3102
3103 c9->cd(3);
3104 //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
3105 kaonspectra1->SetFillColor(5);
3106 kaonspectra1->SetXTitle("log(GeV)");
3107 kaonspectra2->SetFillColor(46);
3108 kaonspectra2->SetXTitle("log(GeV)");
3109 kaonspectra3->SetFillColor(10);
3110 kaonspectra3->SetXTitle("log(GeV)");
3111 //c11->SetLogx();
3112 kaonspectra1->Draw();
3113 kaonspectra2->Draw("same");
3114 kaonspectra3->Draw("same");
3115
3116 c9->cd(4);
3117 //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
3118 chargedspectra1->SetFillColor(5);
3119 chargedspectra1->SetXTitle("log(GeV)");
3120 chargedspectra2->SetFillColor(46);
3121 chargedspectra2->SetXTitle("log(GeV)");
3122 chargedspectra3->SetFillColor(10);
3123 chargedspectra3->SetXTitle("log(GeV)");
3124 //c12->SetLogx();
3125 chargedspectra1->Draw();
3126 chargedspectra2->Draw("same");
3127 chargedspectra3->Draw("same");
3128
3129
3130
3131 printf("*****************************************\n");
3132 printf("* Particle * Counts *\n");
3133 printf("*****************************************\n");
3134
3135 printf("* Pions: * %4d *\n",pion);
3136 printf("* Charged Pions: * %4d *\n",chargedpions);
3137 printf("* Primary Pions: * %4d *\n",primarypions);
3138 printf("* Primary Pions (p>1GeV/c): * %4d *\n",highprimarypions);
3139 printf("* Kaons: * %4d *\n",kaon);
3140 printf("* Charged Kaons: * %4d *\n",chargedkaons);
3141 printf("* Primary Kaons: * %4d *\n",primarykaons);
3142 printf("* Primary Kaons (p>1GeV/c): * %4d *\n",highprimarykaons);
3143 printf("* Muons: * %4d *\n",muon);
3144 printf("* Electrons: * %4d *\n",electron);
3145 printf("* Positrons: * %4d *\n",positron);
3146 printf("* Protons: * %4d *\n",proton);
3147 printf("* All Charged: * %4d *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
3148 printf("*****************************************\n");
3149 //printf("* Photons: * %3.1f *\n",photons);
3150 //printf("* Primary Photons: * %3.1f *\n",primaryphotons);
3151 //printf("* Primary Photons (p>1MeV/c):* %3.1f *\n",highprimaryphotons);
3152 //printf("*****************************************\n");
3153 //printf("* Neutrons: * %3.1f *\n",neutron);
3154 //printf("* Neutrons (p>100keV/c): * %3.1f *\n",highneutrons);
3155 //printf("*****************************************\n");
3156
3157 if (gAlice->TreeD())
3158 {
3159 gAlice->TreeD()->GetEvent(0);
3160
3161 Float_t occ[7];
3162 Float_t sum=0;
3163 Float_t mean=0;
3164 printf("\n*****************************************\n");
3165 printf("* Chamber * Digits * Occupancy *\n");
3166 printf("*****************************************\n");
3167
3168 for (Int_t ich=0;ich<7;ich++)
3169 {
3170 TClonesArray *Digits = DigitsAddress(ich); // Raw clusters branch
3171 Int_t ndigits = Digits->GetEntriesFast();
3172 occ[ich] = Float_t(ndigits)/(160*144);
3173 sum += Float_t(ndigits)/(160*144);
3174 printf("* %d * %d * %3.1f%% *\n",ich,ndigits,occ[ich]*100);
3175 }
3176 mean = sum/7;
3177 printf("*****************************************\n");
3178 printf("* Mean occupancy * %3.1f%% *\n",mean*100);
3179 printf("*****************************************\n");
3180 }
3181
3182 printf("\nEnd of analysis\n");
3183
3184}
3185
3186//_________________________________________________________________________________________________
3187
3188
3189void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
3190{
3191
3192AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
3193 AliRICHSegmentationV0* segmentation;
3194 AliRICHChamber* chamber;
3195
3196 chamber = &(pRICH->Chamber(0));
3197 segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel(0);
3198
3199 Int_t NpadX = segmentation->Npx(); // number of pads on X
3200 Int_t NpadY = segmentation->Npy(); // number of pads on Y
3201
3202 //Int_t Pad[144][160];
3203 /*for (Int_t i=0;i<NpadX;i++) {
3204 for (Int_t j=0;j<NpadY;j++) {
3205 Pad[i][j]=0;
3206 }
3207 } */
3208
3209
3210 Int_t xmin= -NpadX/2;
3211 Int_t xmax= NpadX/2;
3212 Int_t ymin= -NpadY/2;
3213 Int_t ymax= NpadY/2;
3214
3215 TH2F *feedback = 0;
3216 TH2F *mip = 0;
3217 TH2F *cerenkov = 0;
3218 TH2F *h = 0;
3219 TH1F *hitsX = 0;
3220 TH1F *hitsY = 0;
3221
3222 TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-30,30,150,-50,10);
3223
3224 if (diaglevel == 1)
3225 {
3226 printf("Single Ring Hits\n");
3227 feedback = new TH2F("feedback","Feedback hit distribution",150,-30,30,150,-50,10);
3228 mip = new TH2F("mip","Mip hit distribution",150,-30,30,150,-50,10);
3229 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-30,30,150,-50,10);
3230 h = new TH2F("h","Detector hit distribution",150,-30,30,150,-50,10);
3231 hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-30,30);
3232 hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,10);
3233 }
3234 else
3235 {
3236 printf("Full Event Hits\n");
3237
3238 feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
3239 mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
3240 cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
3241 h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
3242 hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
3243 hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
3244 }
3245
3246
3247
3248 TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3249 TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3250 TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3251 TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3252 TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3253 TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3254 TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3255
3256 TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
3257 TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",200,.3,1);
3258 TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
3259 TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
3260 TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
3261 TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
3262 TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
3263 TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
3264 TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
3265 //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
3266 TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
3267 TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
3268 TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
3269 TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
3270 TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
3271 TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
3272 TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
3273 TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
3274 TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
3275 TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
3276 TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
3277 TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,0,360);
3278 TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence",100,0,15);
3279 TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.5,1);
3280 TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",200,0,15);
3281 TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",200,0,360);
3282 TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",200,.3,1);
3283 TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",200,.3,1);
3284 TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
3285 TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
3286
3287// Start loop over events
3288
3289 Int_t Nh=0;
3290 Int_t pads=0;
3291 Int_t Nh1=0;
3292 Int_t mothers[80000];
3293 Int_t mothers2[80000];
3294 Float_t mom[3];
3295 Int_t nraw=0;
3296 Int_t phot=0;
3297 Int_t feed=0;
3298 Int_t padmip=0;
3299 Float_t x=0,y=0;
3300
3301 for (Int_t i=0;i<100;i++) mothers[i]=0;
3302
3303 for (int nev=0; nev<= evNumber2; nev++) {
3304 Int_t nparticles = gAlice->GetEvent(nev);
3305
3306
3307 //cout<<"nev "<<nev<<endl;
3308 printf ("\n**********************************\nProcessing Event: %d\n",nev);
3309 //cout<<"nparticles "<<nparticles<<endl;
3310 printf ("Particles : %d\n\n",nparticles);
3311 if (nev < evNumber1) continue;
3312 if (nparticles <= 0) return;
3313
3314// Get pointers to RICH detector and Hits containers
3315
3316
3317 TTree *TH = gAlice->TreeH();
3318 Stat_t ntracks = TH->GetEntries();
3319
3320 // Start loop on tracks in the hits containers
3321 //Int_t Nc=0;
3322 for (Int_t track=0; track<ntracks;track++) {
3323
3324 printf ("\nProcessing Track: %d\n",track);
3325 gAlice->ResetHits();
3326 TH->GetEvent(track);
3327 Int_t nhits = pRICH->Hits()->GetEntriesFast();
3328 if (nhits) Nh+=nhits;
3329 printf("Hits : %d\n",nhits);
3330 for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
3331 mHit;
3332 mHit=(AliRICHHit*)pRICH->NextHit())
3333 {
3334 //Int_t nch = mHit->fChamber; // chamber number
3335 x = mHit->X(); // x-pos of hit
3336 y = mHit->Z(); // y-pos
3337 Float_t phi = mHit->Phi(); //Phi angle of incidence
3338 Float_t theta = mHit->Theta(); //Theta angle of incidence
3339 Int_t index = mHit->Track();
3340 Int_t particle = (Int_t)(mHit->Particle());
3341 //Int_t freon = (Int_t)(mHit->fLoss);
3342
3343 hitsX->Fill(x,(float) 1);
3344 hitsY->Fill(y,(float) 1);
3345
3346 //printf("Particle:%9d\n",particle);
3347
3348 TParticle *current = (TParticle*)gAlice->Particle(index);
3349 //printf("Particle type: %d\n",sizeoff(Particles));
3350
3351 hitsTheta->Fill(theta,(float) 1);
3352 //hitsPhi->Fill(phi,(float) 1);
3353 //if (pRICH->GetDebugLevel() == -1)
3354 //printf("Theta:%f, Phi:%f\n",theta,phi);
3355
3356 //printf("Debug Level:%d\n",pRICH->GetDebugLevel());
3357
3358 if (current->GetPdgCode() < 10000000)
3359 {
3360 mip->Fill(x,y,(float) 1);
3361 //printf("adding mip\n");
3362 //if (current->Energy() - current->GetCalcMass()>1 && freon==1)
3363 //{
3364 hitsPhi->Fill(TMath::Abs(phi),(float) 1);
3365 //hitsTheta->Fill(theta,(float) 1);
3366 //printf("Theta:%f, Phi:%f\n",theta,phi);
3367 //}
3368 }
3369
3370 if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
3371 {
3372 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3373 }
3374 if (TMath::Abs(particle)==2212)
3375 {
3376 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3377 }
3378 if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
3379 || TMath::Abs(particle)==311)
3380 {
3381 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3382 }
3383 if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3384 {
3385 if (current->Energy() - current->GetCalcMass()>1)
3386 chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3387 }
3388 //printf("Hits:%d\n",hit);
3389 //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3390 // Fill the histograms
3391 Nh1+=nhits;
3392 h->Fill(x,y,(float) 1);
3393 //}
3394 //}
3395 }
3396
3397 Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
3398 //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
3399 //totalphotonsevent->Fill(ncerenkovs,(float) 1);
3400
3401 if (ncerenkovs) {
3402 printf("Cerenkovs : %d\n",ncerenkovs);
3403 totalphotonsevent->Fill(ncerenkovs,(float) 1);
3404 for (Int_t hit=0;hit<ncerenkovs;hit++) {
3405 AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
3406 //Int_t nchamber = cHit->fChamber; // chamber number
3407 Int_t index = cHit->Track();
3408 //Int_t pindex = (Int_t)(cHit->fIndex);
3409 Float_t cx = cHit->X(); // x-position
3410 Float_t cy = cHit->Z(); // y-position
3411 Int_t cmother = cHit->fCMother; // Index of mother particle
3412 Int_t closs = (Int_t)(cHit->fLoss); // How did the particle get lost?
3413 Float_t cherenkov = cHit->fCerenkovAngle; //production cerenkov angle
3414 //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
3415
3416 //printf("Particle:%9d\n",index);
3417
3418 TParticle *current = (TParticle*)gAlice->Particle(index);
3419 Float_t energyckov = current->Energy();
3420
3421 if (current->GetPdgCode() == 50000051)
3422 {
3423 if (closs==4)
3424 {
3425 feedback->Fill(cx,cy,(float) 1);
3426 feed++;
3427 }
3428 }
3429 if (current->GetPdgCode() == 50000050)
3430 {
3431
3432 if (closs !=4)
3433 {
3434 phspectra2->Fill(energyckov*1e9,(float) 1);
3435 }
3436
3437 if (closs==4)
3438 {
3439 cerenkov->Fill(cx,cy,(float) 1);
3440
3441 //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
3442
3443 //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
3444 AliRICHHit* mipHit = (AliRICHHit*) pRICH->Hits()->UncheckedAt(0);
3445 mom[0] = current->Px();
3446 mom[1] = current->Py();
3447 mom[2] = current->Pz();
3448 //mom[0] = cHit->fMomX;
3449 // mom[1] = cHit->fMomZ;
3450 //mom[2] = cHit->fMomY;
3451 //Float_t energymip = MIP->Energy();
3452 //Float_t Mip_px = mipHit->fMomFreoX;
3453 //Float_t Mip_py = mipHit->fMomFreoY;
3454 //Float_t Mip_pz = mipHit->fMomFreoZ;
3455 //Float_t Mip_px = MIP->Px();
3456 //Float_t Mip_py = MIP->Py();
3457 //Float_t Mip_pz = MIP->Pz();
3458
3459
3460
3461 //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
3462 //Float_t rt = TMath::Sqrt(r);
3463 //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
3464 //Float_t Mip_rt = TMath::Sqrt(Mip_r);
3465 //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
3466 //Float_t cherenkov = TMath::ACos(coscerenkov);
3467 ckovangle->Fill(cherenkov,(float) 1); //Cerenkov angle calculus
3468 //printf("Cherenkov: %f\n",cherenkov);
3469 Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
3470 hckphi->Fill(ckphi,(float) 1);
3471
3472
3473 //Float_t mix = MIP->Vx();
3474 //Float_t miy = MIP->Vy();
3475 Float_t mx = mipHit->X();
3476 Float_t my = mipHit->Z();
3477 //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
3478 Float_t dx = cx - mx;
3479 Float_t dy = cy - my;
3480 //printf("Dx:%f, Dy:%f\n",dx,dy);
3481 Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
3482 //printf("Final radius:%f\n",final_radius);
3483 radius->Fill(final_radius,(float) 1);
3484
3485 phspectra1->Fill(energyckov*1e9,(float) 1);
3486 phot++;
3487 }
3488 for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
3489 if (cmother == nmothers){
3490 if (closs == 4)
3491 mothers2[cmother]++;
3492 mothers[cmother]++;
3493 }
3494 }
3495 }
3496 }
3497 }
3498
3499
3500 if(gAlice->TreeR())
3501 {
3502 Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
3503 gAlice->TreeR()->GetEvent(nent-1);
3504 TClonesArray *Rawclusters = pRICH->RawClustAddress(2); // Raw clusters branch
3505 //printf ("Rawclusters:%p",Rawclusters);
3506 Int_t nrawclusters = Rawclusters->GetEntriesFast();
3507
3508 if (nrawclusters) {
3509 printf("Raw Clusters : %d\n",nrawclusters);
3510 for (Int_t hit=0;hit<nrawclusters;hit++) {
3511 AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
3512 //Int_t nchamber = rcHit->fChamber; // chamber number
3513 //Int_t nhit = cHit->fHitNumber; // hit number
3514 Int_t qtot = rcHit->fQ; // charge
3515 Float_t fx = rcHit->fX; // x-position
3516 Float_t fy = rcHit->fY; // y-position
3517 //Int_t type = rcHit->fCtype; // cluster type ?
3518 Int_t mult = rcHit->fMultiplicity; // How many pads form the cluster
3519 pads += mult;
3520 if (qtot > 0) {
3521 //printf ("fx: %d, fy: %d\n",fx,fy);
3522 if (fx>(x-4) && fx<(x+4) && fy>(y-4) && fy<(y+4)) {
3523 //printf("There %d \n",mult);
3524 padmip+=mult;
3525 } else {
3526 padnumber->Fill(mult,(float) 1);
3527 nraw++;
3528 if (mult<4) Clcharge->Fill(qtot,(float) 1);
3529 }
3530
3531 }
3532 }
3533 }
3534
3535
3536 TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3537 Int_t nrechits1D = RecHits1D->GetEntriesFast();
3538 //printf (" nrechits:%d\n",nrechits);
3539
3540 if(nrechits1D)
3541 {
3542 for (Int_t hit=0;hit<nrechits1D;hit++) {
3543 AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
3544 Float_t r_omega = recHit1D->fOmega; // Cerenkov angle
3545 Float_t *cer_pho = recHit1D->fCerPerPhoton; // Cerenkov angle per photon
3546 Int_t *padsx = recHit1D->fPadsUsedX; // Pads Used fo reconstruction (x)
3547 Int_t *padsy = recHit1D->fPadsUsedY; // Pads Used fo reconstruction (y)
3548 Int_t goodPhotons = recHit1D->fGoodPhotons; // Number of pads used for reconstruction
3549
3550 Omega1D->Fill(r_omega,(float) 1);
3551
3552 for (Int_t i=0; i<goodPhotons; i++)
3553 {
3554 PhotonCer->Fill(cer_pho[i],(float) 1);
3555 PadsUsed->Fill(padsx[i],padsy[i],1);
3556 //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
3557 }
3558
3559 //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
3560 }
3561 }
3562
3563
3564 TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3565 Int_t nrechits3D = RecHits3D->GetEntriesFast();
3566 //printf (" nrechits:%d\n",nrechits);
3567
3568 if(nrechits3D)
3569 {
3570 for (Int_t hit=0;hit<nrechits3D;hit++) {
3571 AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(hit);
3572 Float_t r_omega = recHit3D->fOmega; // Cerenkov angle
3573 Float_t r_theta = recHit3D->fTheta; // Theta angle of incidence
3574 Float_t r_phi = recHit3D->fPhi; // Phi angle if incidence
3575 Float_t meanradius = recHit3D->fMeanRadius; // Mean radius for reconstructed point
3576
3577 //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY);
3578
3579 Omega3D->Fill(r_omega,(float) 1);
3580 Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
3581 Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
3582 MeanRadius->Fill(meanradius,(float) 1);
3583 }
3584 }
3585 }
3586 }
3587
3588 for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
3589 totalphotonstrack->Fill(mothers[nmothers],(float) 1);
3590 mother->Fill(mothers2[nmothers],(float) 1);
3591 //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
3592 }
3593
3594 clusev->Fill(nraw,(float) 1);
3595 photev->Fill(phot,(float) 1);
3596 feedev->Fill(feed,(float) 1);
3597 padsmip->Fill(padmip,(float) 1);
3598 padscl->Fill(pads,(float) 1);
3599 //printf("Photons:%d\n",phot);
3600 phot = 0;
3601 feed = 0;
3602 pads = 0;
3603 nraw=0;
3604 padmip=0;
3605
3606
3607
3608 gAlice->ResetDigits();
3609 //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
3610 gAlice->TreeD()->GetEvent(0);
3611
3612 if (diaglevel < 4)
3613 {
3614
3615
3616 TClonesArray *Digits = pRICH->DigitsAddress(2);
3617 Int_t ndigits = Digits->GetEntriesFast();
3618 printf("Digits : %d\n",ndigits);
3619 padsev->Fill(ndigits,(float) 1);
3620 for (Int_t hit=0;hit<ndigits;hit++) {
3621 AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3622 Int_t qtot = dHit->Signal(); // charge
3623 Int_t ipx = dHit->PadX(); // pad number on X
3624 Int_t ipy = dHit->PadY(); // pad number on Y
3625 //printf("%d, %d\n",ipx,ipy);
3626 if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
3627 }
3628 }
3629
3630 if (diaglevel == 5)
3631 {
3632 for (Int_t ich=0;ich<7;ich++)
3633 {
3634 TClonesArray *Digits = pRICH->DigitsAddress(ich); // Raw clusters branch
3635 Int_t ndigits = Digits->GetEntriesFast();
3636 //printf("Digits:%d\n",ndigits);
3637 padsev->Fill(ndigits,(float) 1);
3638 if (ndigits) {
3639 for (Int_t hit=0;hit<ndigits;hit++) {
3640 AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3641 //Int_t nchamber = dHit->GetChamber(); // chamber number
3642 //Int_t nhit = dHit->fHitNumber; // hit number
3643 Int_t qtot = dHit->Signal(); // charge
3644 Int_t ipx = dHit->PadX(); // pad number on X
3645 Int_t ipy = dHit->PadY(); // pad number on Y
3646 //Int_t iqpad = dHit->fQpad; // charge per pad
3647 //Int_t rpad = dHit->fRSec; // R-position of pad
3648 //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
3649 if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
3650 if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
3651 if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
3652 if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
3653 if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
3654 if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
3655 if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
3656 if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
3657 }
3658 }
3659 }
3660 }
3661 }
3662
3663
3664 //Create canvases, set the view range, show histograms
3665
3666 TCanvas *c1 = 0;
3667 TCanvas *c2 = 0;
3668 TCanvas *c3 = 0;
3669 TCanvas *c4 = 0;
3670 TCanvas *c5 = 0;
3671 TCanvas *c6 = 0;
3672 TCanvas *c7 = 0;
3673 TCanvas *c8 = 0;
3674 TCanvas *c9 = 0;
3675 TCanvas *c10 = 0;
3676 TCanvas *c11 = 0;
3677 TCanvas *c12 = 0;
3678 //TF1* expo = 0;
3679 //TF1* gaus = 0;
3680
3681
3682 TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3683 Int_t nrechits3D = RecHits3D->GetEntriesFast();
3684 TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3685 Int_t nrechits1D = RecHits1D->GetEntriesFast();
3686
3687 switch(diaglevel)
3688 {
3689 case 1:
3690
3691 c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
3692 hc0->SetXTitle("ix (npads)");
3693 hc0->Draw("box");
3694
3695//
3696 c2 = new TCanvas("c2","Hits per type",100,100,600,700);
3697 c2->Divide(2,2);
3698 //c4->SetFillColor(42);
3699
3700 c2->cd(1);
3701 feedback->SetXTitle("x (cm)");
3702 feedback->SetYTitle("y (cm)");
3703 feedback->Draw();
3704
3705 c2->cd(2);
3706 //mip->SetFillColor(5);
3707 mip->SetXTitle("x (cm)");
3708 mip->SetYTitle("y (cm)");
3709 mip->Draw();
3710
3711 c2->cd(3);
3712 //cerenkov->SetFillColor(5);
3713 cerenkov->SetXTitle("x (cm)");
3714 cerenkov->SetYTitle("y (cm)");
3715 cerenkov->Draw();
3716
3717 c2->cd(4);
3718 //h->SetFillColor(5);
3719 h->SetXTitle("x (cm)");
3720 h->SetYTitle("y (cm)");
3721 h->Draw();
3722
3723 c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
3724 c3->Divide(2,1);
3725 //c10->SetFillColor(42);
3726
3727 c3->cd(1);
3728 hitsX->SetFillColor(5);
3729 hitsX->SetXTitle("(cm)");
3730 hitsX->Draw();
3731
3732 c3->cd(2);
3733 hitsY->SetFillColor(5);
3734 hitsY->SetXTitle("(cm)");
3735 hitsY->Draw();
3736
3737
3738 break;
3739//
3740 case 2:
3741
3742 c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
3743 c4->Divide(2,1);
3744 //c6->SetFillColor(42);
3745
3746 c4->cd(1);
3747 phspectra2->SetFillColor(5);
3748 phspectra2->SetXTitle("energy (eV)");
3749 phspectra2->Draw();
3750 c4->cd(2);
3751 phspectra1->SetFillColor(5);
3752 phspectra1->SetXTitle("energy (eV)");
3753 phspectra1->Draw();
3754
3755 c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
3756 c5->Divide(2,2);
3757 //c9->SetFillColor(42);
3758
3759 c5->cd(1);
3760 pionspectra->SetFillColor(5);
3761 pionspectra->SetXTitle("(GeV)");
3762 pionspectra->Draw();
3763
3764 c5->cd(2);
3765 protonspectra->SetFillColor(5);
3766 protonspectra->SetXTitle("(GeV)");
3767 protonspectra->Draw();
3768
3769 c5->cd(3);
3770 kaonspectra->SetFillColor(5);
3771 kaonspectra->SetXTitle("(GeV)");
3772 kaonspectra->Draw();
3773
3774 c5->cd(4);
3775 chargedspectra->SetFillColor(5);
3776 chargedspectra->SetXTitle("(GeV)");
3777 chargedspectra->Draw();
3778
3779 break;
3780
3781 case 3:
3782
3783
3784 if(gAlice->TreeR())
3785 {
3786 c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
3787 c6->Divide(2,2);
3788 //c3->SetFillColor(42);
3789
3790 c6->cd(1);
3791 //TPad* c6_1;
3792 //c6_1->SetLogy();
3793 Clcharge->SetFillColor(5);
3794 Clcharge->SetXTitle("ADC counts");
3795 if (evNumber2>10)
3796 {
3797 Clcharge->Fit("expo");
3798 //expo->SetLineColor(2);
3799 //expo->SetLineWidth(3);
3800 }
3801 Clcharge->Draw();
3802
3803 c6->cd(2);
3804 padnumber->SetFillColor(5);
3805 padnumber->SetXTitle("(counts)");
3806 padnumber->Draw();
3807
3808 c6->cd(3);
3809 clusev->SetFillColor(5);
3810 clusev->SetXTitle("(counts)");
3811 if (evNumber2>10)
3812 {
3813 clusev->Fit("gaus");
3814 //gaus->SetLineColor(2);
3815 //gaus->SetLineWidth(3);
3816 }
3817 clusev->Draw();
3818
3819 c6->cd(4);
3820 padsmip->SetFillColor(5);
3821 padsmip->SetXTitle("(counts)");
3822 padsmip->Draw();
3823 }
3824
3825 if(evNumber2<1)
3826 {
3827 c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
3828 mother->SetFillColor(5);
3829 mother->SetXTitle("counts");
3830 mother->Draw();
3831 }
3832
3833 c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
3834 c7->Divide(2,2);
3835 //c7->SetFillColor(42);
3836
3837 c7->cd(1);
3838 totalphotonsevent->SetFillColor(5);
3839 totalphotonsevent->SetXTitle("Photons (counts)");
3840 if (evNumber2>10)
3841 {
3842 totalphotonsevent->Fit("gaus");
3843 //gaus->SetLineColor(2);
3844 //gaus->SetLineWidth(3);
3845 }
3846 totalphotonsevent->Draw();
3847
3848 c7->cd(2);
3849 photev->SetFillColor(5);
3850 photev->SetXTitle("(counts)");
3851 if (evNumber2>10)
3852 {
3853 photev->Fit("gaus");
3854 //gaus->SetLineColor(2);
3855 //gaus->SetLineWidth(3);
3856 }
3857 photev->Draw();
3858
3859 c7->cd(3);
3860 feedev->SetFillColor(5);
3861 feedev->SetXTitle("(counts)");
3862 if (evNumber2>10)
3863 {
3864 feedev->Fit("gaus");
3865 //gaus->SetLineColor(2);
3866 //gaus->SetLineWidth(3);
3867 }
3868 feedev->Draw();
3869
3870 c7->cd(4);
3871 padsev->SetFillColor(5);
3872 padsev->SetXTitle("(counts)");
3873 if (evNumber2>10)
3874 {
3875 padsev->Fit("gaus");
3876 //gaus->SetLineColor(2);
3877 //gaus->SetLineWidth(3);
3878 }
3879 padsev->Draw();
3880
3881 break;
3882
3883 case 4:
3884
3885
3886 if(nrechits3D)
3887 {
3888 c8 = new TCanvas("c8","3D reconstruction",50,50,1100,700);
3889 c8->Divide(4,2);
3890 //c2->SetFillColor(42);
3891
3892 c8->cd(1);
3893 hitsPhi->SetFillColor(5);
3894 hitsPhi->Draw();
3895 c8->cd(2);
3896 hitsTheta->SetFillColor(5);
3897 hitsTheta->Draw();
3898 c8->cd(3);
3899 ckovangle->SetFillColor(5);
3900 ckovangle->SetXTitle("angle (radians)");
3901 ckovangle->Draw();
3902 c8->cd(4);
3903 radius->SetFillColor(5);
3904 radius->SetXTitle("radius (cm)");
3905 radius->Draw();
3906 c8->cd(5);
3907 Phi->SetFillColor(5);
3908 Phi->Draw();
3909 c8->cd(6);
3910 Theta->SetFillColor(5);
3911 Theta->Draw();
3912 c8->cd(7);
3913 Omega3D->SetFillColor(5);
3914 Omega3D->SetXTitle("angle (radians)");
3915 Omega3D->Draw();
3916 c8->cd(8);
3917 MeanRadius->SetFillColor(5);
3918 MeanRadius->SetXTitle("radius (cm)");
3919 MeanRadius->Draw();
3920
3921 }
3922
3923 if(nrechits1D)
3924 {
3925 c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
3926 c9->Divide(3,2);
3927 //c5->SetFillColor(42);
3928
3929 c9->cd(1);
3930 ckovangle->SetFillColor(5);
3931 ckovangle->SetXTitle("angle (radians)");
3932 ckovangle->Draw();
3933
3934 c9->cd(2);
3935 radius->SetFillColor(5);
3936 radius->SetXTitle("radius (cm)");
3937 radius->Draw();
3938
3939 c9->cd(3);
3940 hc0->SetXTitle("pads");
3941 hc0->Draw("box");
3942
3943 c9->cd(5);
3944 Omega1D->SetFillColor(5);
3945 Omega1D->SetXTitle("angle (radians)");
3946 Omega1D->Draw();
3947
3948 c9->cd(4);
3949 PhotonCer->SetFillColor(5);
3950 PhotonCer->SetXTitle("angle (radians)");
3951 PhotonCer->Draw();
3952
3953 c9->cd(6);
3954 PadsUsed->SetXTitle("pads");
3955 PadsUsed->Draw("box");
3956 }
3957
3958 break;
3959
3960 case 5:
3961
3962 printf("Drawing histograms.../n");
3963
3964 //if (ndigits)
3965 //{
3966 c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
3967 c1->Divide(4,2);
3968 //c1->SetFillColor(42);
3969
3970 c10->cd(1);
3971 hc1->SetXTitle("ix (npads)");
3972 hc1->Draw("box");
3973 c10->cd(2);
3974 hc2->SetXTitle("ix (npads)");
3975 hc2->Draw("box");
3976 c10->cd(3);
3977 hc3->SetXTitle("ix (npads)");
3978 hc3->Draw("box");
3979 c10->cd(4);
3980 hc4->SetXTitle("ix (npads)");
3981 hc4->Draw("box");
3982 c10->cd(5);
3983 hc5->SetXTitle("ix (npads)");
3984 hc5->Draw("box");
3985 c10->cd(6);
3986 hc6->SetXTitle("ix (npads)");
3987 hc6->Draw("box");
3988 c10->cd(7);
3989 hc7->SetXTitle("ix (npads)");
3990 hc7->Draw("box");
3991 c10->cd(8);
3992 hc0->SetXTitle("ix (npads)");
3993 hc0->Draw("box");
3994 //}
3995//
3996 c11 = new TCanvas("c11","Hits per type",100,100,600,700);
3997 c11->Divide(2,2);
3998 //c4->SetFillColor(42);
3999
4000 c11->cd(1);
4001 feedback->SetXTitle("x (cm)");
4002 feedback->SetYTitle("y (cm)");
4003 feedback->Draw();
4004
4005 c11->cd(2);
4006 //mip->SetFillColor(5);
4007 mip->SetXTitle("x (cm)");
4008 mip->SetYTitle("y (cm)");
4009 mip->Draw();
4010
4011 c11->cd(3);
4012 //cerenkov->SetFillColor(5);
4013 cerenkov->SetXTitle("x (cm)");
4014 cerenkov->SetYTitle("y (cm)");
4015 cerenkov->Draw();
4016
4017 c11->cd(4);
4018 //h->SetFillColor(5);
4019 h->SetXTitle("x (cm)");
4020 h->SetYTitle("y (cm)");
4021 h->Draw();
4022
4023 c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
4024 c12->Divide(2,1);
4025 //c10->SetFillColor(42);
4026
4027 c12->cd(1);
4028 hitsX->SetFillColor(5);
4029 hitsX->SetXTitle("(cm)");
4030 hitsX->Draw();
4031
4032 c12->cd(2);
4033 hitsY->SetFillColor(5);
4034 hitsY->SetXTitle("(cm)");
4035 hitsY->Draw();
4036
4037 break;
4038
4039 }
4040
4041
4042 // calculate the number of pads which give a signal
4043
4044
4045 //Int_t Np=0;
4046 /*for (Int_t i=0;i< NpadX;i++) {
4047 for (Int_t j=0;j< NpadY;j++) {
4048 if (Pad[i][j]>=6){
4049 Np+=1;
4050 }
4051 }
4052 }*/
4053 //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
4054 printf("\nEnd of analysis\n");
4055 printf("**********************************\n");
4056}
4057
4058////////////////////////////////////////////////////////////////////////
4059void AliRICH::MakeBranchInTreeD(TTree *treeD, const char *file)
4060{
4061 //
4062 // Create TreeD branches for the RICH.
4063 //
4064
4065 const Int_t kBufferSize = 4000;
4066 char branchname[30];
4067
4068 //
4069 // one branch for digits per chamber
4070 //
4071 for (Int_t i=0; i<kNCH ;i++) {
4072 sprintf(branchname,"%sDigits%d",GetName(),i+1);
4073 if (fDchambers && treeD) {
4074 MakeBranchInTree(treeD,
4075 branchname, &((*fDchambers)[i]), kBufferSize, file);
4076// printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
4077 }
4078 }
4079}
4080////////////////////////////////////////////////////////////////////////