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