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