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