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