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