Modifications to allow Cherenkov transport
[u/mrichter/AliRoot.git] / RICH / AliRICH.cxx
1 ///////////////////////////////////////////////////////////////////////////////
2 //                                                                           //
3 //  Ring Imaging Cherenkov                                                   //
4 //  This class contains the basic functions for the Ring Imaging Cherenkov   //
5 //  detector. Functions specific to one particular geometry are              //
6 //  contained in the derived classes                                         //
7 //                                                                           //
8 //Begin_Html
9 /*
10 <img src="picts/AliRICHClass.gif">
11 */
12 //End_Html
13 //                                                                           //
14 //                                                                           //
15 ///////////////////////////////////////////////////////////////////////////////
16
17 #include <TBRIK.h> 
18 #include <TNode.h> 
19 #include <TRandom.h> 
20
21 #include <TClass.h> 
22 #include "AliRICH.h"
23 #include "AliRun.h"
24 #include "TGeant3.h" 
25
26
27 ClassImp(AliRICH)
28
29 //_____________________________________________________________________________
30 AliRICH::AliRICH()
31 {
32   //
33   // Default constructor for RICH
34   //
35   fIshunt   = 0;
36   fHits     = 0;
37   fMips     = 0;
38   fCkovs    = 0;
39   fPadhits  = 0;
40   fNmips    = 0;
41   fNckovs   = 0;
42   fNpadhits = 0;
43   
44   fChslope  = 0;
45   fAlphaFeed= 0;
46   fSxcharge = 0;
47   fIritri   = 0;
48 }
49  
50 //_____________________________________________________________________________
51 AliRICH::AliRICH(const char *name, const char *title)
52        : AliDetector(name,title)
53 {
54   //
55   // Standard constructor for RICH
56   //
57   fIshunt     = 0;
58   fNmips      = 0;
59   fNckovs     = 0;
60   fNpadhits   = 0;
61   //
62   // Allocate space for different components
63   fHits     = new TClonesArray("AliRICHhit", 100);
64   fMips     = new TClonesArray("AliRICHmip", 100);
65   fCkovs    = new TClonesArray("AliRICHckov", 100);
66   fPadhits  = new TClonesArray("AliRICHpadhit", 100);
67   //
68   // Set parameters to default value
69   fChslope  = 40;
70   fAlphaFeed= 0.04;
71   fSxcharge = 0.18;
72   fIritri   = 0;
73   //
74   SetMarkerColor(6);
75   SetMarkerStyle(20);
76   SetMarkerSize(0.5);
77 }
78
79 //_____________________________________________________________________________
80 AliRICH::~AliRICH()
81 {
82   //
83   // Destructor for RICH
84   //
85   fIshunt   = 0;
86   delete fHits;
87   fMips->Delete();    delete fMips;
88   fCkovs->Delete();   delete fCkovs;
89   fPadhits->Delete(); delete fPadhits;
90 }
91
92 //_____________________________________________________________________________
93 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
94 {
95   //
96   // Add a RICH hit
97   //
98   switch ((int) hits[0]) {
99   case 0:
100     {
101       //
102       // Simple hit
103       TClonesArray &lhits = *fHits;
104       new(lhits[fNhits++]) AliHit(fIshunt,track);
105       AliHit *lhit = (AliHit*) fHits->AddrAt(fNhits-1);
106       lhit->fX = hits[19];
107       lhit->fY = hits[20];
108       lhit->fZ = hits[21];
109     }   break;
110   case 1:
111     //
112     // MIP hit
113     AddMipHit(track,vol,hits);
114     break;
115   case 2:
116     //
117     // Cherenkov hit
118     AddCkovHit(track,vol,hits);
119     break;
120   case 3:
121     //
122     // Pad hit
123     AddPadHit(track,vol,hits);
124     break;
125   case 4:
126     // Update a mip hit
127     UpdateMipHit(hits);
128     break; 
129   default:
130     printf("Error: AliRICH::AddHit flag %d not defined./n",(int) hits[0]);
131     return;
132   }    
133 }
134
135 //_____________________________________________________________________________
136 void AliRICH::AddMipHit(Int_t track, Int_t *vol, Float_t *hits)
137
138     // Adds a mip hit in the RICH.
139     TClonesArray &lhits = *fMips;
140     new(lhits[fNmips++]) AliRICHmip(fIshunt,track,vol,hits,
141                                     fNckovs,fNpadhits);
142 }
143
144 //_____________________________________________________________________________
145 void AliRICH::AddCkovHit(Int_t track, Int_t *vol, Float_t *hits)
146
147   //
148   // Adds a cerenkov hit in the RICH.
149   //
150   TClonesArray &lhits = *fCkovs;
151   AliRICHmip *lmip = (AliRICHmip*) fMips->AddrAt(fNmips-1);
152   //
153   // If this ckov come from a mip update the mip lastckov entry.
154   Int_t fmipslocal=-1;
155   if (lmip->GetZ() != -999.0) {
156     fmipslocal    =   fNmips-1;
157     lmip->SetLastCkov(fNckovs);
158   }
159   new(lhits[fNckovs++]) AliRICHckov(fIshunt,track,vol,hits,
160                                     fmipslocal,fNpadhits);
161 }
162
163 //_____________________________________________________________________________
164 void AliRICH::AddPadHit(Int_t track, Int_t *vol, Float_t *hits)
165
166   //
167   // Adds pad hits in the RICH
168   //
169   TClonesArray &lhits = *fPadhits;
170   // Update the last padhit of the respective particle:
171   if ((int) hits[1]==50) { // a ckov
172     ((AliRICHckov *) fCkovs->AddrAt(fNckovs-1))->SetLastpad(fNpadhits);
173     new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,-1,fNckovs-1);
174   }else { // a mip
175     ((AliRICHmip *) fMips->AddrAt(fNmips-1))->SetLastpad(fNpadhits);
176     new(lhits[fNpadhits++]) AliRICHpadhit(fIshunt,track,vol,hits,fNmips-1,-1);
177   }
178 }
179
180 //_____________________________________________________________________________
181 void AliRICH::BuildGeometry()
182 {
183   //
184   // Builds a TNode geometry for event display
185   //
186   TNode *Node, *Top;
187   
188   const int kColorRICH = kGreen;
189   //
190   Top=gAlice->GetGeometry()->GetNode("alice");
191
192   new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
193   new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
194   new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
195   new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
196   new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
197   new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
198   new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
199   new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
200   Top->cd();
201   Node = new TNode("RICH1","RICH1","S_RICH",0,471.8999,165.2599,"rot993");
202   Node->SetLineColor(kColorRICH);
203   fNodes->Add(Node);
204   Top->cd();
205   Node = new TNode("RICH2","RICH2","S_RICH",171,470,0,"rot994");
206   Node->SetLineColor(kColorRICH);
207   fNodes->Add(Node);
208   Top->cd();
209   Node = new TNode("RICH3","RICH3","S_RICH",0,500,0,"rot995");
210   Node->SetLineColor(kColorRICH);
211   fNodes->Add(Node);
212   Top->cd();
213   Node = new TNode("RICH4","RICH4","S_RICH",-171,470,0,"rot996");
214   Node->SetLineColor(kColorRICH);
215   fNodes->Add(Node);
216   Top->cd();
217   Node = new TNode("RICH5","RICH5","S_RICH",161.3999,443.3999,-165.3,"rot997");
218   Node->SetLineColor(kColorRICH);
219   fNodes->Add(Node);
220   Top->cd();
221   Node = new TNode("RICH6","RICH6","S_RICH",0,471.8999,-165.3,"rot998");
222   Node->SetLineColor(kColorRICH);
223   fNodes->Add(Node);
224   Top->cd();
225   Node = new TNode("RICH7","RICH7","S_RICH",-161.399,443.3999,-165.3,"rot999");
226   Node->SetLineColor(kColorRICH);
227   fNodes->Add(Node); 
228 }
229
230 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
231 //
232 //                       Section for Rich Step
233 #define MAXPH 1000
234
235 static Float_t sVloc[3];
236 static Float_t sVectIn[3];
237 static Float_t sDetot;
238 static Float_t sYanode[MAXPH];
239 static Float_t sXpad[MAXPH];
240 static Float_t sYpad[MAXPH];
241 static Int_t sNpx;
242 static Int_t sNpy;
243 static Float_t sDxp;
244 static Float_t sDyp;
245 static Float_t sDlx;
246 static Float_t sDly;
247 static Float_t sDpad;
248
249
250 static Int_t sNsecon;
251 //static Float_t sQint[2];
252
253 #define MAXSEC 2000
254
255 static Int_t sNpads;
256 static Int_t sIpx[MAXSEC];
257 static Int_t sIpy[MAXSEC];
258 static Int_t sIqpad[MAXSEC];
259 static Int_t sNphlink[MAXSEC];
260 static Int_t sNphoton;
261
262 static Int_t sNfeed, sNfeedd, sNdir;
263
264 static Float_t sPparent;
265 static Float_t sThincParent;
266 static Int_t sIloss[MAXPH];
267 static Int_t sIprod[MAXPH];
268 static Float_t sXphit[MAXPH];
269 static Float_t sYphit[MAXPH];
270
271 static Float_t sSycharge;
272 static Float_t sXsox, sYsox, sZsox;
273
274 static Int_t sMckov[MAXPH];
275 static Int_t idpartgx;
276 static Float_t phisx;
277 static Int_t nmodsx;
278 static Float_t psx;
279 static Float_t xsx;
280 static Float_t ysx;
281 static Float_t thetasx;
282
283
284 //_____________________________________________________________________________
285 void AliRICH::Init()
286 {
287   //
288   // Initialise RICH after that it has been built
289   //
290   const Float_t sconv=2*TMath::Sqrt(2*TMath::Log(2));
291   const Float_t yK3=1.20;
292   Float_t ansp;
293   Int_t i;
294   //
295   sNpx=162;
296   sNpy=162;
297   sDxp=0.80;
298   sDyp=0.80;
299   ansp=sDyp/2;
300   sDlx=sNpx*sDxp/2;
301   sDly=sNpy*sDyp/2;
302   sDpad=0.2;
303   //
304   for(i=0;i<sNpx;i++) {
305     sXpad[i]=i*sDxp;
306     sYpad[i]=i*sDyp;
307   }
308   for(i=0;i<2*sNpy+1;i++) sYanode[i]=ansp/2+i*ansp;
309   //
310   
311   sSycharge=4*TMath::ATanH(1/TMath::Sqrt(2*yK3))/TMath::Pi()
312     /(1-0.5*TMath::Sqrt(yK3))/sDpad/sconv;
313   //
314 }
315
316
317 //___________________________________________________________________________
318 void AliRICH::StepManager()
319 {
320   //
321   // Called at every step in the RICH
322   //
323
324   TGeant3 *geant3 = (TGeant3*) gMC;
325
326   const Float_t xshift[3] = { 41.3, 0, -41.3 };
327   const Int_t nrooth = 25;
328   
329   static Int_t ixold=-1, iyold=-1;
330   
331   // System generated locals 
332   Int_t j, i1;
333   Float_t r1, r2;
334   
335   // Local variables 
336   Float_t ranf[2], rrhh[nrooth], phiangle, cost, vmod;
337   //Int_t idpartsx;
338   Int_t i;
339   Float_t t, vxloc[3];
340   Int_t ll, irivol[2];
341   Int_t lcase;
342   //Int_t iprimx;
343   Int_t ix, iy;
344   Float_t stwght;
345   Float_t cophi;
346   Float_t dir[3];
347   Int_t ihitrak;
348   Int_t medprod;
349   
350   Int_t nmult=0;
351   //Float_t xtrig[200], ytrig[200];
352   //Int_t itrig[200];
353
354   
355   
356   //     ILOSS = 0    NOT LOST 
357   //             1    REFLECTED FREON-QUARZ 
358   //             2    REFLECTED QUARZ-METHANE 
359   //             3    REFLECTED METHANE-CSI 
360   //            11    ABSORBED IN FREON 
361   //            12    ABSORBED IN QUARZ 
362   //            13    ABSORBED IN METHANE 
363   //            21    CSI QE 
364   //     IPROD = 1    PRODUCED IN FREON 
365   //     IPROD = 2    PRODUCED IN QUARZ 
366   
367   // new (changed NROOTH from 10 to 25!!!!!!!!!!!!!) 
368   
369   
370   Int_t *idtmed = fIdtmed->GetArray()-999;
371   
372   //--------------------------------------------------------------------------
373   
374   //        MIP inside CsI 
375
376   if (geant3->Gckine()->charge) {
377     
378     //        Charged particles treatment 
379     if (fIritri && !geant3->Gctrak()->upwght) {
380       if (geant3->Gctmed()->numed == idtmed[fIritri-1]) {
381         if (geant3->Gcking()->ngkine > 0) {
382           strncpy((char *)&lcase,"HADR",4);
383           if (geant3->Gcking()->kcase == lcase) {
384             i1 = geant3->Gcking()->ngkine;
385             for (i = 1; i <= i1; ++i) {
386               gMC->Gmtod(geant3->Gckin3()->gpos[i-1], vxloc, 1);
387               gMC->Gmtod(geant3->Gcking()->gkin[i-1], dir, 2);
388               if (geant3->Gcking()->gkin[i-1][4] == 8. || 
389                   geant3->Gcking()->gkin[i-1][4] == 9.) {
390                 ++nmult;
391                 // Computing 2nd power 
392                 r1 = dir[0];
393                 // Computing 2nd power 
394                 r2 = dir[2];
395                 //theta = TMath::ATan2(TMath::Sqrt(r1*r1+r2*r2),dir[1]);
396                 //xtrig[nmult - 1] = theta;
397                 //ytrig[nmult - 1] = vxloc[1] + .25;
398                 //itrig[nmult - 1] = (Int_t) geant3->Gcking()->gkin[i-1][4];
399               }
400             }
401           }
402         }
403       }
404       if ((geant3->Gctmed()->numed == idtmed[1006-1] &&
405            geant3->Gctrak()->inwvol == 2) || 
406           geant3->Gctrak()->istop) {
407         if (!nmult) {
408           printf("NOT TRIGGERED\n");
409           sDetot = 0.;
410           sNsecon = 0;
411           sNpads = 0;
412           sNphoton = 0;
413           sNfeed = 0;
414           sNfeedd = 0;
415           sNdir = 0;
416           geant3->Gctrak()->istory = 0;
417           geant3->Gctrak()->upwght = 0.;
418           geant3->Gcflag()->ieotri = 1;
419           nmult = 0;
420           //sQint[0] = 0.;
421           //sQint[1] = 0.;
422         } else {
423           printf("TRIGGERED %d\n",nmult);
424         }
425       }
426     }
427     //        MIP inside Methane 
428     if (geant3->Gctmed()->numed == idtmed[1009-1]) {
429       
430       // new 
431       //     If particle produced already Cerenkov Photons (istory=1) 
432       //     update the impact point only 
433       if (geant3->Gctrak()->istory == 1) {
434         //     Direction of incidence and where did it hit ? 
435         gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
436         gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
437         phiangle = TMath::ATan2(dir[2], dir[0]);
438         if (phiangle < 0.) phiangle += 2*TMath::Pi();
439         i1 = nrooth;
440         for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
441         irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
442         irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
443         // NMODSX 
444         ihitrak = gAlice->CurrentTrack();
445         rrhh[0] = 4.;
446         // flag to say this is update 
447         rrhh[1] = sVloc[0] + sDlx;
448         // XSX 
449         rrhh[2] = sVloc[2] + sDly;
450         // YSX 
451         rrhh[3] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
452         // NMODSX 
453         // Computing 2nd power 
454         r1 = dir[0];
455         // Computing 2nd power 
456         r2 = dir[2];
457         rrhh[4] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
458         // theta 
459         rrhh[5] = phiangle;
460         //          PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK 
461         AddHit(ihitrak,irivol,rrhh);
462       }
463       // enew 
464       //        Record particle properties 
465       //        If particle produced already Cerenkov Photons (istory=1)
466       
467       //        update the impact point only 
468       if (geant3->Gctrak()->istory != 2) {
469         if (!geant3->Gctrak()->istory) {
470           ++sNsecon;
471           
472           //        Is this a primary particle ? 
473           //iprimx = 1;
474           //if (geant3->Gctrak()->upwght) iprimx = 0;
475           
476           //        Where did it come from ? 
477           sXsox = geant3->Gckine()->vert[0];
478           sYsox = geant3->Gckine()->vert[1];
479           sZsox = geant3->Gckine()->vert[2];
480           
481           //        Momentum 
482           psx = geant3->Gctrak()->vect[6];
483           
484           //        Particle type and parent 
485           //idpartsx = geant3->Gckine()->ipart;
486           r1 = geant3->Gctrak()->upwght / 100.;
487           idpartgx = Int_t(r1+0.5);
488           if (!geant3->Gctrak()->upwght) {
489             sPparent = geant3->Gctrak()->vect[6];
490             sThincParent = thetasx;
491           }
492           
493           //        Direction of incidence and where did it hit ? 
494           gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
495           gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
496           // Computing 2nd power 
497           r1 = dir[0];
498           // Computing 2nd power 
499           r2 = dir[2];
500           thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
501           phisx = TMath::ATan2(dir[2], dir[0]);
502           if (phisx < 0.) phisx += 2*TMath::Pi();
503           ysx = sVloc[2] + sDly;
504           xsx = sVloc[0] + sDlx;
505           nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
506           //     new 
507           i1 = nrooth;
508           for (ll = 0; ll < i1; ++ll) rrhh[ll] = 0;
509           irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
510           irivol[1] = nmodsx;
511           ihitrak = gAlice->CurrentTrack();
512           rrhh[0] = 1.;
513           // Flag to say this is MIP 
514           rrhh[1] = (Float_t) geant3->Gckine()->ipart;
515           rrhh[2] = xsx;
516           rrhh[3] = ysx;
517           rrhh[4] = (Float_t) nmodsx;
518           // Module Number 
519           rrhh[5] = thetasx;
520           rrhh[6] = geant3->Gctrak()->tofg;
521           // in seconds 
522           rrhh[7] = (Float_t) idpartgx;
523           //     mips specific 
524           rrhh[8] = phisx;
525           rrhh[9] = psx;
526           //     charge of current particle in electron charge unit; 
527           rrhh[10] = geant3->Gckine()->charge;
528           rrhh[11] = -999.;
529           // Zo of ckov generation 
530           rrhh[12] = 0.;
531           // no ckov !!! 
532           AddHit(ihitrak, irivol, rrhh);
533           //     end of new 
534           
535           //        Earmark track as being recorded in methane gap 
536           
537           geant3->Gctrak()->istory = 2;
538         }
539       }
540       
541       //        Signal generation in methane gap 
542       gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
543       gMC->Gmtod(geant3->Gckine()->vert, vxloc, 1);
544       ix = (Int_t) ((sVloc[0] + sDlx) /  sDxp);
545       iy = (Int_t) ((sVloc[2] + sDly) /  sDyp);
546       
547       //        Is this the first step? 
548       if (ixold == -1 && iyold == -1) {
549         ixold = ix;
550         iyold = iy;
551         for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
552       }
553       
554       //        Mip left gap 
555       if (geant3->Gctrak()->inwvol == 2 || geant3->Gctrak()->istop) {
556         sDetot += geant3->Gctrak()->destep;
557         if (sDetot > 0.) RichIntegration();
558         sDetot = 0.;
559         ixold = -1;
560         iyold = -1;
561         
562         //        Mip left current pad 
563       } else if (ixold != ix || iyold != iy) {
564         if (sDetot > 0.) RichIntegration();
565         for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
566         sDetot = geant3->Gctrak()->destep;
567         ixold = ix;
568         iyold = iy;
569       } else {
570         sDetot += geant3->Gctrak()->destep;
571       }
572     }
573   }
574   
575   //        End charged particles treatment 
576   //
577   //        Treat photons produced in Freon and Quartz 
578   if (geant3->Gckin2()->ngphot > 0 && 
579       (geant3->Gctmed()->numed == idtmed[1004-1] || 
580        geant3->Gctmed()->numed == idtmed[1003-1])) {
581     if (!geant3->Gctrak()->upwght) {
582       
583       //        If it is a primary, save all generated photons 
584       i1 = geant3->Gckin2()->ngphot;
585       for (i = 1; i <= i1; ++i) {
586         ++sNphoton;
587         if (sNphoton > MAXPH) {
588           sNphoton = MAXPH;
589           printf("ATTENTION NPHOTON %d\n",sNphoton);
590           continue;
591         }
592         
593         //        Production medium 
594         medprod = 1;
595         if (geant3->Gctmed()->numed == idtmed[1003-1]) medprod = 2;
596         //
597         //        Production angle and energy 
598         vmod=0;
599         cost=0;
600         for(j=0;j<3;j++) {
601           cost+=geant3->Gckin2()->xphot[i-1][3+j]*geant3->Gctrak()->vect[3+j];
602           vmod+=geant3->Gckin2()->xphot[i-1][3+j]*
603             geant3->Gckin2()->xphot[i-1][3+j];
604         }
605         cost/=sqrt(vmod);
606         sIloss[sNphoton - 1] = 22;
607         sIprod[sNphoton - 1] = medprod;
608         sXphit[sNphoton - 1] = 0.;
609         sYphit[sNphoton - 1] = 0.;
610         stwght = geant3->Gctrak()->upwght;
611         geant3->Gctrak()->upwght = (Float_t) sNphoton;
612         //      geant3->Gskpho(i);
613         sMckov[sNphoton - 1] = gAlice->CurrentTrack()+i;
614         geant3->Gctrak()->upwght = stwght;
615       }
616     } else {
617       stwght = geant3->Gctrak()->upwght;
618       geant3->Gctrak()->upwght = 0.;
619       //      geant3->Gskpho(0);
620       geant3->Gctrak()->upwght = stwght;
621     }
622     
623     //        Particle did not yet pass the methane gap 
624     if (geant3->Gctrak()->istory == 0) {
625       geant3->Gctrak()->istory = 1;
626       ++sNsecon;
627       //        Is this a primary particle ? 
628       //iprimx = 1;
629       //if (geant3->Gctrak()->upwght) iprimx = 0;
630       
631       //        Where did it come from ? 
632       sXsox = geant3->Gckine()->vert[0];
633       sYsox = geant3->Gckine()->vert[1];
634       sZsox = geant3->Gckine()->vert[2];
635       
636       //        Where did it hit ? 
637       gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
638       gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
639       ysx = sVloc[2] + sDly;
640       if (geant3->Gctmed()->numed == idtmed[1004-1]) {
641         nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
642         xsx = sVloc[0] + sDlx + 
643           xshift[geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 2] - 1];
644       } else if (geant3->Gctmed()->numed == idtmed[1003-1]) {
645         nmodsx = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
646         xsx = sVloc[0] + sDlx;
647       } else {
648         nmodsx = 0;
649       }
650       
651       //        Momentum and direction of incidence 
652       psx = geant3->Gctrak()->vect[6];
653       // Computing 2nd power 
654       r1 = dir[0];
655       // Computing 2nd power 
656       r2 = dir[2];
657       thetasx = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
658       phisx = TMath::ATan2(dir[2], dir[0]);
659       if (phisx < 0.) phisx += 2*TMath::Pi();
660       
661       //        Particle type and parent 
662       //idpartsx = geant3->Gckine()->ipart;
663       r1 = geant3->Gctrak()->upwght / 100.;
664       idpartgx = Int_t(r1+0.5);
665       if (!geant3->Gctrak()->upwght) {
666         sPparent = geant3->Gctrak()->vect[6];
667         sThincParent = thetasx;
668       }
669       // new 
670       for (ll = 0; ll <nrooth; ++ll) rrhh[ll] = 0;
671       irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
672       irivol[1] = nmodsx;
673       ihitrak = gAlice->CurrentTrack();
674       rrhh[0] = 1.;
675       // Flag to say that this is MIP 
676       rrhh[1] = (Float_t) geant3->Gckine()->ipart;
677       rrhh[2] = xsx;
678       rrhh[3] = ysx;
679       rrhh[4] = (Float_t) nmodsx;
680       // Module Number 
681       rrhh[5] = thetasx;
682       rrhh[6] = geant3->Gctrak()->tofg;
683       // in seconds 
684       rrhh[7] = (Float_t) idpartgx;
685       //     mips specific 
686       rrhh[8] = phisx;
687       rrhh[9] = psx;
688       //     charge of current particle in electron charge unit; 
689       rrhh[10] = geant3->Gckine()->charge;
690       rrhh[11] = sVloc[1];
691       // Zo of ckov generation 
692       rrhh[12] = 1.;
693       // ckov generation 
694       AddHit(ihitrak, irivol, rrhh);
695       // enew 
696     }
697   }
698   
699   //        Current particle is cherenkov photon 
700   if (geant3->Gckine()->ipart == 50) {
701     gMC->Gmtod(geant3->Gctrak()->vect, sVloc, 1);
702     //         WRITE(6,* ) UPWGHT, VLOC(2), NUMED, DESTEP 
703     //        Photon crosses ch4-csi boundary 
704     //           take into account fresnel losses with complex refraction index
705     if (geant3->Gctrak()->inwvol == 1 && geant3->Gctmed()->numed == idtmed[1006-1]) {
706       
707       //     fresnel losses commented out for the moment 
708       //    make sure that qe correction for fresnel losses is also switched off !
709       //            CALL FRESNELCSI 
710       //            IF (ISTOP .EQ. 2) RETURN 
711       //        Put transmission of electrodes in by hand 
712       gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
713       cophi = TMath::Cos(TMath::ATan2(dir[0], dir[1]));
714       t = (1. - .025 / cophi) * (1. - .05 /  cophi);
715       gMC->Rndm(ranf, 1);
716       if (ranf[0] > t) {
717         if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5)<MAXPH) 
718           sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 15;
719         geant3->Gctrak()->istop = 2;
720         return;
721       }
722     }
723     
724     //        Photon lost energy in CsI 
725     if (geant3->Gctrak()->destep > 0. && geant3->Gctmed()->numed == idtmed[1006-1]) {
726       geant3->Gctrak()->istop = 2;
727       r1 = geant3->Gctrak()->upwght / 100.;
728       if (Int_t(r1+0.5) > 50) {
729         ++sNfeedd;
730       } else {
731         ++sNdir;
732       }
733       //            WRITE(6,*) 'PHOTON',UPWGHT, MAXPH 
734       if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH) 
735         sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] = 0;
736       //            write(6,*) 'photon detected' 
737       for(j=0;j<3;j++) sVectIn[j]=geant3->Gctrak()->vect[j];
738       // new 
739       // copied from miphit in Freon or Quartz 
740       //        Where did it hit ? 
741       gMC->Gmtod(&geant3->Gctrak()->vect[3], dir, 2);
742       
743       //        Momentum and direction of incidence 
744       for (ll = 0; ll < nrooth; ++ll) rrhh[ll]=0;
745       irivol[0] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
746       // ??? 
747       irivol[1] = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
748       // ??? 
749       ihitrak = gAlice->CurrentTrack();
750       rrhh[0] = 2.;
751       // Flag to say that this is CK 
752       rrhh[1] = (Float_t) geant3->Gckine()->ipart;
753       rrhh[2] = sVloc[0] + sDlx;
754       rrhh[3] = sVloc[2] + sDly;
755       rrhh[4] = (Float_t) geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
756       // ??? Module Number 
757       // Computing 2nd power 
758       r1 = dir[0];
759       // Computing 2nd power 
760       r2 = dir[2];
761       rrhh[5] = TMath::ATan2(TMath::Sqrt(r1 * r1 + r2 * r2), dir[1]);
762       // THETASX 
763       rrhh[6] = geant3->Gctrak()->tofg;
764       // in seconds 
765       r1 = geant3->Gctrak()->upwght / 100.;
766       rrhh[7] = (Float_t) Int_t(r1+0.5);
767       //     ckov specific 
768       // Feedback ??? 
769       rrhh[8] = geant3->Gctrak()->getot;
770       rrhh[9] = 0.;
771       // Stop in CsI 
772       AddHit(ihitrak, irivol, rrhh);
773       //     end of new 
774       RichIntegration();
775       return;
776     }
777     if (geant3->Gctrak()->upwght && Int_t(geant3->Gctrak()->upwght+0.5) < MAXPH) 
778       {
779         //        Losses by reflection and absorption and leaving detector
780         if (sIloss[Int_t(geant3->Gctrak()->upwght+0.5) - 1] == 22) {
781           i1 = geant3->Gctrak()->nmec;
782           for (i = 0; i < i1; ++i) {
783             //        Reflection loss 
784             if (geant3->Gctrak()->lmec[i] == 106) {
785               sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=10;
786               if (geant3->Gctmed()->numed == idtmed[1004-1]) 
787                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=1;
788               
789               if (geant3->Gctmed()->numed == idtmed[1003-1]) 
790                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=2;
791               
792               //        Absorption loss 
793             } else if (geant3->Gctrak()->lmec[i] == 101) {
794               geant3->Gctrak()->istop = 2;
795               sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=20;
796               if (geant3->Gctmed()->numed == idtmed[1004-1]) 
797                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=11;
798               
799               if (geant3->Gctmed()->numed == idtmed[1003-1]) 
800                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=12;
801               
802               if (geant3->Gctmed()->numed == idtmed[1005-1]) 
803                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
804               
805               if (geant3->Gctmed()->numed == idtmed[1009-1]) 
806                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=13;
807               
808               if (geant3->Gctmed()->numed == idtmed[1001-1]) 
809                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=14;
810               
811               //        CsI inefficiency 
812               if (geant3->Gctmed()->numed == idtmed[1006-1]) 
813                 sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=16;
814               
815               
816               //        Photon goes out of tracking scope 
817             } else if (geant3->Gctrak()->lmec[i] == 30) 
818               sIloss[Int_t(geant3->Gctrak()->upwght+0.5)-1]=21;
819             
820           }
821         }
822       }
823   }
824 }
825
826
827 //_____________________________________________________________________________
828 void AliRICH::RichIntegration()
829 {
830   //
831   // Integrates over RICH pads
832   //
833
834   TGeant3 *geant3 = (TGeant3*) gMC;
835   
836   Int_t i1, i2;
837   Float_t r1;
838   
839   // Local variables 
840   Float_t rrhh[25];
841   Float_t qtot;
842   Int_t ifeed;
843   Float_t x0;
844   Int_t ixmin, ixmax, iymin, iymax;
845   Int_t ll, ix, iy;
846   Float_t qp;
847   Float_t source[3];
848   Int_t irivol[2];
849   Float_t y0a, qtot_check, xi1, xi2, yi1, yi2;
850   Int_t nph = 0, iqp;
851   Int_t ihitrak, modulen;
852   //Int_t isec[MAXSEC];
853   //     ILOSS = 0    NOT LOST 
854   //             1    REFLECTED FREON-QUARZ 
855   //             2    REFLECTED QUARZ-METHANE 
856   //             3    REFLECTED METHANE-CSI 
857   //            11    ABSORBED IN FREON 
858   //            12    ABSORBED IN QUARZ 
859   //            13    ABSORBED IN METHANE 
860   //            21    CSI QE 
861   //     IPROD = 1    PRODUCED IN FREON 
862   //     IPROD = 2    PRODUCED IN QUARZ 
863   //       get charge for MIP or cherenkov photon 
864   
865   if (geant3->Gckine()->ipart == 50) {
866     GetCharge(qtot);
867     modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 3];
868     //sQint[1] = qtot;
869   } else {
870     GetChargeMip(qtot);
871     modulen = geant3->Gcvolu()->number[geant3->Gcvolu()->nlevel - 4];
872     //sQint[0] = qtot;
873   }
874   //
875   gMC->Gmtod(sVectIn, sVloc, 1);
876   if (TMath::Abs(sVloc[0]) >= sDlx) return;
877   
878   if (TMath::Abs(sVloc[2]) >= sDly) return;
879   
880   sVloc[0] += sDlx;
881   x0 = sVloc[0];
882   sVloc[2] += sDly;
883   //y0 = sVloc[2];
884   AnodicWires(y0a);
885   ixmin = (Int_t) ((x0 - fSxcharge * 5.) / sDxp) + 1;
886   if (ixmin < 1) ixmin = 1;
887   ixmax = (Int_t) ((x0 + fSxcharge * 5.) / sDxp) + 1;
888   if (ixmax > sNpx) ixmax = sNpx;
889   iymin = (Int_t) ((y0a - sSycharge * 5.) / sDyp) + 1;
890   if (iymin < 1) iymin = 1;
891   iymax = (Int_t) ((y0a + sSycharge * 5.) / sDyp) + 1;
892   if (iymax > sNpy) iymax = sNpy;
893   qtot_check = 0.;
894   i1 = ixmax;
895   for (ix = ixmin; ix <= i1; ++ix) {
896     i2 = iymax;
897     for (iy = iymin; iy <= i2; ++iy) {
898       xi1 = (sXpad[ix - 1] - x0) / sDpad;
899       xi2 = xi1 + sDxp / sDpad;
900       yi1 = (sYpad[iy - 1] - y0a) / sDpad;
901       yi2 = yi1 + sDyp / sDpad;
902       qp = qtot * FMathieson(xi1, xi2) * FMathieson(yi1, yi2);
903       iqp = (Int_t) TMath::Abs(qp);
904       qtot_check += TMath::Abs(qp);
905       
906       //     FILL COMMON FOR PADS 
907       
908       if (iqp) {
909         if (iqp > 4095) {
910           iqp = 4096;
911         }
912         ++sNpads;
913         if (sNpads > MAXSEC) {
914           //                  write(6,*) 'attention npads',npads 
915           sNpads = MAXSEC;
916         }
917         sIpx[sNpads - 1] = ix;
918         sIpy[sNpads - 1] = iy;
919         sIqpad[sNpads - 1] = iqp;
920         //     TAG THE ORIGIN OF THE CHARGE DEPOSITION 
921         r1 = geant3->Gctrak()->upwght / 100.;
922         ifeed = Int_t(r1+0.5);
923         sNphlink[sNpads - 1] = 0;
924         if (geant3->Gckine()->ipart != 50) {
925           //     MIP 
926           //isec[sNpads - 1] = sNsecon;
927         } else {
928           if (ifeed < 50) {
929             //     CERENKOV PHOTON 
930             
931             nph = Int_t(geant3->Gctrak()->upwght+0.5);
932             sNphlink[sNpads - 1] = nph;
933             sXphit[nph - 1] = sVloc[0];
934             sYphit[nph - 1] = sVloc[2];
935             //isec[sNpads - 1] = sNsecon + 100;
936           } else if (ifeed == 51) {
937             //     FEEDBACK FROM CERENKOV 
938             
939             //isec[sNpads - 1] = sNsecon + 300;
940           } else if (ifeed == 52) {
941             //     FEEDBACK FROM MIP 
942             
943             //isec[sNpads - 1] = sNsecon + 200;
944           }
945         }
946         // Generate the hit for Root IO 
947         for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
948         irivol[0] = modulen;
949         irivol[1] = nmodsx;
950         rrhh[0] = 0.;
951         // Flag to say this is a  hit 
952         rrhh[1] = xsx;
953         rrhh[2] = ysx;
954         rrhh[3] = psx;
955         rrhh[4] = thetasx;
956         rrhh[5] = phisx;
957         rrhh[6] = (Float_t) idpartgx;
958         rrhh[7] = sXsox;
959         rrhh[8] = sYsox;
960         rrhh[9] = sZsox;
961         rrhh[10] = (Float_t) sIpx[sNpads - 1];
962         rrhh[11] = (Float_t) sIpy[sNpads - 1];
963         rrhh[12] = (Float_t) sIqpad[sNpads - 1];
964         ihitrak = gAlice->CurrentTrack();
965         if (sNphlink[sNpads - 1] > 0) {
966           rrhh[13] = sPparent;
967           rrhh[14] = sThincParent;
968           rrhh[15] = (Float_t) sIloss[nph - 1];
969           rrhh[16] = (Float_t) sIprod[nph - 1];
970           rrhh[17] = sXphit[nph - 1];
971           rrhh[18] = sYphit[nph - 1];
972           ihitrak = sMckov[nph - 1];
973         }
974         //     This is the current position of the hit 
975         rrhh[19] = geant3->Gctrak()->vect[0];
976         rrhh[20] = geant3->Gctrak()->vect[1];
977         rrhh[21] = geant3->Gctrak()->vect[2];
978         //               PRINT *,' ' 
979         //               PRINT *,'RXAHIT(oldhit)=[',RRHH(20),',',RRHH(21),',', 
980         //     +              RRHH(22),']' 
981         //          PRINT *, 'ITRA = ',ITRA,'ISTAK = ',ISTAK 
982         AddHit(ihitrak, irivol, rrhh);
983         // new Padhits 
984         for (ll = 0; ll < 25; ++ll) rrhh[ll] = 0;
985         rrhh[0] = 3.;
986         rrhh[1] = (Float_t) geant3->Gckine()->ipart;
987         rrhh[2] = (Float_t) sIpx[sNpads - 1];
988         rrhh[3] = (Float_t) sIpy[sNpads - 1];
989         rrhh[4] = (Float_t) modulen;
990         rrhh[5] = -1.;
991         // !!!Prod 
992         rrhh[6] = (Float_t) sIqpad[sNpads - 1];
993         AddHit(ihitrak, irivol, rrhh);
994         // enew 
995       }
996     }
997   }
998   //      IF (IPART .NE. 50) WRITE(6,*) 'CHECK',QTOT,QTOT_CHECK 
999   
1000   //     GENERATE FEEDBACK PHOTONS 
1001   
1002   source[0] = x0 - sDlx;
1003   source[1] = sVloc[1] - .2;
1004   source[2] = y0a - sDly;
1005   gMC->Gdtom(source, source, 1);
1006   FeedBack(source, qtot);
1007   return;
1008 }
1009
1010 //_____________________________________________________________________________
1011 void AliRICH::AnodicWires(Float_t &y0a)
1012 {
1013   //
1014   // Position of anodic wires
1015   //
1016   Int_t i1;
1017   
1018   // Local variables 
1019   Int_t i;
1020   Float_t ass_i, y0, ass_i1;
1021   
1022   y0 = sVloc[2];
1023   y0a = -1.;
1024   i1 = (sNpy << 1) + 1;
1025   for (i = 1; i <= i1; ++i) {
1026     if (y0 > sYanode[i - 1] && y0 <= sYanode[i]) {
1027       ass_i1 = TMath::Abs(sYanode[i] - y0);
1028       ass_i = TMath::Abs(sYanode[i - 1] - y0);
1029       if (ass_i1 <= ass_i) y0a = sYanode[i];
1030       if (ass_i1 > ass_i) y0a = sYanode[i - 1];
1031       return;
1032     }
1033   }
1034
1035
1036 //_____________________________________________________________________________
1037 void AliRICH::GetChargeMip(Float_t &qtot)
1038 {
1039   //
1040   // Calculates the charge deposited by a MIP
1041   //
1042
1043
1044   Int_t i1;
1045   
1046   // Local variables 
1047   Float_t ranf[2];
1048   Int_t i;
1049   Int_t nel;
1050   
1051   //     NUMBER OF ELECTRONS 
1052   
1053   nel = (Int_t) (sDetot * 1e9 / 26.);
1054   qtot = 0.;
1055   if (!nel) nel = 1;
1056   i1 = nel;
1057   for (i = 1; i <= i1; ++i) {
1058     gMC->Rndm(ranf, 1);
1059     qtot -= fChslope * TMath::Log(ranf[0]);
1060   }
1061
1062
1063 //_____________________________________________________________________________
1064 void AliRICH::GetCharge(Float_t &qtot)
1065 {
1066   //
1067   // Charge deposited
1068   //
1069
1070
1071   Float_t ranf[1];
1072   
1073   gMC->Rndm(ranf, 1);
1074   qtot = -fChslope * TMath::Log(ranf[0]);
1075
1076
1077 //_____________________________________________________________________________
1078 void AliRICH::FeedBack(Float_t *source, Float_t qtot)
1079 {
1080   //
1081   // Generate FeedBack photons
1082   //
1083
1084   TGeant3 *geant3 = (TGeant3*) gMC;
1085
1086   Int_t i1, j;
1087   Float_t r1, r2;
1088   
1089   // Local variables 
1090   Float_t cthf, ranf[2], phif, enfp = 0, sthf;
1091   Int_t i, ifeed;
1092   Float_t e1[3], e2[3], e3[3];
1093   Float_t vmod, uswop;
1094   Float_t fp, random;
1095   Float_t dir[3], phi;
1096   Int_t nfp;
1097   Float_t pol[3], supwght;
1098   
1099   //     DETERMINE NUMBER OF FEEDBACK PHOTONS 
1100   
1101   // Function Body 
1102   fp = fAlphaFeed * qtot;
1103   nfp = gRandom->Poisson(fp);
1104   
1105   //     GENERATE PHOTONS 
1106   
1107   geant3->Gckin2()->ngphot = 0;
1108   i1 = nfp;
1109   for (i = 1; i <= i1; ++i) {
1110     
1111     //     DIRECTION 
1112     gMC->Rndm(ranf, 2);
1113     cthf = ranf[0] * 2. - 1.;
1114     if (cthf < 0.)  continue;
1115     sthf = TMath::Sqrt((1. - cthf) * (cthf + 1.));
1116     phif = ranf[1] * 2 * TMath::Pi();
1117     ++geant3->Gckin2()->ngphot;
1118     if (geant3->Gckin2()->ngphot > 800) {
1119       printf("ATTENTION NGPHOT ! %d %f %d\n",
1120              geant3->Gckin2()->ngphot,fp,nfp);
1121       break;
1122     }
1123     
1124     //     POSITION 
1125     for(j=0;j<3;j++) 
1126       geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][j]=source[j];
1127     
1128     //     ENERGY 
1129     gMC->Rndm(&random, 1);
1130     if (random <= .57) {
1131       enfp = 7.5e-9;
1132     } else if (random > .57 && random <= .7) {
1133       enfp = 6.4e-9;
1134     } else if (random > .7) {
1135       enfp = 7.9e-9;
1136     }
1137     geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][6] = enfp;
1138     dir[0] = sthf * TMath::Sin(phif);
1139     dir[1] = cthf;
1140     dir[2] = sthf * TMath::Cos(phif);
1141     gMC->Gdtom(dir, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][3], 2);
1142     
1143     //     POLARISATION 
1144     e1[0] = 0.;
1145     e1[1] = -dir[2];
1146     e1[2] = dir[1];
1147     
1148     e2[0] = -dir[1];
1149     e2[1] = dir[0];
1150     e2[2] = 0.;
1151     
1152     e3[0] = dir[1];
1153     e3[1] = 0.;
1154     e3[2] = -dir[0];
1155     
1156     vmod=0;
1157     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
1158     if (!vmod) for(j=0;j<3;j++) {
1159       uswop=e1[j];
1160       e1[j]=e3[j];
1161       e3[j]=uswop;
1162     }
1163     vmod=0;
1164     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
1165     if (!vmod) for(j=0;j<3;j++) {
1166       uswop=e2[j];
1167       e2[j]=e3[j];
1168       e3[j]=uswop;
1169     }
1170     
1171     vmod=0;
1172     for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
1173     vmod=TMath::Sqrt(1/vmod);
1174     for(j=0;j<3;j++) e1[j]*=vmod;
1175     
1176     vmod=0;
1177     for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
1178     vmod=TMath::Sqrt(1/vmod);
1179     for(j=0;j<3;j++) e2[j]*=vmod;
1180     
1181     gMC->Rndm(ranf, 1);
1182     phi = ranf[0] * 2 * TMath::Pi();
1183     r1 = TMath::Sin(phi);
1184     r2 = TMath::Cos(phi);
1185     for(j=0;j<3;j++) pol[j]=e1[j]*r1+e2[j]*r2;
1186     gMC->Gdtom(pol, &geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][7], 2);
1187     
1188     //     TIME OF FLIGHT 
1189     geant3->Gckin2()->xphot[geant3->Gckin2()->ngphot-1][10] = 0.;
1190     
1191     //     PUT PHOTON ON THE STACK AND LABEL IT AS FEEDBACK (51,52) 
1192     supwght = geant3->Gctrak()->upwght;
1193     r1 = geant3->Gctrak()->upwght / 100.;
1194     ifeed = Int_t(r1+0.5);
1195     ++sNfeed;
1196     if (geant3->Gckine()->ipart == 50 && ifeed != 52) {
1197       geant3->Gctrak()->upwght = 5100.;
1198     } else {
1199       geant3->Gctrak()->upwght = 5200.;
1200     }
1201     //    geant3->Gskpho(geant3->Gckin2()->ngphot);
1202     geant3->Gctrak()->upwght = supwght;
1203     
1204   }
1205   geant3->Gckin2()->ngphot = 0;
1206 }
1207
1208 //_____________________________________________________________________________
1209 Float_t AliRICH::FMathieson(Float_t lambda1, Float_t lambda2)
1210 {
1211   //
1212   // Mathieson charge distribution
1213   //
1214   // CALCULATES INTEGRAL OF GATTI/MATHIESON CHARGE DISTRIBUTION 
1215   // FOR CPC AND CSC CHAMBERS 
1216   // INTEGRATION LIMITS LAMBDA1,LAMBDA2=  X/D WHERE D IS THE ANODE CATHODE
1217   //                                          SEPARATION 
1218   // K3: GEOMETRY PARAMETER 
1219   //
1220   Float_t u1, u2;
1221   //    const Float_t k1=0.2828;
1222   const Float_t k2=0.962;
1223   const Float_t k3=0.6;
1224   const Float_t k4=0.379;
1225   const Float_t sqrtk3=TMath::Sqrt(k3);
1226   
1227   
1228   u1 = tanh(lambda1 * k2) * sqrtk3;
1229   u2 = tanh(lambda2 * k2) * sqrtk3;
1230   
1231   return (atan(u2) - atan(u1)) * 2 * k4;
1232   
1233 }
1234
1235
1236 //_____________________________________________________________________________
1237 void AliRICH::UpdateMipHit(Float_t *hits)
1238 {
1239   //
1240   // Update some parameters when the current mip hits the detector.
1241   //
1242   TClonesArray &lhits = *fMips;
1243   AliRICHmip *lmip = (AliRICHmip *) lhits.AddrAt(fNmips-1);
1244   lmip->SetX     (      hits[1]);
1245   lmip->SetY     (      hits[2]);
1246   lmip->SetModule((int) hits[3]);
1247   lmip->SetTheta (      hits[4]);
1248   lmip->SetPhi   (      hits[5]);
1249 }
1250
1251 //*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
1252
1253 //_____________________________________________________________________________
1254 void AliRICH::MakeBranch(Option_t *option)
1255 {
1256   //
1257   // Create a new branch in the current Root Tree.
1258   // The branch of fHits is automatically split
1259   //
1260   AliDetector::MakeBranch(option);
1261   
1262   Int_t buffersize = 4000;
1263   if (gAlice->TreeH()) {
1264     if(fMips){
1265       gAlice->TreeH()->Branch("RICHmip",&fMips, buffersize);
1266       printf("Making Branch %s for mips\n","RICHmip");
1267     }
1268     if(fCkovs){
1269       gAlice->TreeH()->Branch("RICHckov",&fCkovs, buffersize);
1270       printf("Making Branch %s for ckovs\n","RICHckov");
1271     }
1272     if(fPadhits){
1273       gAlice->TreeH()->Branch("RICHpadhit",&fPadhits, buffersize);
1274       printf("Making Branch %s for padhits\n","RICHpadhit");
1275     }
1276   }     
1277 }
1278
1279 //_____________________________________________________________________________
1280 void AliRICH::SetTreeAddress()
1281 {
1282   //
1283   // Set branch address for the Hits and Digits Tree.
1284   //
1285   AliDetector::SetTreeAddress();
1286   TBranch *branch;
1287   TTree *treeH = gAlice->TreeH();
1288   if (treeH) {
1289     if (fMips) {
1290       branch = treeH->GetBranch("RICHmip");
1291       if (branch) branch->SetAddress(&fMips);
1292     }
1293     if (fCkovs) {
1294       branch = treeH->GetBranch("RICHckov");
1295       if (branch) branch->SetAddress(&fCkovs);
1296     }
1297     if (fPadhits) {
1298       branch = treeH->GetBranch("RICHpadhit");
1299       if (branch) branch->SetAddress(&fPadhits);
1300     }
1301   }
1302 }
1303
1304 //_____________________________________________________________________________
1305 void AliRICH::ResetHits()
1306 {
1307   //
1308   // Reset number of mips and the mips array for RICH
1309   //
1310   AliDetector::ResetHits();
1311   fNmips    = 0;
1312   if (fMips)    fMips->Clear();
1313   // Reset number of Ckovs and the Ckovs array for RICH
1314   fNckovs   = 0;
1315   if (fCkovs)   fCkovs->Clear();
1316   // Reset number of Padhits and the Padhits array for RICH
1317   fNpadhits = 0;
1318   if (fPadhits) fPadhits->Clear();
1319 }
1320
1321 //_____________________________________________________________________________
1322 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1323 {
1324   //
1325   // Distance from mouse RICH on the display
1326   // Dummy routine
1327   //
1328   return 9999;
1329 }
1330  
1331 //_____________________________________________________________________________
1332 void AliRICH::PreTrack()
1333 {
1334   //
1335   // Called before transporting a track
1336   //
1337   sDetot=0;
1338   sNsecon=0;
1339   sNpads=0;
1340   sNphoton=0;
1341   sNfeed=0;
1342   sNfeedd=0;
1343   sNdir=0;
1344 }
1345
1346 //_____________________________________________________________________________
1347 void AliRICH::Streamer(TBuffer &R__b)
1348 {
1349   //
1350   // Stream an object of class AliRICH.
1351   //
1352   if (R__b.IsReading()) {
1353     Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1354     AliDetector::Streamer(R__b);
1355     R__b >> fNmips;
1356     R__b >> fNckovs;
1357     R__b >> fNpadhits;
1358     R__b >> fMips; //diff
1359     R__b >> fCkovs; //diff
1360     R__b >> fPadhits; //diff
1361     R__b >> fChslope;
1362     R__b >> fAlphaFeed;
1363     R__b >> fSxcharge;
1364     R__b >> fIritri;
1365   } else {
1366     R__b.WriteVersion(AliRICH::IsA());
1367     AliDetector::Streamer(R__b);
1368     R__b << fNmips;
1369     R__b << fNckovs;
1370     R__b << fNpadhits;
1371     R__b << fMips; //diff
1372     R__b << fCkovs; //diff
1373     R__b << fPadhits; //diff
1374     R__b << fChslope;
1375     R__b << fAlphaFeed;
1376     R__b << fSxcharge;
1377     R__b << fIritri;
1378   }
1379 }
1380
1381 ClassImp(AliRICHv1)
1382
1383 ///////////////////////////////////////////////////////////////////////////////
1384 //                                                                           //
1385 //  Ring Imaging Cherenkov                                                   //
1386 //  This class contains version 1 of the Ring Imaging Cherenkov              //
1387 //                                                                           //
1388 //Begin_Html
1389 /*
1390 <img src="picts/AliRICHv1Class.gif">
1391 */
1392 //End_Html
1393 //                                                                           //
1394 ///////////////////////////////////////////////////////////////////////////////
1395
1396 //_____________________________________________________________________________
1397 AliRICHv1::AliRICHv1() : AliRICH()
1398 {
1399   //
1400   // Default Constructor RICH for version 1
1401   //
1402 }
1403  
1404 //_____________________________________________________________________________
1405 AliRICHv1::AliRICHv1(const char *name, const char *title)
1406        : AliRICH(name,title)
1407 {
1408   //
1409   // Standard constructor for RICH version 1
1410   //
1411 }
1412
1413 //_____________________________________________________________________________
1414 AliRICHv1::~AliRICHv1()
1415 {
1416   //
1417   // Standard destructor for RICH version 1
1418   //
1419 }
1420  
1421 //_____________________________________________________________________________
1422 void AliRICHv1::CreateGeometry()
1423 {
1424   //
1425   // Create the geometry for RICH version 1
1426   //
1427   // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
1428   //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
1429   //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
1430   //
1431   //Begin_Html
1432   /*
1433     <img src="picts/AliRICHv1.gif">
1434   */
1435   //End_Html
1436   //Begin_Html
1437   /*
1438     <img src="picts/AliRICHv1Tree.gif">
1439   */
1440   //End_Html
1441
1442
1443   Int_t *idtmed = fIdtmed->GetArray()-999;
1444   
1445   Int_t i;
1446   Float_t zs;
1447   Int_t idrotm[1099];
1448   Float_t par[3];
1449   
1450   // --- Define the RICH detector 
1451   //     External aluminium box 
1452   par[0] = 71.1;
1453   par[1] = 11.5;
1454   par[2] = 73.15;
1455   gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
1456   
1457   //     Sensitive part of the whole RICH 
1458   par[0] = 64.8;
1459   par[1] = 11.5;
1460   par[2] = 66.55;
1461   gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
1462   
1463   //     Honeycomb 
1464   par[0] = 63.1;
1465   par[1] = .188;
1466   par[2] = 66.55;
1467   gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
1468   
1469   //     Aluminium sheet 
1470   par[0] = 63.1;
1471   par[1] = .025;
1472   par[2] = 66.55;
1473   gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
1474   
1475   //     Quartz 
1476   par[0] = 63.1;
1477   par[1] = .25;
1478   par[2] = 65.5;
1479   gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
1480   
1481   //     Spacers (cylinders) 
1482   par[0] = 0.;
1483   par[1] = .5;
1484   par[2] = .5;
1485   gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
1486   
1487   //     Opaque quartz 
1488   par[0] = 61.95;
1489   par[1] = .2;
1490   par[2] = 66.5;
1491   gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
1492   
1493   //     Frame of opaque quartz 
1494   par[0] = 20.65;
1495   par[1] = .5;
1496   par[2] = 66.5;
1497   gMC->Gsvolu("OQUF", "BOX ", idtmed[1007], par, 3);
1498   
1499   //     Little bar of opaque quartz 
1500   par[0] = 63.1;
1501   par[1] = .25;
1502   par[2] = .275;
1503   gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
1504   
1505   //     Freon 
1506   par[0] = 20.15;
1507   par[1] = .5;
1508   par[2] = 65.5;
1509   gMC->Gsvolu("FREO", "BOX ", idtmed[1003], par, 3);
1510   
1511   //     Methane 
1512   par[0] = 64.8;
1513   par[1] = 5.;
1514   par[2] = 64.8;
1515   gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
1516   
1517   //     Methane gap 
1518   par[0] = 64.8;
1519   par[1] = .2;
1520   par[2] = 64.8;
1521   gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1522   
1523   //     CsI photocathode 
1524   par[0] = 64.8;
1525   par[1] = .25;
1526   par[2] = 64.8;
1527   gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1528   
1529   //     Anode grid 
1530   par[0] = 0.;
1531   par[1] = .0025;
1532   par[2] = 20.;
1533   gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1534   
1535   // --- Places the detectors defined with GSVOLU 
1536   //     Place material inside RICH 
1537   gMC->Gspos("SRIC", 1, "RICH", 0., 0., 0., 0, "ONLY");
1538   
1539   gMC->Gspos("ALUM", 1, "SRIC", 0., -6.075, 0., 0, "ONLY");
1540   gMC->Gspos("HONE", 1, "SRIC", 0., -5.862, 0., 0, "ONLY");
1541   gMC->Gspos("ALUM", 2, "SRIC", 0., -5.649, 0., 0, "ONLY");
1542   gMC->Gspos("OQUA", 1, "SRIC", 0., -5.424, 0., 0, "ONLY");
1543   
1544   AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1545   
1546   for (i = 1; i <= 9; ++i) {
1547     zs = (5 - i) * 14.4;
1548     gMC->Gspos("SPAC", i, "FREO", 6.7, 0., zs, idrotm[1019], "ONLY");
1549   }
1550   for (i = 10; i <= 18; ++i) {
1551     zs = (14 - i) * 14.4;
1552     gMC->Gspos("SPAC", i, "FREO", -6.7, 0., zs, idrotm[1019], "ONLY");
1553   }
1554   
1555   gMC->Gspos("FREO", 1, "OQUF", 0., 0., 0., 0, "ONLY");
1556   gMC->Gspos("OQUF", 1, "SRIC", 41.3, -4.724, 0., 0, "ONLY");
1557   gMC->Gspos("OQUF", 2, "SRIC", 0., -4.724, 0., 0, "ONLY");
1558   gMC->Gspos("OQUF", 3, "SRIC", -41.3, -4.724, 0., 0, "ONLY");
1559   gMC->Gspos("BARR", 1, "QUAR", 0., 0., -21.65, 0, "ONLY");
1560   gMC->Gspos("BARR", 2, "QUAR", 0., 0., 21.65, 0, "ONLY");
1561   gMC->Gspos("QUAR", 1, "SRIC", 0., -3.974, 0., 0, "ONLY");
1562   gMC->Gspos("GAP ", 1, "META", 0., 4.8, 0., 0, "ONLY");
1563   gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1564   gMC->Gspos("CSI ", 1, "SRIC", 0., 6.526, 0., 0, "ONLY");
1565   
1566   //     Place RICH inside ALICE apparatus 
1567   
1568   AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1569   AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1570   AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1571   AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1572   AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1573   AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1574   AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1575   
1576   gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
1577   gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
1578   gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
1579   gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
1580   gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
1581   gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
1582   gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
1583   
1584 }
1585
1586 //_____________________________________________________________________________
1587 void AliRICHv1::DrawModule()
1588 {
1589   //
1590   // Draw a shaded view of the Ring Imaging Cherenkov
1591   //
1592
1593   TGeant3 *geant3 = (TGeant3*) gMC;
1594
1595   // Set everything unseen
1596   gMC->Gsatt("*", "seen", -1);
1597   // 
1598   // Set ALIC mother transparent
1599   gMC->Gsatt("ALIC","SEEN",0);
1600   //
1601   // Set the volumes visible
1602
1603   gMC->Gsatt("RICH","seen",0);
1604   gMC->Gsatt("SRIC","seen",0);
1605   gMC->Gsatt("HONE","seen",1);
1606   gMC->Gsatt("ALUM","seen",1);
1607   gMC->Gsatt("QUAR","seen",1);
1608   gMC->Gsatt("SPAC","seen",1);
1609   gMC->Gsatt("OQUA","seen",1);
1610   gMC->Gsatt("OQUF","seen",1);
1611   gMC->Gsatt("BARR","seen",1);
1612   gMC->Gsatt("FREO","seen",1);
1613   gMC->Gsatt("META","seen",1);
1614   gMC->Gsatt("GAP ","seen",1);
1615   gMC->Gsatt("CSI ","seen",1);
1616   gMC->Gsatt("GRID","seen",1);
1617   
1618   //
1619   geant3->Gdopt("hide", "on");
1620   geant3->Gdopt("shad", "on");
1621   geant3->Gsatt("*", "fill", 7);
1622   geant3->SetClipBox(".");
1623   geant3->Gdopt("hide", "on");
1624   geant3->Gdopt("shad", "on");
1625   geant3->Gsatt("*", "fill", 7);
1626   geant3->SetClipBox(".");
1627   //  geant3->SetClipBox("*", 0, 2000, -2000, 2000, -2000, 2000);
1628   geant3->DefaultRange();
1629   gMC->Gdraw("alic", 60, 50, 0, 10, 0, .03, .03);
1630   geant3->Gdhead(1111, "Ring Imaging Cherenkov version 1");
1631   geant3->Gdman(16, 6, "MAN");
1632   geant3->Gdopt("hide", "off");
1633 }
1634
1635 //_____________________________________________________________________________
1636 void AliRICHv1::CreateMaterials()
1637 {
1638   //
1639   // *** DEFINITION OF AVAILABLE RICH MATERIALS *** 
1640   // ORIGIN    : NICK VAN EIJNDHOVEN 
1641   // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it) 
1642   //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it) 
1643   //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it) 
1644   //
1645   Int_t   ISXFLD = gAlice->Field()->Integ();
1646   Float_t SXMGMX = gAlice->Field()->Max();
1647   
1648   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,
1649                          6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1650   Float_t rindex_quarz[14] = { 1.528309,1.533333,
1651                                1.538243,1.544223,1.550568,1.55777,
1652                                1.565463,1.574765,1.584831,1.597027,
1653                                1.611858,1.6277,1.6472,1.6724 };
1654   Float_t rindex_quarzo[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1655   Float_t rindex_methane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1656   Float_t rindex_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1657   Float_t absco_freon[14] = { 179.0987,179.0987,
1658                               179.0987,179.0987,179.0987,35.7,12.54,5.92,4.92,3.86,1.42,.336,.134,0. };
1659   Float_t absco_quarz[14] = { 20.126,16.27,13.49,11.728,9.224,8.38,7.44,7.17,
1660                               6.324,4.483,1.6,.323,.073,0. };
1661   Float_t absco_quarzo[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1662                                1e-5,1e-5,1e-5,1e-5,1e-5 };
1663   Float_t absco_csi[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1664                             1e-4,1e-4,1e-4,1e-4 };
1665   Float_t absco_methane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1666                                 1e6,1e6,1e6 };
1667   Float_t absco_gri[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1668                             1e-4,1e-4,1e-4,1e-4 };
1669   Float_t effic_all[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1670   Float_t effic_csi[14] = { 4.74e-4,.00438,.009,.0182,.0282,.0653,.1141,.163,
1671                             .2101,.2554,.293,.376,.3861,.418 };
1672   Float_t effic_gri[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1673   
1674   Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon, 
1675     zmet[2], zqua[2];
1676   Int_t nlmatfre;
1677   Float_t densquao;
1678   Int_t nlmatmet, nlmatqua;
1679   Float_t wmatquao[2], rindex_freon[14];
1680   Int_t i;
1681   Float_t aquao[2], epsil, stmin, zquao[2];
1682   Int_t nlmatquao;
1683   Float_t radlal, densal, tmaxfd, deemax, stemax;
1684   Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1685   
1686   Int_t *idtmed = fIdtmed->GetArray()-999;
1687   
1688   TGeant3 *geant3 = (TGeant3*) gMC;
1689
1690   // --- Photon energy (GeV) 
1691   // --- Refraction indexes 
1692   for (i = 0; i < 14; ++i) {
1693     rindex_freon[i] = ppckov[i] * .01095 * 1e9 + 1.2177;
1694   }
1695   // need to be changed 
1696   
1697   // --- Absorbtion lenghts (in cm) 
1698   //      DATA ABSCO_QUARZ / 
1699   //     &    5 * 1000000., 
1700   //     &    29.85,    7.34,     4.134,    1.273,    0.722, 
1701   //     &    0.365, 0.365, 0.365, 0.  / 
1702   // need to be changed 
1703   
1704   // --- Detection efficiencies (quantum efficiency for CsI) 
1705   // --- Define parameters for honeycomb. 
1706   //     Used carbon of equivalent rad. lenght 
1707   
1708   ahon    = 12.01;
1709   zhon    = 6.;
1710   denshon = 2.265;
1711   radlhon = 18.8;
1712   
1713   // --- Parameters to include in GSMIXT, relative to Quarz (SiO2) 
1714   
1715   aqua[0]    = 28.09;
1716   aqua[1]    = 16.;
1717   zqua[0]    = 14.;
1718   zqua[1]    = 8.;
1719   densqua    = 2.64;
1720   nlmatqua   = -2;
1721   wmatqua[0] = 1.;
1722   wmatqua[1] = 2.;
1723   
1724   // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2) 
1725   
1726   aquao[0]    = 28.09;
1727   aquao[1]    = 16.;
1728   zquao[0]    = 14.;
1729   zquao[1]    = 8.;
1730   densquao    = 2.64;
1731   nlmatquao   = -2;
1732   wmatquao[0] = 1.;
1733   wmatquao[1] = 2.;
1734   
1735   // --- Parameters to include in GSMIXT, relative to Freon (C6F14) 
1736   
1737   afre[0]    = 12.;
1738   afre[1]    = 19.;
1739   zfre[0]    = 6.;
1740   zfre[1]    = 9.;
1741   densfre    = 1.7;
1742   nlmatfre   = -2;
1743   wmatfre[0] = 6.;
1744   wmatfre[1] = 14.;
1745   
1746   // --- Parameters to include in GSMIXT, relative to methane (CH4) 
1747   
1748   amet[0]    = 12.01;
1749   amet[1]    = 1.;
1750   zmet[0]    = 6.;
1751   zmet[1]    = 1.;
1752   densmet    = 7.17e-4;
1753   nlmatmet   = -2;
1754   wmatmet[0] = 1.;
1755   wmatmet[1] = 4.;
1756   
1757   // --- Parameters to include in GSMIXT, relative to anode grid (Cu) 
1758   
1759   agri    = 63.54;
1760   zgri    = 29.;
1761   densgri = 8.96;
1762   radlgri = 1.43;
1763   
1764   // --- Parameters to include in GSMATE related to aluminium sheet 
1765   
1766   aal    = 26.98;
1767   zal    = 13.;
1768   densal = 2.7;
1769   radlal = 8.9;
1770   
1771   AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
1772   AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1773   AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1774   AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1775   AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1776   AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1777   AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1778   AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1779   AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1780   AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1781   
1782   tmaxfd = -10.;
1783   stemax = -.1;
1784   deemax = -.2;
1785   epsil  = .001;
1786   stmin  = -.001;
1787   
1788   AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1789   AliMedium(2, "HONEYCOMB$", 6, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1790   AliMedium(3, "QUARZO$", 20, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1791   AliMedium(4, "FREON$", 30, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1792   AliMedium(5, "METANO$", 40, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1793   AliMedium(6, "CSI$", 16, 1, ISXFLD, SXMGMX,tmaxfd, stemax, deemax, epsil, stmin);
1794   AliMedium(7, "GRIGLIA$", 11, 0, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1795   AliMedium(8, "QUARZOO$", 21, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1796   AliMedium(9, "GAP$", 41, 1, ISXFLD, SXMGMX,tmaxfd, .1, -deemax, epsil, -stmin);
1797   AliMedium(10, "ALUMINUM$", 50, 1, ISXFLD, SXMGMX, tmaxfd, stemax, deemax, epsil, stmin);
1798   
1799   
1800   //     Switch on delta-ray production in the methane and freon gaps 
1801   
1802   gMC->Gstpar(idtmed[1002], "LOSS", 1.);
1803   gMC->Gstpar(idtmed[1003], "LOSS", 1.);
1804   gMC->Gstpar(idtmed[1004], "LOSS", 1.);
1805   gMC->Gstpar(idtmed[1008], "LOSS", 1.);
1806   gMC->Gstpar(idtmed[1005], "LOSS", 1.);
1807   gMC->Gstpar(idtmed[1002], "HADR", 1.);
1808   gMC->Gstpar(idtmed[1003], "HADR", 1.);
1809   gMC->Gstpar(idtmed[1004], "HADR", 1.);
1810   gMC->Gstpar(idtmed[1008], "HADR", 1.);
1811   gMC->Gstpar(idtmed[1005], "HADR", 1.);
1812   gMC->Gstpar(idtmed[1002], "DCAY", 1.);
1813   gMC->Gstpar(idtmed[1003], "DCAY", 1.);
1814   gMC->Gstpar(idtmed[1004], "DCAY", 1.);
1815   gMC->Gstpar(idtmed[1008], "DCAY", 1.);
1816   gMC->Gstpar(idtmed[1005], "DCAY", 1.);
1817   geant3->Gsckov(idtmed[1000], 14, ppckov, absco_methane, effic_all, rindex_methane);
1818   geant3->Gsckov(idtmed[1001], 14, ppckov, absco_methane, effic_all, rindex_methane);
1819   geant3->Gsckov(idtmed[1002], 14, ppckov, absco_quarz, effic_all,rindex_quarz);
1820   geant3->Gsckov(idtmed[1003], 14, ppckov, absco_freon, effic_all,rindex_freon);
1821   geant3->Gsckov(idtmed[1004], 14, ppckov, absco_methane, effic_all, rindex_methane);
1822   geant3->Gsckov(idtmed[1005], 14, ppckov, absco_csi, effic_csi, rindex_methane);
1823   geant3->Gsckov(idtmed[1006], 14, ppckov, absco_gri, effic_gri, rindex_gri);
1824   geant3->Gsckov(idtmed[1007], 14, ppckov, absco_quarzo, effic_all, rindex_quarzo);
1825   geant3->Gsckov(idtmed[1008], 14, ppckov, absco_methane, effic_all, rindex_methane);
1826   geant3->Gsckov(idtmed[1009], 14, ppckov, absco_gri, effic_gri, rindex_gri);
1827 }
1828
1829 ClassImp(AliRICHhit)
1830   
1831 //_____________________________________________________________________________
1832 AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol,
1833                        Float_t *hits, Int_t fNpadhits) :
1834     AliHit(shunt,track)
1835 {
1836   //
1837   // Standard constructor for RICH hit
1838   //
1839   Int_t i;
1840   for (i=0;i<2;i++) fVolume[i] = vol[i];
1841   fTrack      = (int)  track;
1842   //Hit information
1843   fPart       = (int)  hits[ 1]; 
1844   //AliHit information, position of the hit
1845   fX          =        hits[ 2];       
1846   fY          =        hits[ 3];
1847   //Pad information
1848   fFirstpad   = fNpadhits;
1849   fLastpad    = fNpadhits-1;
1850   
1851   //Hit information
1852   fModule     = (int)  hits[ 4];
1853   fTheta      =        hits[ 5];
1854   fArrivaltime=        hits[ 6];
1855   fFeed       = (int)  hits[ 7];
1856 }
1857
1858 ClassImp(AliRICHmip)
1859
1860 //_____________________________________________________________________________
1861 AliRICHmip::AliRICHmip(Int_t shunt, Int_t track, Int_t *vol, 
1862                        Float_t *hits, Int_t fNckovs, Int_t fNpadhits) :
1863     AliRICHhit(shunt,track,vol,hits,fNpadhits)
1864 {
1865   //
1866   // Standard constructor for RICH Mip hit
1867   //
1868   fPhi        =        hits[ 8];  
1869   fPs         =        hits[ 9];
1870   fQ          =        hits[10];    
1871   fZ          =        hits[11];
1872   //Ckov information
1873   if ((int) hits[12]){
1874     fFirstCkov  =    fNckovs;
1875     fLastCkov   =    fNckovs-1;
1876   } else {
1877     fFirstCkov  =    -1;
1878     fLastCkov   =    -1;
1879   }
1880 }
1881
1882 ClassImp(AliRICHckov)
1883   
1884 //_____________________________________________________________________________
1885 AliRICHckov::AliRICHckov(Int_t shunt, Int_t track, Int_t *vol, 
1886                          Float_t *hits, Int_t fNmips, Int_t fNpadhits) :
1887   AliRICHhit(shunt,track,vol,hits,fNpadhits)
1888 {
1889   //
1890   // Standard creator for RICH Cherenkov information
1891   //
1892   fEnergy     =        hits[8];
1893   fStop       = (int)  hits[9];
1894   
1895   //Parent info
1896   fParent     =         fNmips;
1897 }
1898
1899 ClassImp(AliRICHpadhit)
1900
1901
1902 //_____________________________________________________________________________
1903 AliRICHpadhit::AliRICHpadhit(Int_t shunt, Int_t track, Int_t *vol,
1904                              Float_t *hits, Int_t fNmips, Int_t fNckovs):
1905   AliHit(shunt,track)
1906 {
1907   //
1908   // Standard constructor for RICH pad hit
1909   //
1910   Int_t i;
1911   for (i=0;i<2;i++) fVolume[i] = vol[i];
1912   
1913   // Hit information
1914   fX          = (int)  hits[ 2];
1915   fY          = (int)  hits[ 3];
1916   fModule     = (int)  hits[ 4];
1917   fParentMip  =        fNmips;
1918   fParentCkov =        fNckovs;
1919   fProd       = (int)  hits[ 5];
1920   fCharge     =        hits[ 6];
1921   fZ          =          -999.0; // Not implemented
1922 }
1923
1924