RICH/AliRICH.cxx

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