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