RICH/AliRICH.cxx

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