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