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