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