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