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