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