Setters fro models moved to header
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
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
199 static Int_t sMaxIterPad=0;
200 static Int_t sCurIterPad=0;
201  
202 ClassImp(AliRICH)
203     
204 //___________________________________________
205 AliRICH::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 //___________________________________________
233 AliRICH::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
301 AliRICH::AliRICH(const AliRICH& RICH)
302 {
303 // Copy Constructor
304 }
305
306
307 //___________________________________________
308 AliRICH::~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 //_____________________________________________________________________________
344 Int_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     
390 return nnew;
391 }
392 //___________________________________________
393 void 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 //___________________________________________
406 void 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 //___________________________________________
432 void AliRICH::SDigits2Digits()
433 {
434   SDigits2Digits(0,0);
435 }
436 //___________________________________________
437 void 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 //___________________________________________
450 void 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 //_____________________________________________________________________________
460 void 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 //___________________________________________
472 void 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 //_____________________________________________________________________________
484 void 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 //_____________________________________________________________________________
498 void 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 //_____________________________________________________________________________
510 void 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 //_____________________________________________________________________________
523 void 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 //___________________________________________
536 void 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 //___________________________________________
793 void 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 //___________________________________________
1289 void 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
1614 Float_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 //__________________________________________
1675 Float_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 //___________________________________________
1821 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1822 {
1823
1824 // Default value
1825
1826     return 9999;
1827 }
1828
1829 //___________________________________________
1830 void 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 //___________________________________________
1927 void 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 //___________________________________________
2005 void 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 //____________________________________________
2017 void 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 //____________________________________________
2030 void 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 //____________________________________________
2043 void 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 //____________________________________________
2057 void 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 //___________________________________________
2072 void 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
2580 void 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
2634 AliRICHSDigit* 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
2651 AliRICHSDigit* 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
2664 AliRICH& AliRICH::operator=(const AliRICH& rhs)
2665 {
2666 // Assignment operator
2667     return *this;
2668     
2669 }
2670
2671 void 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
3189 void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
3190 {
3191
3192 AliRICH *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++) {