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