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