RICH/AliRICH.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17   $Log$
  18   Revision 1.47  2001/03/14 18:13:56  jbarbosa
  19   Several changes to adapt to new IO.
  20   Removed digitising function, using AliRICHMerger::Digitise from now on.
  21
  22   Revision 1.46  2001/03/12 17:46:33  hristov
  23   Changes needed on Sun with CC 5.0
  24
  25   Revision 1.45  2001/02/27 22:11:46  jbarbosa
  26   Testing TreeS, removing of output.
  27
  28   Revision 1.44  2001/02/27 15:19:12  jbarbosa
  29   Transition to SDigits.
  30
  31   Revision 1.43  2001/02/23 17:19:06  jbarbosa
  32   Corrected photocathode definition in BuildGeometry().
  33
  34   Revision 1.42  2001/02/13 20:07:23  jbarbosa
  35   Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum
  36   when entering the freon. Corrected calls to particle stack.
  37
  38   Revision 1.41  2001/01/26 20:00:20  hristov
  39   Major upgrade of AliRoot code
  40
  41   Revision 1.40  2001/01/24 20:58:03  jbarbosa
  42   Enhanced BuildGeometry. Now the photocathodes are drawn.
  43
  44   Revision 1.39  2001/01/22 21:40:24  jbarbosa
  45   Removing magic numbers
  46
  47   Revision 1.37  2000/12/20 14:07:25  jbarbosa
  48   Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
  49
  50   Revision 1.36  2000/12/18 17:45:54  jbarbosa
  51   Cleaned up PadHits object.
  52
  53   Revision 1.35  2000/12/15 16:49:40  jbarbosa
  54   Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
  55
  56   Revision 1.34  2000/11/10 18:12:12  jbarbosa
  57   Bug fix for AliRICHCerenkov (thanks to P. Hristov)
  58
  59   Revision 1.33  2000/11/02 10:09:01  jbarbosa
  60   Minor bug correction (some pointers were not initialised in the default constructor)
  61
  62   Revision 1.32  2000/11/01 15:32:55  jbarbosa
  63   Updated to handle both reconstruction algorithms.
  64
  65   Revision 1.31  2000/10/26 20:18:33  jbarbosa
  66   Supports for methane and freon vessels
  67
  68   Revision 1.30  2000/10/24 13:19:12  jbarbosa
  69   Geometry updates.
  70
  71   Revision 1.29  2000/10/19 19:39:25  jbarbosa
  72   Some more changes to geometry. Further correction of digitisation "per part. type"
  73
  74   Revision 1.28  2000/10/17 20:50:57  jbarbosa
  75   Inversed digtise by particle type (now, only the selected particle type is not digitsed).
  76   Corrected several geometry minor bugs.
  77   Added new parameter (opaque quartz thickness).
  78
  79   Revision 1.27  2000/10/11 10:33:55  jbarbosa
  80   Corrected bug introduced by earlier revisions  (CerenkovData array cannot be reset to zero on wach call of StepManager)
  81
  82   Revision 1.26  2000/10/03 21:44:08  morsch
  83   Use AliSegmentation and AliHit abstract base classes.
  84
  85   Revision 1.25  2000/10/02 21:28:12  fca
  86   Removal of useless dependecies via forward declarations
  87
  88   Revision 1.24  2000/10/02 15:43:17  jbarbosa
  89   Fixed forward declarations.
  90   Fixed honeycomb density.
  91   Fixed cerenkov storing.
  92   New electronics.
  93
  94   Revision 1.23  2000/09/13 10:42:14  hristov
  95   Minor corrections for HP, DEC and Sun; strings.h included
  96
  97   Revision 1.22  2000/09/12 18:11:13  fca
  98   zero hits area before using
  99
 100   Revision 1.21  2000/07/21 10:21:07  morsch
 101   fNrawch   = 0; and  fNrechits = 0; in the default constructor.
 102
 103   Revision 1.20  2000/07/10 15:28:39  fca
 104   Correction of the inheritance scheme
 105
 106   Revision 1.19  2000/06/30 16:29:51  dibari
 107   Added kDebugLevel variable to control output size on demand
 108
 109   Revision 1.18  2000/06/12 15:15:46  jbarbosa
 110   Cleaned up version.
 111
 112   Revision 1.17  2000/06/09 14:58:37  jbarbosa
 113   New digitisation per particle type
 114
 115   Revision 1.16  2000/04/19 12:55:43  morsch
 116   Newly structured and updated version (JB, AM)
 117
 118 */
 119
 120
 121 ////////////////////////////////////////////////
 122 //  Manager and hits classes for set:RICH     //
 123 ////////////////////////////////////////////////
 124
 125 #include <TBRIK.h>
 126 #include <TTUBE.h>
 127 #include <TNode.h>
 128 #include <TRandom.h>
 129 #include <TObject.h>
 130 #include <TVector.h>
 131 #include <TObjArray.h>
 132 #include <TArrayF.h>
 133 #include <TFile.h>
 134 #include <TParticle.h>
 135 #include <TGeometry.h>
 136 #include <TTree.h>
 137
 138 #include <iostream.h>
 139 #include <strings.h>
 140
 141 #include "AliRICH.h"
 142 #include "AliSegmentation.h"
 143 #include "AliRICHSegmentationV0.h"
 144 #include "AliRICHHit.h"
 145 #include "AliRICHCerenkov.h"
 146 #include "AliRICHSDigit.h"
 147 #include "AliRICHDigit.h"
 148 #include "AliRICHTransientDigit.h"
 149 #include "AliRICHRawCluster.h"
 150 #include "AliRICHRecHit1D.h"
 151 #include "AliRICHRecHit3D.h"
 152 #include "AliRICHHitMapA1.h"
 153 #include "AliRICHClusterFinder.h"
 154 #include "AliRICHMerger.h"
 155 #include "AliRun.h"
 156 #include "AliMC.h"
 157 #include "AliMagF.h"
 158 #include "AliConst.h"
 159 #include "AliPDG.h"
 160 #include "AliPoints.h"
 161 #include "AliCallf77.h"
 162
 163
 164 // Static variables for the pad-hit iterator routines
 165 static Int_t sMaxIterPad=0;
 166 static Int_t sCurIterPad=0;
 167
 168 ClassImp(AliRICH)
 169
 170 //___________________________________________
 171 AliRICH::AliRICH()
 172 {
 173 // Default constructor for RICH manager class
 174
 175     fIshunt     = 0;
 176     fHits       = 0;
 177     fSDigits    = 0;
 178     fNSDigits   = 0;
 179     fNcerenkovs = 0;
 180     fDchambers  = 0;
 181     fRecHits1D = 0;
 182     fRecHits3D = 0;
 183     fRawClusters = 0;
 184     fChambers = 0;
 185     fCerenkovs  = 0;
 186     for (Int_t i=0; i<7; i++)
 187       {
 188         fNdch[i]       = 0;
 189         fNrawch[i]   = 0;
 190         fNrechits1D[i] = 0;
 191         fNrechits3D[i] = 0;
 192       }
 193
 194     fFileName = 0;
 195 }
 196
 197 //___________________________________________
 198 AliRICH::AliRICH(const char *name, const char *title)
 199     : AliDetector(name,title)
 200 {
 201 //Begin_Html
 202 /*
 203   <img src="gif/alirich.gif">
 204 */
 205 //End_Html
 206
 207     fHits       = new TClonesArray("AliRICHHit",1000  );
 208     gAlice->AddHitList(fHits);
 209     fSDigits    = new TClonesArray("AliRICHSDigit",100000);
 210     fCerenkovs  = new TClonesArray("AliRICHCerenkov",1000);
 211     gAlice->AddHitList(fCerenkovs);
 212     //gAlice->AddHitList(fHits);
 213     fNSDigits   = 0;
 214     fNcerenkovs = 0;
 215     fIshunt     = 0;
 216
 217     //fNdch      = new Int_t[kNCH];
 218
 219     fDchambers = new TObjArray(kNCH);
 220
 221     fRecHits1D = new TObjArray(kNCH);
 222     fRecHits3D = new TObjArray(kNCH);
 223
 224     Int_t i;
 225
 226     for (i=0; i<kNCH ;i++) {
 227         (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
 228         fNdch[i]=0;
 229     }
 230
 231     //fNrawch      = new Int_t[kNCH];
 232
 233     fRawClusters = new TObjArray(kNCH);
 234     //printf("Created fRwClusters with adress:%p",fRawClusters);
 235
 236     for (i=0; i<kNCH ;i++) {
 237       (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
 238       fNrawch[i]=0;
 239     }
 240
 241     //fNrechits      = new Int_t[kNCH];
 242
 243     for (i=0; i<kNCH ;i++) {
 244       (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
 245     }
 246     for (i=0; i<kNCH ;i++) {
 247       (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
 248     }
 249     //printf("Created fRecHits with adress:%p",fRecHits);
 250
 251
 252     SetMarkerColor(kRed);
 253
 254     /*fChambers = new TObjArray(kNCH);
 255     for (i=0; i<kNCH; i++)
 256       (*fChambers)[i] = new AliRICHChamber();*/
 257
 258     fFileName = 0;
 259 }
 260
 261 AliRICH::AliRICH(const AliRICH& RICH)
 262 {
 263 // Copy Constructor
 264 }
 265
 266
 267 //___________________________________________
 268 AliRICH::~AliRICH()
 269 {
 270
 271 // Destructor of RICH manager class
 272
 273     fIshunt  = 0;
 274     delete fHits;
 275     delete fSDigits;
 276     delete fCerenkovs;
 277
 278     //PH Delete TObjArrays
 279     if (fChambers) {
 280       fChambers->Delete();
 281       delete fChambers;
 282     }
 283     if (fDchambers) {
 284       fDchambers->Delete();
 285       delete fDchambers;
 286     }
 287     if (fRawClusters) {
 288       fRawClusters->Delete();
 289       delete fRawClusters;
 290     }
 291     if (fRecHits1D) {
 292       fRecHits1D->Delete();
 293       delete fRecHits1D;
 294     }
 295     if (fRecHits3D) {
 296       fRecHits3D->Delete();
 297       delete fRecHits3D;
 298     }
 299
 300 }
 301
 302
 303 //_____________________________________________________________________________
 304 Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
 305 {
 306 //
 307 //  Calls the charge disintegration method of the current chamber and adds
 308 //  the simulated cluster to the root treee
 309 //
 310     Int_t clhits[5];
 311     Float_t newclust[4][500];
 312     Int_t nnew;
 313
 314 //
 315 //  Integrated pulse height on chamber
 316
 317     clhits[0]=fNhits+1;
 318
 319     ((AliRICHChamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
 320     Int_t ic=0;
 321
 322 //
 323 //  Add new clusters
 324     for (Int_t i=0; i<nnew; i++) {
 325         if (Int_t(newclust[0][i]) > 0) {
 326             ic++;
 327 //  Cluster Charge
 328             clhits[1] = Int_t(newclust[0][i]);
 329 //  Pad: ix
 330             clhits[2] = Int_t(newclust[1][i]);
 331 //  Pad: iy
 332             clhits[3] = Int_t(newclust[2][i]);
 333 //  Pad: chamber sector
 334             clhits[4] = Int_t(newclust[3][i]);
 335
 336             //printf(" %d %d %d %d %d\n",  clhits[0],  clhits[1],  clhits[2],  clhits[3],  clhits[4]);
 337
 338             AddSDigit(clhits);
 339         }
 340     }
 341
 342     if (gAlice->TreeS())
 343       {
 344         gAlice->TreeS()->Fill();
 345         gAlice->TreeS()->Write(0,TObject::kOverwrite);
 346         //printf("Filled SDigits...\n");
 347       }
 348
 349 return nnew;
 350 }
 351 //___________________________________________
 352 void AliRICH::Hits2SDigits()
 353 {
 354
 355 // Dummy: sdigits are created during transport.
 356 // Called from alirun.
 357
 358   int nparticles = gAlice->GetNtrack();
 359   cout << "Particles (RICH):" <<nparticles<<endl;
 360   if (nparticles > 0) printf("SDigits were already generated.\n");
 361
 362 }
 363
 364 //___________________________________________
 365 void AliRICH::SDigits2Digits(Int_t nev, Int_t flag)
 366 {
 367
 368 //
 369 // Generate digits.
 370 // Called from macro. Multiple events, more functionality.
 371
 372   AliRICHChamber*       iChamber;
 373
 374   printf("Generating tresholds...\n");
 375
 376   for(Int_t i=0;i<7;i++) {
 377     iChamber = &(Chamber(i));
 378     iChamber->GenerateTresholds();
 379   }
 380
 381   int nparticles = gAlice->GetNtrack();
 382   if (nparticles > 0)
 383     {
 384       if (fMerger) {
 385         fMerger->Init();
 386         fMerger->Digitise(nev,flag);
 387       }
 388     }
 389   //Digitise(nev,flag);
 390 }
 391 //___________________________________________
 392 void AliRICH::SDigits2Digits()
 393 {
 394
 395 //
 396 // Generate digits
 397 // Called from alirun, single event only.
 398
 399   AliRICHChamber*       iChamber;
 400
 401   printf("Generating tresholds...\n");
 402
 403   for(Int_t i=0;i<7;i++) {
 404     iChamber = &(Chamber(i));
 405     iChamber->GenerateTresholds();
 406   }
 407
 408   int nparticles = gAlice->GetNtrack();
 409   cout << "Particles (RICH):" <<nparticles<<endl;
 410   if (nparticles > 0)
 411     {
 412       if (fMerger) {
 413         fMerger->Init();
 414         fMerger->Digitise(0,0);
 415       }
 416     }
 417 }
 418 //___________________________________________
 419 void AliRICH::Digits2Reco()
 420 {
 421
 422 // Generate clusters
 423 // Called from alirun, single event only.
 424
 425   int nparticles = gAlice->GetNtrack();
 426   cout << "Particles (RICH):" <<nparticles<<endl;
 427   if (nparticles > 0) FindClusters(0,0);
 428
 429 }
 430
 431 //___________________________________________
 432 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 433 {
 434
 435 //
 436 // Adds a hit to the Hits list
 437 //
 438
 439     TClonesArray &lhits = *fHits;
 440     new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
 441 }
 442 //_____________________________________________________________________________
 443 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
 444 {
 445
 446 //
 447 // Adds a RICH cerenkov hit to the Cerenkov Hits list
 448 //
 449
 450     TClonesArray &lcerenkovs = *fCerenkovs;
 451     new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
 452     //printf ("Done for Cerenkov %d\n\n\n\n",fNcerenkovs);
 453 }
 454 //___________________________________________
 455 void AliRICH::AddSDigit(Int_t *clhits)
 456 {
 457
 458 //
 459 // Add a RICH pad hit to the list
 460 //
 461
 462   //printf("fsdigits:%p, data: %d\n",fSDigits,clhits[2]);
 463   TClonesArray &lSDigits = *fSDigits;
 464   new(lSDigits[fNSDigits++]) AliRICHSDigit(clhits);
 465 }
 466 //_____________________________________________________________________________
 467 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
 468 {
 469
 470   //
 471   // Add a RICH digit to the list
 472   //
 473
 474   //printf("fdigits:%p, data: %d\n",((TClonesArray*)(*fDchambers)[id]),digits[0]);
 475   TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
 476   new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
 477 }
 478
 479 //_____________________________________________________________________________
 480 void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
 481 {
 482     //
 483     // Add a RICH digit to the list
 484     //
 485
 486     TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
 487     new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
 488 }
 489
 490 //_____________________________________________________________________________
 491 void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
 492 {
 493
 494   //
 495   // Add a RICH reconstructed hit to the list
 496   //
 497
 498     TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
 499     new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
 500 }
 501
 502 //_____________________________________________________________________________
 503 void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit)
 504 {
 505
 506   //
 507   // Add a RICH reconstructed hit to the list
 508   //
 509
 510     TClonesArray &lrec3D = *((TClonesArray*)(*fRecHits3D)[id]);
 511     new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit);
 512 }
 513
 514 //___________________________________________
 515 void AliRICH::BuildGeometry()
 516
 517 {
 518
 519   //
 520   // Builds a TNode geometry for event display
 521   //
 522     TNode *node, *subnode, *top;
 523
 524     const int kColorRICH = kRed;
 525     //
 526     top=gAlice->GetGeometry()->GetNode("alice");
 527
 528     AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
 529     AliRICHSegmentationV0*  segmentation;
 530     AliRICHChamber*       iChamber;
 531
 532     iChamber = &(pRICH->Chamber(0));
 533     segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
 534
 535     new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
 536
 537     Float_t padplane_width = segmentation->GetPadPlaneWidth();
 538     Float_t padplane_length = segmentation->GetPadPlaneLength();
 539
 540     //printf("\n\n\n\n\n In BuildGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f\n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy());
 541
 542     new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2);
 543
 544     //printf("\n\n\n\n\n Padplane   w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2);
 545     //printf("\n\n\n\n\n Padplane   w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength());
 546
 547
 548     top->cd();
 549     Float_t pos1[3]={0,471.8999,165.2599};
 550     //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
 551     new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
 552     node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
 553     node->SetLineColor(kColorRICH);
 554     node->cd();
 555     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 556     subnode->SetLineColor(kGreen);
 557     fNodes->Add(subnode);
 558     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 559     subnode->SetLineColor(kGreen);
 560     fNodes->Add(subnode);
 561     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 562     subnode->SetLineColor(kGreen);
 563     fNodes->Add(subnode);
 564     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 565     subnode->SetLineColor(kGreen);
 566     fNodes->Add(subnode);
 567     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 568     subnode->SetLineColor(kGreen);
 569     fNodes->Add(subnode);
 570     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 571     subnode->SetLineColor(kGreen);
 572     fNodes->Add(subnode);
 573     fNodes->Add(node);
 574
 575
 576     top->cd();
 577     Float_t pos2[3]={171,470,0};
 578     //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
 579     new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
 580     node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
 581     node->SetLineColor(kColorRICH);
 582     node->cd();
 583     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 584     subnode->SetLineColor(kGreen);
 585     fNodes->Add(subnode);
 586     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 587     subnode->SetLineColor(kGreen);
 588     fNodes->Add(subnode);
 589     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 590     subnode->SetLineColor(kGreen);
 591     fNodes->Add(subnode);
 592     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 593     subnode->SetLineColor(kGreen);
 594     fNodes->Add(subnode);
 595     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 596     subnode->SetLineColor(kGreen);
 597     fNodes->Add(subnode);
 598     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 599     subnode->SetLineColor(kGreen);
 600     fNodes->Add(subnode);
 601     fNodes->Add(node);
 602
 603
 604     top->cd();
 605     Float_t pos3[3]={0,500,0};
 606     //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
 607     new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
 608     node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
 609     node->SetLineColor(kColorRICH);
 610     node->cd();
 611     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 612     subnode->SetLineColor(kGreen);
 613     fNodes->Add(subnode);
 614     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 615     subnode->SetLineColor(kGreen);
 616     fNodes->Add(subnode);
 617     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 618     subnode->SetLineColor(kGreen);
 619     fNodes->Add(subnode);
 620     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 621     subnode->SetLineColor(kGreen);
 622     fNodes->Add(subnode);
 623     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 624     subnode->SetLineColor(kGreen);
 625     fNodes->Add(subnode);
 626     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 627     subnode->SetLineColor(kGreen);
 628     fNodes->Add(subnode);
 629     fNodes->Add(node);
 630
 631     top->cd();
 632     Float_t pos4[3]={-171,470,0};
 633     //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
 634     new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
 635     node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
 636     node->SetLineColor(kColorRICH);
 637     node->cd();
 638     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 639     subnode->SetLineColor(kGreen);
 640     fNodes->Add(subnode);
 641     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 642     subnode->SetLineColor(kGreen);
 643     fNodes->Add(subnode);
 644     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 645     subnode->SetLineColor(kGreen);
 646     fNodes->Add(subnode);
 647     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 648     subnode->SetLineColor(kGreen);
 649     fNodes->Add(subnode);
 650     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 651     subnode->SetLineColor(kGreen);
 652     fNodes->Add(subnode);
 653     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 654     subnode->SetLineColor(kGreen);
 655     fNodes->Add(subnode);
 656     fNodes->Add(node);
 657
 658
 659     top->cd();
 660     Float_t pos5[3]={161.3999,443.3999,-165.3};
 661     //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
 662     new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
 663     node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
 664     node->SetLineColor(kColorRICH);
 665     node->cd();
 666     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 667     subnode->SetLineColor(kGreen);
 668     fNodes->Add(subnode);
 669     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 670     subnode->SetLineColor(kGreen);
 671     fNodes->Add(subnode);
 672     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 673     subnode->SetLineColor(kGreen);
 674     fNodes->Add(subnode);
 675     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 676     subnode->SetLineColor(kGreen);
 677     fNodes->Add(subnode);
 678     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 679     subnode->SetLineColor(kGreen);
 680     fNodes->Add(subnode);
 681     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 682     subnode->SetLineColor(kGreen);
 683     fNodes->Add(subnode);
 684     fNodes->Add(node);
 685
 686
 687     top->cd();
 688     Float_t pos6[3]={0., 471.9, -165.3,};
 689     //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
 690     new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
 691     node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
 692     node->SetLineColor(kColorRICH);
 693
 694     fNodes->Add(node);node->cd();
 695     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 696     subnode->SetLineColor(kGreen);
 697     fNodes->Add(subnode);
 698     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 699     subnode->SetLineColor(kGreen);
 700     fNodes->Add(subnode);
 701     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 702     subnode->SetLineColor(kGreen);
 703     fNodes->Add(subnode);
 704     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 705     subnode->SetLineColor(kGreen);
 706     fNodes->Add(subnode);
 707     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 708     subnode->SetLineColor(kGreen);
 709     fNodes->Add(subnode);
 710     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 711     subnode->SetLineColor(kGreen);
 712     fNodes->Add(subnode);
 713
 714
 715     top->cd();
 716     Float_t pos7[3]={-161.399,443.3999,-165.3};
 717     //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
 718     new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
 719     node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
 720     node->SetLineColor(kColorRICH);
 721     node->cd();
 722     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 723     subnode->SetLineColor(kGreen);
 724     fNodes->Add(subnode);
 725     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,padplane_length/2 + segmentation->DeadZone()/2,"");
 726     subnode->SetLineColor(kGreen);
 727     fNodes->Add(subnode);
 728     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,padplane_length/2 + segmentation->DeadZone()/2,"");
 729     subnode->SetLineColor(kGreen);
 730     fNodes->Add(subnode);
 731     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 732     subnode->SetLineColor(kGreen);
 733     fNodes->Add(subnode);
 734     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 735     subnode->SetLineColor(kGreen);
 736     fNodes->Add(subnode);
 737     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),0,-padplane_length/2 - segmentation->DeadZone()/2,"");
 738     subnode->SetLineColor(kGreen);
 739     fNodes->Add(subnode);
 740     fNodes->Add(node);
 741
 742 }
 743
 744 //___________________________________________
 745 void AliRICH::CreateGeometry()
 746 {
 747     //
 748     // Create the geometry for RICH version 1
 749     //
 750     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
 751     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
 752     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
 753     //
 754     //Begin_Html
 755     /*
 756       <img src="picts/AliRICHv1.gif">
 757     */
 758     //End_Html
 759     //Begin_Html
 760     /*
 761       <img src="picts/AliRICHv1Tree.gif">
 762     */
 763     //End_Html
 764
 765   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
 766   AliRICHSegmentationV0*  segmentation;
 767   AliRICHGeometry*  geometry;
 768   AliRICHChamber*       iChamber;
 769
 770   iChamber = &(pRICH->Chamber(0));
 771   segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
 772   geometry=iChamber->GetGeometryModel();
 773
 774   Float_t distance;
 775   distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
 776   geometry->SetRadiatorToPads(distance);
 777
 778   //Opaque quartz thickness
 779   Float_t oqua_thickness = .5;
 780   //CsI dimensions
 781
 782   //Float_t csi_length = 160*.8 + 2.6;
 783   //Float_t csi_width = 144*.84 + 2*2.6;
 784
 785   Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
 786   Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
 787
 788   //printf("\n\n\n\n\n In CreateGeometry() npx: %d, npy: %d, dpx: %f, dpy:%f  deadzone: %f \n\n\n\n\n\n",segmentation->Npx(),segmentation->Npy(),segmentation->Dpx(),segmentation->Dpy(),segmentation->DeadZone());
 789
 790   Int_t *idtmed = fIdtmed->GetArray()-999;
 791
 792     Int_t i;
 793     Float_t zs;
 794     Int_t idrotm[1099];
 795     Float_t par[3];
 796
 797     // --- Define the RICH detector
 798     //     External aluminium box
 799     par[0] = 68.8;
 800     par[1] = 13;                 //Original Settings
 801     par[2] = 70.86;
 802     /*par[0] = 73.15;
 803     par[1] = 11.5;
 804     par[2] = 71.1;*/
 805     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 806
 807     //     Air
 808     par[0] = 66.3;
 809     par[1] = 13;                 //Original Settings
 810     par[2] = 68.35;
 811     /*par[0] = 66.55;
 812     par[1] = 11.5;
 813     par[2] = 64.8;*/
 814     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
 815
 816     //    Air 2 (cutting the lower part of the box)
 817
 818     par[0] = 1.25;
 819     par[1] = 3;                 //Original Settings
 820     par[2] = 70.86;
 821     gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
 822
 823     //    Air 3 (cutting the lower part of the box)
 824
 825     par[0] = 66.3;
 826     par[1] = 3;                 //Original Settings
 827     par[2] = 1.2505;
 828     gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
 829
 830     //     Honeycomb
 831     par[0] = 66.3;
 832     par[1] = .188;                 //Original Settings
 833     par[2] = 68.35;
 834     /*par[0] = 66.55;
 835     par[1] = .188;
 836     par[2] = 63.1;*/
 837     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
 838
 839     //     Aluminium sheet
 840     par[0] = 66.3;
 841     par[1] = .025;                 //Original Settings
 842     par[2] = 68.35;
 843     /*par[0] = 66.5;
 844     par[1] = .025;
 845     par[2] = 63.1;*/
 846     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
 847
 848     //     Quartz
 849     par[0] = geometry->GetQuartzWidth()/2;
 850     par[1] = geometry->GetQuartzThickness()/2;
 851     par[2] = geometry->GetQuartzLength()/2;
 852     /*par[0] = 63.1;
 853     par[1] = .25;                  //Original Settings
 854     par[2] = 65.5;*/
 855     /*par[0] = geometry->GetQuartzWidth()/2;
 856     par[1] = geometry->GetQuartzThickness()/2;
 857     par[2] = geometry->GetQuartzLength()/2;*/
 858     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f %f %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[0],par[1],par[2]);
 859     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
 860
 861     //     Spacers (cylinders)
 862     par[0] = 0.;
 863     par[1] = .5;
 864     par[2] = geometry->GetFreonThickness()/2;
 865     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
 866
 867     //     Feet (freon slabs supports)
 868
 869     par[0] = .7;
 870     par[1] = .3;
 871     par[2] = 1.9;
 872     gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
 873
 874     //     Opaque quartz
 875     par[0] = geometry->GetQuartzWidth()/2;
 876     par[1] = .2;
 877     par[2] = geometry->GetQuartzLength()/2;
 878     /*par[0] = 61.95;
 879     par[1] = .2;                   //Original Settings
 880     par[2] = 66.5;*/
 881     /*par[0] = 66.5;
 882     par[1] = .2;
 883     par[2] = 61.95;*/
 884     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
 885
 886     //     Frame of opaque quartz
 887     par[0] = geometry->GetOuterFreonWidth()/2;
 888     //+ oqua_thickness;
 889     par[1] = geometry->GetFreonThickness()/2;
 890     par[2] = geometry->GetOuterFreonLength()/2;
 891     //+ oqua_thickness;
 892     /*par[0] = 20.65;
 893     par[1] = .5;                   //Original Settings
 894     par[2] = 66.5;*/
 895     /*par[0] = 66.5;
 896     par[1] = .5;
 897     par[2] = 20.65;*/
 898     gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
 899
 900     par[0] = geometry->GetInnerFreonWidth()/2;
 901     par[1] = geometry->GetFreonThickness()/2;
 902     par[2] = geometry->GetInnerFreonLength()/2;
 903     gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
 904
 905     //     Little bar of opaque quartz
 906     //par[0] = .275;
 907     //par[1] = geometry->GetQuartzThickness()/2;
 908     //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
 909     //par[2] = geometry->GetInnerFreonLength()/2;
 910     //+ oqua_thickness;
 911     /*par[0] = .275;
 912     par[1] = .25;                   //Original Settings
 913     par[2] = 63.1;*/
 914     /*par[0] = 63.1;
 915     par[1] = .25;
 916     par[2] = .275;*/
 917     //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
 918
 919     //     Freon
 920     par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
 921     par[1] = geometry->GetFreonThickness()/2;
 922     par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
 923     /*par[0] = 20.15;
 924     par[1] = .5;                   //Original Settings
 925     par[2] = 65.5;*/
 926     /*par[0] = 65.5;
 927     par[1] = .5;
 928     par[2] = 20.15;*/
 929     gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
 930
 931     par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
 932     par[1] = geometry->GetFreonThickness()/2;
 933     par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
 934     gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
 935
 936     //     Methane
 937     //par[0] = 64.8;
 938     par[0] = csi_width/2;
 939     par[1] = geometry->GetGapThickness()/2;
 940     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
 941     //par[2] = 64.8;
 942     par[2] = csi_length/2;
 943     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
 944
 945     //     Methane gap
 946     //par[0] = 64.8;
 947     par[0] = csi_width/2;
 948     par[1] = geometry->GetProximityGapThickness()/2;
 949     //printf("\n\n\n\n\n\n\n\\n\n\n\n Gap Thickness: %f\n\n\n\n\n\n\n\n\n\n\n\n\n\n",par[1]);
 950     //par[2] = 64.8;
 951     par[2] = csi_length/2;
 952     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
 953
 954     //     CsI photocathode
 955     //par[0] = 64.8;
 956     par[0] = csi_width/2;
 957     par[1] = .25;
 958     //par[2] = 64.8;
 959     par[2] = csi_length/2;
 960     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
 961
 962     //     Anode grid
 963     par[0] = 0.;
 964     par[1] = .001;
 965     par[2] = 20.;
 966     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
 967
 968     // Wire supports
 969     // Bar of metal
 970
 971     par[0] = csi_width/2;
 972     par[1] = 1.05;
 973     par[2] = 1.05;
 974     gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
 975
 976     // Ceramic pick up (base)
 977
 978     par[0] =  csi_width/2;
 979     par[1] = .25;
 980     par[2] = 1.05;
 981     gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
 982
 983     // Ceramic pick up (head)
 984
 985     par[0] = csi_width/2;
 986     par[1] = .1;
 987     par[2] = .1;
 988     gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
 989
 990     // Aluminium supports for methane and CsI
 991     // Short bar
 992
 993     par[0] = csi_width/2;
 994     par[1] = geometry->GetGapThickness()/2 + .25;
 995     par[2] = (68.35 - csi_length/2)/2;
 996     gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
 997
 998     // Long bar
 999
1000     par[0] = (66.3 - csi_width/2)/2;
1001     par[1] = geometry->GetGapThickness()/2 + .25;
1002     par[2] = csi_length/2 + 68.35 - csi_length/2;
1003     gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
1004
1005     // Aluminium supports for freon
1006     // Short bar
1007
1008     par[0] = geometry->GetQuartzWidth()/2;
1009     par[1] = .3;
1010     par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
1011     gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
1012
1013     // Long bar
1014
1015     par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
1016     par[1] = .3;
1017     par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
1018     gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
1019
1020     // PCB backplane
1021
1022     par[0] = csi_width/2;
1023     par[1] = .25;
1024     par[2] = csi_length/4 -.5025;
1025     gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
1026
1027
1028     // Backplane supports
1029
1030     // Aluminium slab
1031
1032     par[0] = 33.15;
1033     par[1] = 2;
1034     par[2] = 21.65;
1035     gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
1036
1037     // Big hole
1038
1039     par[0] = 9.05;
1040     par[1] = 2;
1041     par[2] = 4.4625;
1042     gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
1043
1044     // Small hole
1045
1046     par[0] = 5.7;
1047     par[1] = 2;
1048     par[2] = 4.4625;
1049     gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
1050
1051     // Place holes inside backplane support
1052
1053     gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
1054     gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
1055     gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
1056     gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
1057     gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1058     gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1059     gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1060     gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1061     gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
1062     gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
1063     gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1064     gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1065     gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1066     gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1067     gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1068     gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1069
1070
1071
1072     // --- Places the detectors defined with GSVOLU
1073     //     Place material inside RICH
1074     gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
1075     gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
1076     gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, 0., 0, "ONLY");
1077     gMC->Gspos("AIR3", 1, "RICH", 0.,  1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5, -68.35 - 1.25, 0, "ONLY");
1078     gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.5,  68.35 + 1.25, 0, "ONLY");
1079
1080
1081     gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
1082     gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2  - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
1083     gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
1084     gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1085     gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
1086     gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1087     gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1088     gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1089     gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1090     gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1091     gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1092     gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
1093
1094     // Supports placing
1095
1096     // Methane supports
1097     gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1098     gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1099     gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
1100     gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
1101
1102     //Freon supports
1103
1104     Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
1105
1106     gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1107     gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1108     gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1109     gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1110
1111     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1112
1113      //Placing of the spacers inside the freon slabs
1114
1115     Int_t nspacers = 30;
1116     //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n Spacers:%d\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n",nspacers);
1117
1118     //printf("Nspacers: %d", nspacers);
1119
1120     for (i = 0; i < nspacers/3; i++) {
1121         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1122         gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
1123     }
1124
1125     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1126         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1127         gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
1128     }
1129
1130     for (i = (nspacers*2)/3; i < nspacers; ++i) {
1131         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1132         gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1133     }
1134
1135     for (i = 0; i < nspacers/3; i++) {
1136         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1137         gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
1138     }
1139
1140     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1141         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1142         gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
1143     }
1144
1145     for (i = (nspacers*2)/3; i < nspacers; ++i) {
1146         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1147         gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1148     }
1149
1150
1151     gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
1152     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
1153     gMC->Gspos("OQF1", 1, "SRIC", geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2 + 2, 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
1154 //    printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2);
1155     gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");          //Original settings
1156     gMC->Gspos("OQF1", 3, "SRIC", - (geometry->GetOuterFreonWidth()/2 + geometry->GetInnerFreonWidth()/2) - 2, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");       //Original settings (-31.3)
1157     //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY");           //Original settings (-21.65)
1158     //gMC->Gspos("BARR", 2, "QUAR",  geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY");            //Original settings (21.65)
1159     gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
1160     gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
1161     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1162     gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
1163
1164     // Wire support placing
1165
1166     gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
1167     gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
1168     gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
1169
1170     // Backplane placing
1171
1172     gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
1173     gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
1174     gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1175     gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1176     gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1177     gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1178
1179     // PCB placing
1180
1181     gMC->Gspos("PCB ", 1, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 1.2, 0, "ONLY");
1182     gMC->Gspos("PCB ", 2, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 1.2, 0, "ONLY");
1183
1184
1185
1186     //printf("Position of the gap: %f to %f\n", 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - .2, 1.276 + geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 + .2);
1187
1188     //     Place RICH inside ALICE apparatus
1189
1190     // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm)
1191
1192     //Float_t offset =  1.276 - geometry->GetGapThickness()/2;
1193
1194     AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1195     AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1196     AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1197     AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1198     AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1199     AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1200     AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1201
1202     gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
1203     gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
1204     gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
1205     gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
1206     gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
1207     gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
1208     gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");
1209
1210 }
1211
1212
1213 //___________________________________________
1214 void AliRICH::CreateMaterials()
1215 {
1216     //
1217     // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1218     // ORIGIN    : NICK VAN EIJNDHOVEN
1219     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1220     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
1221     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1222     //
1223     Int_t   isxfld = gAlice->Field()->Integ();
1224     Float_t sxmgmx = gAlice->Field()->Max();
1225     Int_t i;
1226
1227     /************************************Antonnelo's Values (14-vectors)*****************************************/
1228     /*
1229     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,
1230                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1231     Float_t rIndexQuarz[14] = { 1.528309,1.533333,
1232                                  1.538243,1.544223,1.550568,1.55777,
1233                                  1.565463,1.574765,1.584831,1.597027,
1234                                1.611858,1.6277,1.6472,1.6724 };
1235     Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1236     Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1237     Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1238     Float_t abscoFreon[14] = { 179.0987,179.0987,
1239                                 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
1240     //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1241         //                       1e-5,1e-5,1e-5,1e-5,1e-5 };
1242     Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
1243                                 14.177,9.282,4.0925,1.149,.3627,.10857 };
1244     Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1245                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
1246     Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1247                               1e-4,1e-4,1e-4,1e-4 };
1248     Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1249                                   1e6,1e6,1e6 };
1250     Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1251                               1e-4,1e-4,1e-4,1e-4 };
1252     Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1253     Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
1254                               .17425,.1785,.1836,.1904,.1938,.221 };
1255     Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1256     */
1257
1258
1259     /**********************************End of Antonnelo's Values**********************************/
1260
1261     /**********************************Values from rich_media.f (31-vectors)**********************************/
1262
1263
1264     //Photons energy intervals
1265     Float_t ppckov[26];
1266     for (i=0;i<26;i++)
1267     {
1268         ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
1269         //printf ("Energy intervals: %e\n",ppckov[i]);
1270     }
1271
1272
1273     //Refraction index for quarz
1274     Float_t rIndexQuarz[26];
1275     Float_t  e1= 10.666;
1276     Float_t  e2= 18.125;
1277     Float_t  f1= 46.411;
1278     Float_t  f2= 228.71;
1279     for (i=0;i<26;i++)
1280     {
1281         Float_t ene=ppckov[i]*1e9;
1282         Float_t a=f1/(e1*e1 - ene*ene);
1283         Float_t b=f2/(e2*e2 - ene*ene);
1284         rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
1285         //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
1286     }
1287
1288     //Refraction index for opaque quarz, methane and grid
1289     Float_t rIndexOpaqueQuarz[26];
1290     Float_t rIndexMethane[26];
1291     Float_t rIndexGrid[26];
1292     for (i=0;i<26;i++)
1293     {
1294         rIndexOpaqueQuarz[i]=1;
1295         rIndexMethane[i]=1.000444;
1296         rIndexGrid[i]=1;
1297         //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
1298     }
1299
1300     //Absorption index for freon
1301     Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987,  179.0987, 179.0987, 179.0987,
1302                                179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
1303                                7.069031, 4.461292, 2.028366, 1.293013, .577267,   .40746,  .334964, 0., 0., 0.};
1304
1305     //Absorption index for quarz
1306     /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1307                         .906,.907,.907,.907};
1308     Float_t Wavl2[] = {150.,155.,160.0,165.0,170.0,175.0,180.0,185.0,190.0,195.0,200.0,205.0,210.0,
1309                        215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1310     Float_t abscoQuarz[31];
1311     for (Int_t i=0;i<31;i++)
1312     {
1313         Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1314         if (Xlam <= 160) abscoQuarz[i] = 0;
1315         if (Xlam > 250) abscoQuarz[i] = 1;
1316         else
1317         {
1318             for (Int_t j=0;j<21;j++)
1319             {
1320                 //printf ("Passed\n");
1321                 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1322                 {
1323                     Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1324                     Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1325                     abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1326                 }
1327             }
1328         }
1329         printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1330     }*/
1331
1332     /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1333                                39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1334                                17.94531, 11.88753, 5.99128,  3.83503,  2.36661,  1.53155, 1.30582, 1.08574, .8779708,
1335                                .675275, 0., 0., 0.};
1336
1337     for (Int_t i=0;i<31;i++)
1338     {
1339         abscoQuarz[i] = abscoQuarz[i]/10;
1340     }*/
1341
1342     Float_t abscoQuarz [26] = {105.8, 65.52, 48.58, 42.85, 35.79, 31.262, 28.598, 27.527, 25.007, 22.815, 21.004,
1343                                 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1344                                 .192, .1497, .10857};
1345
1346     //Absorption index for methane
1347     Float_t abscoMethane[26];
1348     for (i=0;i<26;i++)
1349     {
1350         abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1351         //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1352     }
1353
1354     //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1355     Float_t abscoOpaqueQuarz[26];
1356     Float_t abscoCsI[26];
1357     Float_t abscoGrid[26];
1358     Float_t efficAll[26];
1359     Float_t efficGrid[26];
1360     for (i=0;i<26;i++)
1361     {
1362         abscoOpaqueQuarz[i]=1e-5;
1363         abscoCsI[i]=1e-4;
1364         abscoGrid[i]=1e-4;
1365         efficAll[i]=1;
1366         efficGrid[i]=1;
1367         //printf ("All must be 1: %e,  %e,  %e,  %e,  %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1368     }
1369
1370     //Efficiency for csi
1371
1372     Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1373                              0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1374                              0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1375                              0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1376
1377
1378
1379     //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1380     //UNPOLARIZED PHOTONS
1381
1382     for (i=0;i<26;i++)
1383     {
1384         efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1385         //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1386     }
1387
1388     /*******************************************End of rich_media.f***************************************/
1389
1390
1391
1392
1393
1394
1395     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1396     zmet[2], zqua[2];
1397     Int_t nlmatfre;
1398     Float_t densquao;
1399     Int_t nlmatmet, nlmatqua;
1400     Float_t wmatquao[2], rIndexFreon[26];
1401     Float_t aquao[2], epsil, stmin, zquao[2];
1402     Int_t nlmatquao;
1403     Float_t radlal, densal, tmaxfd, deemax, stemax;
1404     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1405
1406     Int_t *idtmed = fIdtmed->GetArray()-999;
1407
1408     // --- Photon energy (GeV)
1409     // --- Refraction indexes
1410     for (i = 0; i < 26; ++i) {
1411       rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1412       //rIndexFreon[i] = 1;
1413         //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1414     }
1415
1416     // --- Detection efficiencies (quantum efficiency for CsI)
1417     // --- Define parameters for honeycomb.
1418     //     Used carbon of equivalent rad. lenght
1419
1420     ahon    = 12.01;
1421     zhon    = 6.;
1422     denshon = 0.1;
1423     radlhon = 18.8;
1424
1425     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1426
1427     aqua[0]    = 28.09;
1428     aqua[1]    = 16.;
1429     zqua[0]    = 14.;
1430     zqua[1]    = 8.;
1431     densqua    = 2.64;
1432     nlmatqua   = -2;
1433     wmatqua[0] = 1.;
1434     wmatqua[1] = 2.;
1435
1436     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1437
1438     aquao[0]    = 28.09;
1439     aquao[1]    = 16.;
1440     zquao[0]    = 14.;
1441     zquao[1]    = 8.;
1442     densquao    = 2.64;
1443     nlmatquao   = -2;
1444     wmatquao[0] = 1.;
1445     wmatquao[1] = 2.;
1446
1447     // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1448
1449     afre[0]    = 12.;
1450     afre[1]    = 19.;
1451     zfre[0]    = 6.;
1452     zfre[1]    = 9.;
1453     densfre    = 1.7;
1454     nlmatfre   = -2;
1455     wmatfre[0] = 6.;
1456     wmatfre[1] = 14.;
1457
1458     // --- Parameters to include in GSMIXT, relative to methane (CH4)
1459
1460     amet[0]    = 12.01;
1461     amet[1]    = 1.;
1462     zmet[0]    = 6.;
1463     zmet[1]    = 1.;
1464     densmet    = 7.17e-4;
1465     nlmatmet   = -2;
1466     wmatmet[0] = 1.;
1467     wmatmet[1] = 4.;
1468
1469     // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1470
1471     agri    = 63.54;
1472     zgri    = 29.;
1473     densgri = 8.96;
1474     radlgri = 1.43;
1475
1476     // --- Parameters to include in GSMATE related to aluminium sheet
1477
1478     aal    = 26.98;
1479     zal    = 13.;
1480     densal = 2.7;
1481     radlal = 8.9;
1482
1483     // --- Glass parameters
1484
1485     Float_t aglass[5]={12.01, 28.09, 16.,   10.8,  23.};
1486     Float_t zglass[5]={ 6.,   14.,    8.,    5.,   11.};
1487     Float_t wglass[5]={ 0.5,  0.105, 0.355, 0.03,  0.01};
1488     Float_t dglass=1.74;
1489
1490
1491     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
1492     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1493     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1494     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1495     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1496     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1497     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1498     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1499     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1500     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1501     AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1502     AliMaterial(31, "COPPER$",   63.54,    29.,   8.96,  1.4, 0.);
1503
1504     tmaxfd = -10.;
1505     stemax = -.1;
1506     deemax = -.2;
1507     epsil  = .001;
1508     stmin  = -.001;
1509
1510     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1511     AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1512     AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1513     AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1514     AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1515     AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1516     AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1517     AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1518     AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1519     AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1520     AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1521     AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1522
1523
1524     gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1525     gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1526     gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1527     gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1528     gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1529     gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1530     gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1531     gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1532     gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1533     gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1534     gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1535 }
1536
1537 //___________________________________________
1538
1539 Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1540 {
1541
1542     //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1543
1544     Float_t en[36] = {5.0,5.1,5.2,5.3,5.4,5.5,5.6,5.7,5.8,5.9,6.0,6.1,6.2,
1545                       6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,
1546                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1547
1548
1549     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1550                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1551                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1552                         1.72,1.53};
1553
1554     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1555                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1556                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1557                         1.714,1.498};
1558     Float_t xe=ene;
1559     Int_t  j=Int_t(xe*10)-49;
1560     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1561     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1562
1563     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1564     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1565
1566     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1567     Float_t tanin=sinin/pdoti;
1568
1569     Float_t c1=cn*cn-ck*ck-sinin*sinin;
1570     Float_t c2=4*cn*cn*ck*ck;
1571     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1572     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1573
1574     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1575     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1576
1577
1578     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1579     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
1580
1581     Float_t sigraf=18.;
1582     Float_t lamb=1240/ene;
1583     Float_t fresn;
1584
1585     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1586
1587     if(pola)
1588     {
1589         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
1590         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1591     }
1592     else
1593         fresn=0.5*(rp+rs);
1594
1595     fresn = fresn*rO;
1596     return(fresn);
1597 }
1598
1599 //__________________________________________
1600 Float_t AliRICH::AbsoCH4(Float_t x)
1601 {
1602
1603     //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1604     Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
1605     //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1606     Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1607     const Float_t kLosch=2.686763E19;                                      // LOSCHMIDT NUMBER IN CM-3
1608     const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1609     Float_t pn=kPressure/760.;
1610     Float_t tn=kTemperature/273.16;
1611
1612
1613 // ------- METHANE CROSS SECTION -----------------
1614 // ASTROPH. J. 214, L47 (1978)
1615
1616     Float_t sm=0;
1617     if (x<7.75)
1618         sm=.06e-22;
1619
1620     if(x>=7.75 && x<=8.1)
1621     {
1622         Float_t c0=-1.655279e-1;
1623         Float_t c1=6.307392e-2;
1624         Float_t c2=-8.011441e-3;
1625         Float_t c3=3.392126e-4;
1626         sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1627     }
1628
1629     if (x> 8.1)
1630     {
1631         Int_t j=0;
1632         while (x<=em[j] && x>=em[j+1])
1633         {
1634             j++;
1635             Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1636             sm=(sch4[j]+a*(x-em[j]))*1e-22;
1637         }
1638     }
1639
1640     Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1641     Float_t abslm=1./sm/dm;
1642
1643 //    ------- ISOBUTHANE CROSS SECTION --------------
1644 //     i-C4H10 (ai) abs. length from curves in
1645 //     Lu-McDonald paper for BARI RICH workshop .
1646 //     -----------------------------------------------------------
1647
1648     Float_t ai;
1649     Float_t absli;
1650     if (kIgas2 != 0)
1651     {
1652         if (x<7.25)
1653             ai=100000000.;
1654
1655         if(x>=7.25 && x<7.375)
1656             ai=24.3;
1657
1658         if(x>=7.375)
1659             ai=.0000000001;
1660
1661         Float_t si = 1./(ai*kLosch*273.16/293.);                    // ISOB. CRO.SEC.IN CM2
1662         Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1663         absli =1./si/di;
1664     }
1665     else
1666         absli=1.e18;
1667 //    ---------------------------------------------------------
1668 //
1669 //       transmission of O2
1670 //
1671 //       y= path in cm, x=energy in eV
1672 //       so= cross section for UV absorption in cm2
1673 //       do= O2 molecular density in cm-3
1674 //    ---------------------------------------------------------
1675
1676     Float_t abslo;
1677     Float_t so=0;
1678     if(x>=6.0)
1679     {
1680         if(x>=6.0 && x<6.5)
1681         {
1682             so=3.392709e-13 * TMath::Exp(2.864104 *x);
1683             so=so*1e-18;
1684         }
1685
1686         if(x>=6.5 && x<7.0)
1687         {
1688             so=2.910039e-34 * TMath::Exp(10.3337*x);
1689             so=so*1e-18;
1690         }
1691
1692
1693         if (x>=7.0)
1694         {
1695             Float_t a0=-73770.76;
1696             Float_t a1=46190.69;
1697             Float_t a2=-11475.44;
1698             Float_t a3=1412.611;
1699             Float_t a4=-86.07027;
1700             Float_t a5=2.074234;
1701             so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1702             so=so*1e-18;
1703         }
1704
1705         Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1706         abslo=1./so/dox;
1707     }
1708     else
1709         abslo=1.e18;
1710 //     ---------------------------------------------------------
1711 //
1712 //       transmission of H2O
1713 //
1714 //       y= path in cm, x=energy in eV
1715 //       sw= cross section for UV absorption in cm2
1716 //       dw= H2O molecular density in cm-3
1717 //     ---------------------------------------------------------
1718
1719     Float_t abslw;
1720
1721     Float_t b0=29231.65;
1722     Float_t b1=-15807.74;
1723     Float_t b2=3192.926;
1724     Float_t b3=-285.4809;
1725     Float_t b4=9.533944;
1726
1727     if(x>6.75)
1728     {
1729         Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1730         sw=sw*1e-18;
1731         Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1732         abslw=1./sw/dw;
1733     }
1734     else
1735         abslw=1.e18;
1736
1737 //    ---------------------------------------------------------
1738
1739     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1740     return (alength);
1741 }
1742
1743
1744
1745 //___________________________________________
1746 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1747 {
1748
1749 // Default value
1750
1751     return 9999;
1752 }
1753
1754 //___________________________________________
1755 void AliRICH::MakeBranch(Option_t* option, char *file)
1756 {
1757   // Create Tree branches for the RICH.
1758
1759     const Int_t kBufferSize = 4000;
1760     char branchname[20];
1761
1762     AliDetector::MakeBranch(option,file);
1763
1764     const char *cH = strstr(option,"H");
1765     const char *cD = strstr(option,"D");
1766     const char *cR = strstr(option,"R");
1767     const char *cS = strstr(option,"S");
1768
1769
1770     if (cH) {
1771       sprintf(branchname,"%sCerenkov",GetName());
1772       if (fCerenkovs   && gAlice->TreeH()) {
1773         //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1774         gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1775         //branch->SetAutoDelete(kFALSE);
1776       }
1777       sprintf(branchname,"%sSDigits",GetName());
1778       if (fSDigits   && gAlice->TreeH()) {
1779         //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1780         gAlice->MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1781         //branch->SetAutoDelete(kFALSE);
1782         //printf("Making branch %sSDigits in TreeH\n",GetName());
1783       }
1784     }
1785
1786     if (cS) {
1787       sprintf(branchname,"%sSDigits",GetName());
1788       if (fSDigits   && gAlice->TreeS()) {
1789         //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1790         gAlice->MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1791         //branch->SetAutoDelete(kFALSE);
1792         //printf("Making branch %sSDigits in TreeS\n",GetName());
1793       }
1794     }
1795
1796     if (cD) {
1797     //
1798     // one branch for digits per chamber
1799     //
1800       Int_t i;
1801
1802       for (i=0; i<kNCH ;i++) {
1803         sprintf(branchname,"%sDigits%d",GetName(),i+1);
1804         if (fDchambers   && gAlice->TreeD()) {
1805            //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1806            gAlice->MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1807            //branch->SetAutoDelete(kFALSE);
1808           //printf("Making Branch %sDigits%d\n",GetName(),i+1);
1809         }
1810       }
1811     }
1812
1813     if (cR) {
1814     //
1815     // one branch for raw clusters per chamber
1816     //
1817
1818       //printf("Called MakeBranch for TreeR\n");
1819
1820       Int_t i;
1821
1822       for (i=0; i<kNCH ;i++) {
1823         sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1824         if (fRawClusters && gAlice->TreeR()) {
1825            //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1826           gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1827            //branch->SetAutoDelete(kFALSE);
1828         }
1829       }
1830      //
1831      // one branch for rec hits per chamber
1832      //
1833      for (i=0; i<kNCH ;i++) {
1834        sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1835        if (fRecHits1D   && gAlice->TreeR()) {
1836          //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1837          gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1838          //branch->SetAutoDelete(kFALSE);
1839        }
1840      }
1841      for (i=0; i<kNCH ;i++) {
1842        sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1843        if (fRecHits3D   && gAlice->TreeR()) {
1844          //TBranch* branch = gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
1845          gAlice->MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1846          //branch->SetAutoDelete(kFALSE);
1847       }
1848     }
1849   }
1850 }
1851
1852 //___________________________________________
1853 void AliRICH::SetTreeAddress()
1854 {
1855   // Set branch address for the Hits and Digits Tree.
1856   char branchname[20];
1857   Int_t i;
1858
1859     AliDetector::SetTreeAddress();
1860
1861     TBranch *branch;
1862     TTree *treeH = gAlice->TreeH();
1863     TTree *treeD = gAlice->TreeD();
1864     TTree *treeR = gAlice->TreeR();
1865     TTree *treeS = gAlice->TreeS();
1866
1867     if (treeH) {
1868       if (fCerenkovs) {
1869             branch = treeH->GetBranch("RICHCerenkov");
1870             if (branch) branch->SetAddress(&fCerenkovs);
1871         }
1872     if (fSDigits) {
1873         branch = treeH->GetBranch("RICHSDigits");
1874         if (branch)
1875           {
1876             branch->SetAddress(&fSDigits);
1877             //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits);
1878           }
1879       }
1880     }
1881
1882     if (treeS) {
1883       if (fSDigits) {
1884         branch = treeS->GetBranch("RICHSDigits");
1885         if (branch)
1886           {
1887             branch->SetAddress(&fSDigits);
1888             //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits);
1889           }
1890       }
1891     }
1892
1893
1894     if (treeD) {
1895         for (int i=0; i<kNCH; i++) {
1896             sprintf(branchname,"%sDigits%d",GetName(),i+1);
1897             if (fDchambers) {
1898               branch = treeD->GetBranch(branchname);
1899               if (branch) branch->SetAddress(&((*fDchambers)[i]));
1900             }
1901         }
1902     }
1903   if (treeR) {
1904       for (i=0; i<kNCH; i++) {
1905           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1906           if (fRawClusters) {
1907               branch = treeR->GetBranch(branchname);
1908               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
1909           }
1910       }
1911
1912       for (i=0; i<kNCH; i++) {
1913         sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1914         if (fRecHits1D) {
1915           branch = treeR->GetBranch(branchname);
1916           if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
1917           }
1918       }
1919
1920      for (i=0; i<kNCH; i++) {
1921         sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1922         if (fRecHits3D) {
1923           branch = treeR->GetBranch(branchname);
1924           if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
1925           }
1926       }
1927
1928   }
1929 }
1930 //___________________________________________
1931 void AliRICH::ResetHits()
1932 {
1933   // Reset number of clusters and the cluster array for this detector
1934     AliDetector::ResetHits();
1935     fNSDigits   = 0;
1936     fNcerenkovs = 0;
1937     if (fSDigits)  fSDigits->Clear();
1938     if (fCerenkovs) fCerenkovs->Clear();
1939 }
1940
1941
1942 //____________________________________________
1943 void AliRICH::ResetDigits()
1944 {
1945   //
1946   // Reset number of digits and the digits array for this detector
1947   //
1948     for ( int i=0;i<kNCH;i++ ) {
1949         if (fDchambers && (*fDchambers)[i])   (*fDchambers)[i]->Clear();
1950         if (fNdch)  fNdch[i]=0;
1951     }
1952 }
1953
1954 //____________________________________________
1955 void AliRICH::ResetRawClusters()
1956 {
1957   //
1958   // Reset number of raw clusters and the raw clust array for this detector
1959   //
1960     for ( int i=0;i<kNCH;i++ ) {
1961         if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
1962         if (fNrawch)  fNrawch[i]=0;
1963     }
1964 }
1965
1966 //____________________________________________
1967 void AliRICH::ResetRecHits1D()
1968 {
1969   //
1970   // Reset number of raw clusters and the raw clust array for this detector
1971   //
1972
1973   for ( int i=0;i<kNCH;i++ ) {
1974         if ((*fRecHits1D)[i])    ((TClonesArray*)(*fRecHits1D)[i])->Clear();
1975         if (fNrechits1D)  fNrechits1D[i]=0;
1976     }
1977 }
1978
1979 //____________________________________________
1980 void AliRICH::ResetRecHits3D()
1981 {
1982   //
1983   // Reset number of raw clusters and the raw clust array for this detector
1984   //
1985
1986   for ( int i=0;i<kNCH;i++ ) {
1987         if ((*fRecHits3D)[i])    ((TClonesArray*)(*fRecHits3D)[i])->Clear();
1988         if (fNrechits3D)  fNrechits3D[i]=0;
1989     }
1990 }
1991
1992 //___________________________________________
1993 void   AliRICH::SetGeometryModel(Int_t id, AliRICHGeometry *geometry)
1994 {
1995
1996 //
1997 // Setter for the RICH geometry model
1998 //
1999
2000
2001     ((AliRICHChamber*) (*fChambers)[id])->GeometryModel(geometry);
2002 }
2003
2004 //___________________________________________
2005 void   AliRICH::SetSegmentationModel(Int_t id, AliSegmentation *segmentation)
2006 {
2007
2008 //
2009 // Setter for the RICH segmentation model
2010 //
2011
2012     ((AliRICHChamber*) (*fChambers)[id])->SetSegmentationModel(segmentation);
2013 }
2014
2015 //___________________________________________
2016 void   AliRICH::SetResponseModel(Int_t id, AliRICHResponse *response)
2017 {
2018
2019 //
2020 // Setter for the RICH response model
2021 //
2022
2023     ((AliRICHChamber*) (*fChambers)[id])->ResponseModel(response);
2024 }
2025
2026 void   AliRICH::SetReconstructionModel(Int_t id, AliRICHClusterFinder *reconst)
2027 {
2028
2029 //
2030 // Setter for the RICH reconstruction model (clusters)
2031 //
2032
2033     ((AliRICHChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
2034 }
2035
2036 //___________________________________________
2037 void AliRICH::StepManager()
2038 {
2039
2040 // Full Step Manager
2041
2042     Int_t          copy, id;
2043     static Int_t   idvol;
2044     static Int_t   vol[2];
2045     Int_t          ipart;
2046     static Float_t hits[22];
2047     static Float_t ckovData[19];
2048     TLorentzVector position;
2049     TLorentzVector momentum;
2050     Float_t        pos[3];
2051     Float_t        mom[4];
2052     Float_t        localPos[3];
2053     Float_t        localMom[4];
2054     Float_t        localTheta,localPhi;
2055     Float_t        theta,phi;
2056     Float_t        destep, step;
2057     Float_t        ranf[2];
2058     Int_t          nPads;
2059     Float_t        coscerenkov;
2060     static Float_t eloss, xhit, yhit, tlength;
2061     const  Float_t kBig=1.e10;
2062
2063     TClonesArray &lhits = *fHits;
2064     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
2065
2066  //if (current->Energy()>1)
2067    //{
2068
2069     // Only gas gap inside chamber
2070     // Tag chambers and record hits when track enters
2071
2072     idvol=-1;
2073     id=gMC->CurrentVolID(copy);
2074     Float_t cherenkovLoss=0;
2075     //gAlice->KeepTrack(gAlice->CurrentTrack());
2076
2077     gMC->TrackPosition(position);
2078     pos[0]=position(0);
2079     pos[1]=position(1);
2080     pos[2]=position(2);
2081     //bzero((char *)ckovData,sizeof(ckovData)*19);
2082     ckovData[1] = pos[0];                 // X-position for hit
2083     ckovData[2] = pos[1];                 // Y-position for hit
2084     ckovData[3] = pos[2];                 // Z-position for hit
2085     ckovData[6] = 0;                      // dummy track length
2086     //ckovData[11] = gAlice->CurrentTrack();
2087
2088     //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
2089
2090     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
2091
2092     /********************Store production parameters for Cerenkov photons************************/
2093 //is it a Cerenkov photon?
2094     if (gMC->TrackPid() == 50000050) {
2095
2096       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
2097         //{
2098           Float_t ckovEnergy = current->Energy();
2099           //energy interval for tracking
2100           if  (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
2101             //if (ckovEnergy > 0)
2102             {
2103               if (gMC->IsTrackEntering()){        //is track entering?
2104                 //printf("Track entered (1)\n");
2105                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2106                   {                                                          //is it in freo?
2107                     if (gMC->IsNewTrack()){                          //is it the first step?
2108                       //printf("I'm in!\n");
2109                       Int_t mother = current->GetFirstMother();
2110
2111                       //printf("Second Mother:%d\n",current->GetSecondMother());
2112
2113                       ckovData[10] = mother;
2114                       ckovData[11] = gAlice->CurrentTrack();
2115                       ckovData[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
2116                       //printf("Produced in FREO\n");
2117                       fCkovNumber++;
2118                       fFreonProd=1;
2119                       //printf("Index: %d\n",fCkovNumber);
2120                     }    //first step question
2121                   }        //freo question
2122
2123                 if (gMC->IsNewTrack()){                                  //is it first step?
2124                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
2125                     {
2126                       ckovData[12] = 2;
2127                       //printf("Produced in QUAR\n");
2128                     }    //quarz question
2129                 }        //first step question
2130
2131                 //printf("Before %d\n",fFreonProd);
2132               }   //track entering question
2133
2134               if (ckovData[12] == 1)                                        //was it produced in Freon?
2135                 //if (fFreonProd == 1)
2136                 {
2137                   if (gMC->IsTrackEntering()){                                     //is track entering?
2138                     //printf("Track entered (2)\n");
2139                     //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
2140                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
2141                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
2142                       {
2143                         //printf("Got in META\n");
2144                         gMC->TrackMomentum(momentum);
2145                         mom[0]=momentum(0);
2146                         mom[1]=momentum(1);
2147                         mom[2]=momentum(2);
2148                         mom[3]=momentum(3);
2149                         // Z-position for hit
2150
2151
2152                         /**************** Photons lost in second grid have to be calculated by hand************/
2153
2154                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2155                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
2156                         gMC->Rndm(ranf, 1);
2157                         //printf("grid calculation:%f\n",t);
2158                         if (ranf[0] > t) {
2159                           gMC->StopTrack();
2160                           ckovData[13] = 5;
2161                           AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2162                           //printf("Added One (1)!\n");
2163                           //printf("Lost one in grid\n");
2164                         }
2165                         /**********************************************************************************/
2166                       }    //gap
2167
2168                     //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
2169                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2170                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
2171                       {
2172                         //printf("Got in CSI\n");
2173                         gMC->TrackMomentum(momentum);
2174                         mom[0]=momentum(0);
2175                         mom[1]=momentum(1);
2176                         mom[2]=momentum(2);
2177                         mom[3]=momentum(3);
2178
2179                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
2180                         /***********************Cerenkov phtons (always polarised)*************************/
2181
2182                         Float_t cophi = TMath::Cos(TMath::ATan2(mom[0], mom[1]));
2183                         Float_t t = Fresnel(ckovEnergy*1e9,cophi,1);
2184                         gMC->Rndm(ranf, 1);
2185                         if (ranf[0] < t) {
2186                           gMC->StopTrack();
2187                           ckovData[13] = 6;
2188                           AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2189                           //printf("Added One (2)!\n");
2190                           //printf("Lost by Fresnel\n");
2191                         }
2192                         /**********************************************************************************/
2193                       }
2194                   } //track entering?
2195
2196
2197                   /********************Evaluation of losses************************/
2198                   /******************still in the old fashion**********************/
2199
2200                   TArrayI procs;
2201                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
2202                   for (Int_t i = 0; i < i1; ++i) {
2203                     //        Reflection loss
2204                     if (procs[i] == kPLightReflection) {        //was it reflected
2205                       ckovData[13]=10;
2206                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2207                         ckovData[13]=1;
2208                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2209                         ckovData[13]=2;
2210                       //gMC->StopTrack();
2211                       //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2212                     } //reflection question
2213
2214                     //        Absorption loss
2215                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
2216                       //printf("Got in absorption\n");
2217                       ckovData[13]=20;
2218                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2219                         ckovData[13]=11;
2220                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2221                         ckovData[13]=12;
2222                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
2223                         ckovData[13]=13;
2224                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
2225                         ckovData[13]=13;
2226
2227                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
2228                         ckovData[13]=15;
2229
2230                       //        CsI inefficiency
2231                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
2232                         ckovData[13]=16;
2233                       }
2234                       gMC->StopTrack();
2235                       AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2236                       //printf("Added One (3)!\n");
2237                       //printf("Added cerenkov %d\n",fCkovNumber);
2238                     } //absorption question
2239
2240
2241                     //        Photon goes out of tracking scope
2242                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
2243                       ckovData[13]=21;
2244                       gMC->StopTrack();
2245                       AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2246                       //printf("Added One (4)!\n");
2247                     }   // energy treshold question
2248                   }  //number of mechanisms cycle
2249                   /**********************End of evaluation************************/
2250                 } //freon production question
2251             } //energy interval question
2252         //}//inside the proximity gap question
2253     } //cerenkov photon question
2254
2255     /**************************************End of Production Parameters Storing*********************/
2256
2257
2258     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
2259
2260     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
2261       //printf("Cerenkov\n");
2262         if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
2263         {
2264           //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
2265           //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2266           //printf("Got in CSI\n");
2267           if (gMC->Edep() > 0.){
2268                 gMC->TrackPosition(position);
2269                 gMC->TrackMomentum(momentum);
2270                 pos[0]=position(0);
2271                 pos[1]=position(1);
2272                 pos[2]=position(2);
2273                 mom[0]=momentum(0);
2274                 mom[1]=momentum(1);
2275                 mom[2]=momentum(2);
2276                 mom[3]=momentum(3);
2277                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2278                 Double_t rt = TMath::Sqrt(tc);
2279                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2280                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2281                 gMC->Gmtod(pos,localPos,1);
2282                 gMC->Gmtod(mom,localMom,2);
2283
2284                 gMC->CurrentVolOffID(2,copy);
2285                 vol[0]=copy;
2286                 idvol=vol[0]-1;
2287
2288                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2289                         //->Sector(localPos[0], localPos[2]);
2290                 //printf("Sector:%d\n",sector);
2291
2292                 /*if (gMC->TrackPid() == 50000051){
2293                   fFeedbacks++;
2294                   printf("Feedbacks:%d\n",fFeedbacks);
2295                 }*/
2296
2297                 ((AliRICHChamber*) (*fChambers)[idvol])
2298                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2299                 if(idvol<kNCH) {
2300                     ckovData[0] = gMC->TrackPid();        // particle type
2301                     ckovData[1] = pos[0];                 // X-position for hit
2302                     ckovData[2] = pos[1];                 // Y-position for hit
2303                     ckovData[3] = pos[2];                 // Z-position for hit
2304                     ckovData[4] = theta;                      // theta angle of incidence
2305                     ckovData[5] = phi;                      // phi angle of incidence
2306                     ckovData[8] = (Float_t) fNSDigits;      // first sdigit
2307                     ckovData[9] = -1;                       // last pad hit
2308                     ckovData[13] = 4;                       // photon was detected
2309                     ckovData[14] = mom[0];
2310                     ckovData[15] = mom[1];
2311                     ckovData[16] = mom[2];
2312
2313                     destep = gMC->Edep();
2314                     gMC->SetMaxStep(kBig);
2315                     cherenkovLoss  += destep;
2316                     ckovData[7]=cherenkovLoss;
2317
2318                     nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
2319
2320                     if (fNSDigits > (Int_t)ckovData[8]) {
2321                         ckovData[8]= ckovData[8]+1;
2322                         ckovData[9]= (Float_t) fNSDigits;
2323                     }
2324
2325                     //printf("Cerenkov loss: %f\n", cherenkovLoss);
2326
2327                     ckovData[17] = nPads;
2328                     //printf("nPads:%d",nPads);
2329
2330                     //TClonesArray *Hits = RICH->Hits();
2331                     AliRICHHit *mipHit =  (AliRICHHit*) (fHits->UncheckedAt(0));
2332                     if (mipHit)
2333                       {
2334                         mom[0] = current->Px();
2335                         mom[1] = current->Py();
2336                         mom[2] = current->Pz();
2337                         Float_t mipPx = mipHit->fMomX;
2338                         Float_t mipPy = mipHit->fMomY;
2339                         Float_t mipPz = mipHit->fMomZ;
2340
2341                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
2342                         Float_t rt = TMath::Sqrt(r);
2343                         Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
2344                         Float_t mipRt = TMath::Sqrt(mipR);
2345                         if ((rt*mipRt) > 0)
2346                           {
2347                             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
2348                           }
2349                         else
2350                           {
2351                             coscerenkov = 0;
2352                           }
2353                         Float_t cherenkov = TMath::ACos(coscerenkov);
2354                         ckovData[18]=cherenkov;
2355                       }
2356                     //if (sector != -1)
2357                     //{
2358                     AddHit(gAlice->CurrentTrack(),vol,ckovData);
2359                     AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2360                     //printf("Added One (5)!\n");
2361                     //}
2362                 }
2363             }
2364         }
2365     }
2366
2367     /***********************************************End of photon hits*********************************************/
2368
2369
2370     /**********************************************Charged particles treatment*************************************/
2371
2372     else if (gMC->TrackCharge())
2373     //else if (1 == 1)
2374       {
2375 //If MIP
2376         /*if (gMC->IsTrackEntering())
2377           {
2378             hits[13]=20;//is track entering?
2379           }*/
2380         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2381           {
2382             gMC->TrackMomentum(momentum);
2383             mom[0]=momentum(0);
2384             mom[1]=momentum(1);
2385             mom[2]=momentum(2);
2386             mom[3]=momentum(3);
2387             hits [19] = mom[0];
2388             hits [20] = mom[1];
2389             hits [21] = mom[2];
2390             fFreonProd=1;
2391           }
2392
2393         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2394 // Get current particle id (ipart), track position (pos)  and momentum (mom)
2395
2396             gMC->CurrentVolOffID(3,copy);
2397             vol[0]=copy;
2398             idvol=vol[0]-1;
2399
2400             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2401                         //->Sector(localPos[0], localPos[2]);
2402             //printf("Sector:%d\n",sector);
2403
2404             gMC->TrackPosition(position);
2405             gMC->TrackMomentum(momentum);
2406             pos[0]=position(0);
2407             pos[1]=position(1);
2408             pos[2]=position(2);
2409             mom[0]=momentum(0);
2410             mom[1]=momentum(1);
2411             mom[2]=momentum(2);
2412             mom[3]=momentum(3);
2413             gMC->Gmtod(pos,localPos,1);
2414             gMC->Gmtod(mom,localMom,2);
2415
2416             ipart  = gMC->TrackPid();
2417             //
2418             // momentum loss and steplength in last step
2419             destep = gMC->Edep();
2420             step   = gMC->TrackStep();
2421
2422             //
2423             // record hits when track enters ...
2424             if( gMC->IsTrackEntering()) {
2425 //              gMC->SetMaxStep(fMaxStepGas);
2426                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2427                 Double_t rt = TMath::Sqrt(tc);
2428                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2429                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2430
2431
2432                 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2433                 Double_t localRt = TMath::Sqrt(localTc);
2434                 localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2435                 localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2436
2437                 hits[0] = Float_t(ipart);         // particle type
2438                 hits[1] = localPos[0];                 // X-position for hit
2439                 hits[2] = localPos[1];                 // Y-position for hit
2440                 hits[3] = localPos[2];                 // Z-position for hit
2441                 hits[4] = localTheta;                  // theta angle of incidence
2442                 hits[5] = localPhi;                    // phi angle of incidence
2443                 hits[8] = (Float_t) fNSDigits;    // first sdigit
2444                 hits[9] = -1;                     // last pad hit
2445                 hits[13] = fFreonProd;           // did id hit the freon?
2446                 hits[14] = mom[0];
2447                 hits[15] = mom[1];
2448                 hits[16] = mom[2];
2449                 hits[18] = 0;               // dummy cerenkov angle
2450
2451                 tlength = 0;
2452                 eloss   = 0;
2453                 fFreonProd = 0;
2454
2455                 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2456
2457
2458                 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2459                 xhit    = localPos[0];
2460                 yhit    = localPos[2];
2461                 // Only if not trigger chamber
2462                 if(idvol<kNCH) {
2463                     //
2464                     //  Initialize hit position (cursor) in the segmentation model
2465                     ((AliRICHChamber*) (*fChambers)[idvol])
2466                         ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2467                 }
2468             }
2469
2470             //
2471             // Calculate the charge induced on a pad (disintegration) in case
2472             //
2473             // Mip left chamber ...
2474             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2475                 gMC->SetMaxStep(kBig);
2476                 eloss   += destep;
2477                 tlength += step;
2478
2479
2480                 // Only if not trigger chamber
2481                 if(idvol<kNCH) {
2482                   if (eloss > 0)
2483                     {
2484                       if(gMC->TrackPid() == kNeutron)
2485                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2486                       nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2487                       hits[17] = nPads;
2488                       //printf("nPads:%d",nPads);
2489                     }
2490                 }
2491
2492                 hits[6]=tlength;
2493                 hits[7]=eloss;
2494                 if (fNSDigits > (Int_t)hits[8]) {
2495                     hits[8]= hits[8]+1;
2496                     hits[9]= (Float_t) fNSDigits;
2497                 }
2498
2499                 //if(sector !=-1)
2500                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2501                 eloss = 0;
2502                 //
2503                 // Check additional signal generation conditions
2504                 // defined by the segmentation
2505                 // model (boundary crossing conditions)
2506             } else if
2507                 (((AliRICHChamber*) (*fChambers)[idvol])
2508                  ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2509             {
2510                 ((AliRICHChamber*) (*fChambers)[idvol])
2511                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2512                 if (eloss > 0)
2513                   {
2514                     if(gMC->TrackPid() == kNeutron)
2515                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2516                     nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2517                     hits[17] = nPads;
2518                     //printf("Npads:%d",NPads);
2519                   }
2520                 xhit     = localPos[0];
2521                 yhit     = localPos[2];
2522                 eloss    = destep;
2523                 tlength += step ;
2524                 //
2525                 // nothing special  happened, add up energy loss
2526             } else {
2527                 eloss   += destep;
2528                 tlength += step ;
2529             }
2530         }
2531       }
2532     /*************************************************End of MIP treatment**************************************/
2533    //}
2534 }
2535
2536 void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
2537 {
2538
2539 //
2540 // Loop on chambers and on cathode planes
2541 //
2542     for (Int_t icat=1;icat<2;icat++) {
2543         gAlice->ResetDigits();
2544         gAlice->TreeD()->GetEvent(0);
2545         for (Int_t ich=0;ich<kNCH;ich++) {
2546           AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
2547           TClonesArray *pRICHdigits  = this->DigitsAddress(ich);
2548           if (pRICHdigits == 0)
2549               continue;
2550           //
2551           // Get ready the current chamber stuff
2552           //
2553           AliRICHResponse* response = iChamber->GetResponseModel();
2554           AliSegmentation*  seg = iChamber->GetSegmentationModel();
2555           AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2556           if (seg) {
2557               rec->SetSegmentation(seg);
2558               rec->SetResponse(response);
2559               rec->SetDigits(pRICHdigits);
2560               rec->SetChamber(ich);
2561               if (nev==0) rec->CalibrateCOG();
2562               rec->FindRawClusters();
2563           }
2564           TClonesArray *fRch;
2565           fRch=RawClustAddress(ich);
2566           fRch->Sort();
2567         } // for ich
2568
2569         gAlice->TreeR()->Fill();
2570         TClonesArray *fRch;
2571         for (int i=0;i<kNCH;i++) {
2572             fRch=RawClustAddress(i);
2573             int nraw=fRch->GetEntriesFast();
2574             printf ("Chamber %d, raw clusters %d\n",i,nraw);
2575         }
2576
2577         ResetRawClusters();
2578
2579     } // for icat
2580
2581     char hname[30];
2582     sprintf(hname,"TreeR%d",nev);
2583     gAlice->TreeR()->Write(hname,kOverwrite,0);
2584     gAlice->TreeR()->Reset();
2585
2586     //gObjectTable->Print();
2587 }
2588
2589 AliRICHSDigit* AliRICH::FirstPad(AliRICHHit*  hit,TClonesArray *clusters )
2590 {
2591 //
2592     // Initialise the pad iterator
2593     // Return the address of the first sdigit for hit
2594     TClonesArray *theClusters = clusters;
2595     Int_t nclust = theClusters->GetEntriesFast();
2596     if (nclust && hit->fPHlast > 0) {
2597         sMaxIterPad=Int_t(hit->fPHlast);
2598         sCurIterPad=Int_t(hit->fPHfirst);
2599         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2600     } else {
2601         return 0;
2602     }
2603
2604 }
2605
2606 AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters)
2607 {
2608
2609   // Iterates over pads
2610
2611     sCurIterPad++;
2612     if (sCurIterPad <= sMaxIterPad) {
2613         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2614     } else {
2615         return 0;
2616     }
2617 }
2618
2619 AliRICH& AliRICH::operator=(const AliRICH& rhs)
2620 {
2621 // Assignment operator
2622     return *this;
2623
2624 }
2625
2626
2627