RICH/AliRICH.cxx

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