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