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