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