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.62  2002/10/31 08:44:04  morsch
  19   Problems with rotated RICH solved:
  20   Detector response (fresnel reflection, grid absorption ...) has to be
  21   determined using local coordinates.
  22
  23   Revision 1.61  2002/10/29 15:00:08  morsch
  24   - Diagnostics updated.
  25   - RecHits structure synchronized.
  26   - Digitizer method using AliRICHDigitizer.
  27   (J. Barbosa)
  28
  29
  30   Revision 1.60  2002/10/22 16:28:21  alibrary
  31   Introducing Riostream.h
  32
  33   Revision 1.59  2002/10/14 14:57:31  hristov
  34   Merging the VirtualMC branch to the main development branch (HEAD)
  35
  36   Revision 1.58.6.1  2002/06/10 15:12:46  hristov
  37   Merged with v3-08-02
  38
  39   Revision 1.58  2001/11/14 09:49:37  dibari
  40   Use debug methods
  41
  42   Revision 1.57  2001/11/09 17:29:31  dibari
  43   Setters fro models moved to header
  44
  45   Revision 1.56  2001/11/02 15:37:25  hristov
  46   Digitizer class created. Code cleaning and bug fixes (J.Chudoba)
  47
  48   Revision 1.55  2001/10/23 13:03:35  hristov
  49   The access to several data members was changed from public to protected. The digitisation was adapted to the multi-event case (J.Chudoba)
  50
  51   Revision 1.54  2001/09/07 08:38:10  hristov
  52   Pointers initialised to 0 in the default constructors
  53
  54   Revision 1.53  2001/08/30 09:51:23  hristov
  55   The operator[] is replaced by At() or AddAt() in case of TObjArray.
  56
  57   Revision 1.52  2001/05/16 14:57:20  alibrary
  58   New files for folders and Stack
  59
  60   Revision 1.51  2001/05/14 10:18:55  hristov
  61   Default arguments declared once
  62
  63   Revision 1.50  2001/05/10 14:44:16  jbarbosa
  64   Corrected some overlaps (thanks I. Hrivnacovna).
  65
  66   Revision 1.49  2001/05/10 12:23:49  jbarbosa
  67   Repositioned the RICH modules.
  68   Eliminated magic numbers.
  69   Incorporated diagnostics (from macros).
  70
  71   Revision 1.48  2001/03/15 10:35:00  jbarbosa
  72   Corrected bug in MakeBranch (was using a different version of STEER)
  73
  74   Revision 1.47  2001/03/14 18:13:56  jbarbosa
  75   Several changes to adapt to new IO.
  76   Removed digitising function, using AliRICHMerger::Digitise from now on.
  77
  78   Revision 1.46  2001/03/12 17:46:33  hristov
  79   Changes needed on Sun with CC 5.0
  80
  81   Revision 1.45  2001/02/27 22:11:46  jbarbosa
  82   Testing TreeS, removing of output.
  83
  84   Revision 1.44  2001/02/27 15:19:12  jbarbosa
  85   Transition to SDigits.
  86
  87   Revision 1.43  2001/02/23 17:19:06  jbarbosa
  88   Corrected photocathode definition in BuildGeometry().
  89
  90   Revision 1.42  2001/02/13 20:07:23  jbarbosa
  91   Parametrised definition of photcathode dimensions. New spacers. New data members in AliRICHHit to store particle momentum
  92   when entering the freon. Corrected calls to particle stack.
  93
  94   Revision 1.41  2001/01/26 20:00:20  hristov
  95   Major upgrade of AliRoot code
  96
  97   Revision 1.40  2001/01/24 20:58:03  jbarbosa
  98   Enhanced BuildGeometry. Now the photocathodes are drawn.
  99
 100   Revision 1.39  2001/01/22 21:40:24  jbarbosa
 101   Removing magic numbers
 102
 103   Revision 1.37  2000/12/20 14:07:25  jbarbosa
 104   Removed dependencies on TGeant3 (thanks to F. Carminati and I. Hrivnacova)
 105
 106   Revision 1.36  2000/12/18 17:45:54  jbarbosa
 107   Cleaned up PadHits object.
 108
 109   Revision 1.35  2000/12/15 16:49:40  jbarbosa
 110   Geometry and materials updates (wire supports, pcbs, backplane supports, frame).
 111
 112   Revision 1.34  2000/11/10 18:12:12  jbarbosa
 113   Bug fix for AliRICHCerenkov (thanks to P. Hristov)
 114
 115   Revision 1.33  2000/11/02 10:09:01  jbarbosa
 116   Minor bug correction (some pointers were not initialised in the default constructor)
 117
 118   Revision 1.32  2000/11/01 15:32:55  jbarbosa
 119   Updated to handle both reconstruction algorithms.
 120
 121   Revision 1.31  2000/10/26 20:18:33  jbarbosa
 122   Supports for methane and freon vessels
 123
 124   Revision 1.30  2000/10/24 13:19:12  jbarbosa
 125   Geometry updates.
 126
 127   Revision 1.29  2000/10/19 19:39:25  jbarbosa
 128   Some more changes to geometry. Further correction of digitisation "per part. type"
 129
 130   Revision 1.28  2000/10/17 20:50:57  jbarbosa
 131   Inversed digtise by particle type (now, only the selected particle type is not digitsed).
 132   Corrected several geometry minor bugs.
 133   Added new parameter (opaque quartz thickness).
 134
 135   Revision 1.27  2000/10/11 10:33:55  jbarbosa
 136   Corrected bug introduced by earlier revisions  (CerenkovData array cannot be reset to zero on wach call of StepManager)
 137
 138   Revision 1.26  2000/10/03 21:44:08  morsch
 139   Use AliSegmentation and AliHit abstract base classes.
 140
 141   Revision 1.25  2000/10/02 21:28:12  fca
 142   Removal of useless dependecies via forward declarations
 143
 144   Revision 1.24  2000/10/02 15:43:17  jbarbosa
 145   Fixed forward declarations.
 146   Fixed honeycomb density.
 147   Fixed cerenkov storing.
 148   New electronics.
 149
 150   Revision 1.23  2000/09/13 10:42:14  hristov
 151   Minor corrections for HP, DEC and Sun; strings.h included
 152
 153   Revision 1.22  2000/09/12 18:11:13  fca
 154   zero hits area before using
 155
 156   Revision 1.21  2000/07/21 10:21:07  morsch
 157   fNrawch   = 0; and  fNrechits = 0; in the default constructor.
 158
 159   Revision 1.20  2000/07/10 15:28:39  fca
 160   Correction of the inheritance scheme
 161
 162   Revision 1.19  2000/06/30 16:29:51  dibari
 163   Added kDebugLevel variable to control output size on demand
 164
 165   Revision 1.18  2000/06/12 15:15:46  jbarbosa
 166   Cleaned up version.
 167
 168   Revision 1.17  2000/06/09 14:58:37  jbarbosa
 169   New digitisation per particle type
 170
 171   Revision 1.16  2000/04/19 12:55:43  morsch
 172   Newly structured and updated version (JB, AM)
 173
 174 */
 175
 176
 177 ////////////////////////////////////////////////
 178 //  Manager and hits classes for set:RICH     //
 179 ////////////////////////////////////////////////
 180
 181 #include <TBRIK.h>
 182 #include <TTUBE.h>
 183 #include <TNode.h>
 184 #include <TRandom.h>
 185 #include <TObject.h>
 186 #include <TVector.h>
 187 #include <TObjArray.h>
 188 #include <TArrayF.h>
 189 #include <TFile.h>
 190 #include <TParticle.h>
 191 #include <TGeometry.h>
 192 #include <TTree.h>
 193 #include <TH1.h>
 194 #include <TH2.h>
 195 #include <TCanvas.h>
 196 #include <TStyle.h>
 197 #include <TF1.h>
 198
 199 #include <Riostream.h>
 200 #include <strings.h>
 201
 202 #include "AliRICH.h"
 203 #include "AliSegmentation.h"
 204 #include "AliRICHSegmentationV0.h"
 205 #include "AliRICHHit.h"
 206 #include "AliRICHCerenkov.h"
 207 #include "AliRICHSDigit.h"
 208 #include "AliRICHDigit.h"
 209 #include "AliRICHTransientDigit.h"
 210 #include "AliRICHRawCluster.h"
 211 #include "AliRICHRecHit1D.h"
 212 #include "AliRICHRecHit3D.h"
 213 #include "AliRICHHitMapA1.h"
 214 #include "AliRICHClusterFinder.h"
 215 #include "AliRICHMerger.h"
 216 #include "AliRICHDigitizer.h"
 217 #include "AliRun.h"
 218 #include "AliRunDigitizer.h"
 219 #include "AliMC.h"
 220 #include "AliMagF.h"
 221 #include "AliConst.h"
 222 #include "AliPDG.h"
 223 #include "AliPoints.h"
 224
 225
 226
 227
 228 static Int_t sMaxIterPad=0;    // Static variables for the pad-hit iterator routines
 229 static Int_t sCurIterPad=0;
 230
 231 ClassImp(AliRICH)
 232
 233 //___________________________________________
 234 // RICH manager class
 235 //Begin_Html
 236 /*
 237   <img src="gif/alirich.gif">
 238 */
 239 //End_Html
 240
 241 AliRICH::AliRICH()
 242 {
 243 // Default ctor should not contain any new operators
 244    cout<<ClassName()<<"::named ctor(sName,sTitle)>\n"; // no way to control it as ctor is called before call to SetDebugXXXX()
 245
 246     fIshunt     = 0;
 247     fHits       = 0;
 248     fSDigits    = 0;
 249     fNSDigits   = 0;
 250     fNcerenkovs = 0;
 251     fDchambers  = 0;
 252     fRecHits1D = 0;
 253     fRecHits3D = 0;
 254     fRawClusters = 0;
 255     fChambers = 0;
 256     fCerenkovs  = 0;
 257     for (Int_t i=0; i<7; i++){
 258         fNdch[i]       = 0;
 259         fNrawch[i]     = 0;
 260         fNrechits1D[i] = 0;
 261         fNrechits3D[i] = 0;
 262     }
 263
 264     fFileName = 0;
 265     fMerger = 0;
 266 }//AliRICH::AliRICH()
 267
 268 AliRICH::AliRICH(const char *name, const char *title)
 269     : AliDetector(name,title)
 270 {
 271 // Named ctor
 272    cout<<ClassName()<<"::named ctor(sName,sTitle)>\n"; // no way to control it as ctor is called before call to SetDebugXXXX()
 273
 274     fHits       = new TClonesArray("AliRICHHit",1000  );
 275     gAlice->AddHitList(fHits);
 276     fSDigits    = new TClonesArray("AliRICHSDigit",100000);
 277     fCerenkovs  = new TClonesArray("AliRICHCerenkov",1000);
 278     gAlice->AddHitList(fCerenkovs);
 279     fNSDigits   = 0;
 280     fNcerenkovs = 0;
 281     fIshunt     = 0;
 282
 283     //fNdch      = new Int_t[kNCH];
 284
 285     fDchambers = new TObjArray(kNCH);
 286
 287     fRecHits1D = new TObjArray(kNCH);
 288     fRecHits3D = new TObjArray(kNCH);
 289
 290     Int_t i;
 291
 292     for (i=0; i<kNCH ;i++) {
 293       //PH      (*fDchambers)[i] = new TClonesArray("AliRICHDigit",10000);
 294         fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i);
 295         fNdch[i]=0;
 296     }
 297
 298     //fNrawch      = new Int_t[kNCH];
 299
 300     fRawClusters = new TObjArray(kNCH);
 301     //printf("Created fRwClusters with adress:%p",fRawClusters);
 302
 303     for (i=0; i<kNCH ;i++) {
 304       //PH      (*fRawClusters)[i] = new TClonesArray("AliRICHRawCluster",10000);
 305       fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i);
 306       fNrawch[i]=0;
 307     }
 308
 309     //fNrechits      = new Int_t[kNCH];
 310
 311     for (i=0; i<kNCH ;i++) {
 312       //PH      (*fRecHits1D)[i] = new TClonesArray("AliRICHRecHit1D",1000);
 313       fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
 314     }
 315     for (i=0; i<kNCH ;i++) {
 316       //PH      (*fRecHits3D)[i] = new TClonesArray("AliRICHRecHit3D",1000);
 317       fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
 318     }
 319     //printf("Created fRecHits with adress:%p",fRecHits);
 320
 321
 322     SetMarkerColor(kRed);
 323
 324     /*fChambers = new TObjArray(kNCH);
 325     for (i=0; i<kNCH; i++)
 326       (*fChambers)[i] = new AliRICHChamber();*/
 327
 328     fFileName = 0;
 329     fMerger = 0;
 330 }
 331
 332 AliRICH::AliRICH(const AliRICH& RICH)
 333 {
 334 // Copy ctor
 335 }
 336
 337
 338 AliRICH::~AliRICH()
 339 {
 340 // Dtor of RICH manager class
 341    if(IsDebugStart()) cout<<ClassName()<<"::default dtor()>\n";
 342
 343     fIshunt  = 0;
 344     delete fHits;
 345     delete fSDigits;
 346     delete fCerenkovs;
 347
 348     //PH Delete TObjArrays
 349     if (fChambers) {
 350       fChambers->Delete();
 351       delete fChambers;
 352     }
 353     if (fDchambers) {
 354       fDchambers->Delete();
 355       delete fDchambers;
 356     }
 357     if (fRawClusters) {
 358       fRawClusters->Delete();
 359       delete fRawClusters;
 360     }
 361     if (fRecHits1D) {
 362       fRecHits1D->Delete();
 363       delete fRecHits1D;
 364     }
 365     if (fRecHits3D) {
 366       fRecHits3D->Delete();
 367       delete fRecHits3D;
 368     }
 369
 370 }
 371
 372
 373 //_____________________________________________________________________________
 374 Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
 375 {
 376 //  Calls the charge disintegration method of the current chamber and adds
 377 //  the simulated cluster to the root tree
 378    if(IsDebugHit()||IsDebugDigit()) cout<<ClassName()<<"::Hits2SDigits(...)>\n";
 379
 380    Int_t clhits[5];
 381    Float_t newclust[4][500];
 382    Int_t nnew;
 383
 384 //
 385 //  Integrated pulse height on chamber
 386
 387     clhits[0]=fNhits+1;
 388
 389     ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, nnew, newclust, res);
 390     Int_t ic=0;
 391
 392 //
 393 //  Add new clusters
 394     for (Int_t i=0; i<nnew; i++) {
 395         if (Int_t(newclust[0][i]) > 0) {
 396             ic++;
 397 //  Cluster Charge
 398             clhits[1] = Int_t(newclust[0][i]);
 399 //  Pad: ix
 400             clhits[2] = Int_t(newclust[1][i]);
 401 //  Pad: iy
 402             clhits[3] = Int_t(newclust[2][i]);
 403 //  Pad: chamber sector
 404             clhits[4] = Int_t(newclust[3][i]);
 405
 406             //printf(" %d %d %d %d %d\n",  clhits[0],  clhits[1],  clhits[2],  clhits[3],  clhits[4]);
 407
 408             AddSDigit(clhits);
 409         }
 410     }
 411
 412    if (gAlice->TreeS()){
 413         gAlice->TreeS()->Fill();
 414         gAlice->TreeS()->Write(0,TObject::kOverwrite);
 415         //printf("Filled SDigits...\n");
 416    }
 417
 418    return nnew;
 419 }//Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
 420
 421 void AliRICH::Hits2SDigits()
 422 {
 423 // Dummy: sdigits are created during transport.
 424 // Called from alirun.
 425    if(IsDebugHit()||IsDebugDigit()) cout<<ClassName()<<"::Hits2SDigits()>\n";
 426
 427
 428   int nparticles = gAlice->GetNtrack();
 429   cout << "Particles (RICH):" <<nparticles<<endl;
 430   if (nparticles > 0) printf("SDigits were already generated.\n");
 431
 432 }
 433
 434 //___________________________________________
 435 void AliRICH::SDigits2Digits(Int_t nev, Int_t flag)
 436 {
 437 // Generate digits.
 438 // Called from macro. Multiple events, more functionality.
 439    if(IsDebugDigit()) cout<<ClassName()<<"::SDigits2Digits()>\n";
 440
 441    //AliRICHChamber*       iChamber;
 442
 443    //printf("Generating tresholds...\n");
 444
 445    //for(Int_t i=0;i<7;i++) {
 446    //iChamber = &(Chamber(i));
 447    //iChamber->GenerateTresholds();
 448    //}
 449
 450    //int nparticles = gAlice->GetNtrack();
 451    //cout << "Particles (RICH):" <<nparticles<<endl;
 452    //if (nparticles <= 0) return;
 453    //if (!fMerger) {
 454    //fMerger = new AliRICHMerger();
 455    //}
 456
 457
 458    //fMerger->Init();
 459    //fMerger->Digitise(nev,flag);
 460
 461    AliRunDigitizer * manager = new AliRunDigitizer(1,1);
 462    manager->SetInputStream(0,"galice.root");
 463    //AliRICHDigitizer *dRICH  = new AliRICHDigitizer(manager);
 464    manager->Exec("deb");
 465
 466 }
 467 //___________________________________________
 468 void AliRICH::SDigits2Digits()
 469 {
 470   SDigits2Digits(0,0);
 471 }
 472 //___________________________________________
 473 void AliRICH::Digits2Reco()
 474 {
 475 // Generate clusters
 476 // Called from alirun, single event only.
 477    if(IsDebugDigit()||IsDebugReco()) cout<<ClassName()<<"::Digits2Reco()>\n";
 478
 479
 480   int nparticles = gAlice->GetNtrack();
 481   cout << "Particles (RICH):" <<nparticles<<endl;
 482   if (nparticles > 0) FindClusters(0,0);
 483
 484 }
 485
 486 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 487 {
 488 // Adds the current hit to the RICH hits list
 489    if(IsDebugHit()) cout<<ClassName()<<"::AddHit(...)>\n";
 490
 491     TClonesArray &lhits = *fHits;
 492     new(lhits[fNhits++]) AliRICHHit(fIshunt,track,vol,hits);
 493 }
 494
 495 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
 496 {
 497 // Adds a RICH cerenkov hit to the Cerenkov Hits list
 498    if(IsDebugHit()) cout<<ClassName()<<"::AddCerenkov()>\n";
 499
 500     TClonesArray &lcerenkovs = *fCerenkovs;
 501     new(lcerenkovs[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
 502 }
 503
 504 void AliRICH::AddSDigit(Int_t *aSDigit)
 505 {
 506 // Adds the current S digit to the RICH list of S digits
 507    if(IsDebugDigit()) cout<<ClassName()<<"::AddSDigit()>\n";
 508
 509   TClonesArray &lSDigits = *fSDigits;
 510   new(lSDigits[fNSDigits++]) AliRICHSDigit(aSDigit);
 511 }
 512
 513
 514 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
 515 {
 516 // Add a RICH digit to the list
 517    if(IsDebugDigit()) cout<<ClassName()<<"::AddDigit()>\n";
 518
 519    TClonesArray &ldigits = *((TClonesArray*)fDchambers->At(id));
 520    new(ldigits[fNdch[id]++]) AliRICHDigit(tracks,charges,digits);
 521 }
 522
 523 void AliRICH::AddRawCluster(Int_t id, const AliRICHRawCluster& c)
 524 {
 525 // Add a RICH digit to the list
 526
 527    if(IsDebugStart())
 528       cout<<ClassName()<<"::AddRawCluster()>\n";
 529
 530   //PH    TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
 531     TClonesArray &lrawcl = *((TClonesArray*)fRawClusters->At(id));
 532     new(lrawcl[fNrawch[id]++]) AliRICHRawCluster(c);
 533 }
 534
 535 //_____________________________________________________________________________
 536 void AliRICH::AddRecHit1D(Int_t id, Float_t *rechit, Float_t *photons, Int_t *padsx, Int_t* padsy)
 537 {
 538
 539   //
 540   // Add a RICH reconstructed hit to the list
 541   //
 542
 543   //PH    TClonesArray &lrec1D = *((TClonesArray*)(*fRecHits1D)[id]);
 544     TClonesArray &lrec1D = *((TClonesArray*)fRecHits1D->At(id));
 545     new(lrec1D[fNrechits1D[id]++]) AliRICHRecHit1D(id,rechit,photons,padsx,padsy);
 546 }
 547
 548 //_____________________________________________________________________________
 549 void AliRICH::AddRecHit3D(Int_t id, Float_t *rechit, Float_t omega, Float_t theta, Float_t phi)
 550 {
 551 // Add a RICH reconstructed hit to the list
 552
 553     TClonesArray &lrec3D = *((TClonesArray*)fRecHits3D->At(id));
 554     new(lrec3D[fNrechits3D[id]++]) AliRICHRecHit3D(id,rechit,omega,theta,phi);
 555 }
 556
 557 //___________________________________________
 558 void AliRICH::BuildGeometry()
 559
 560 {
 561
 562   //
 563   // Builds a TNode geometry for event display
 564   //
 565     TNode *node, *subnode, *top;
 566
 567     const int kColorRICH = kRed;
 568     //
 569     top=gAlice->GetGeometry()->GetNode("alice");
 570
 571     AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
 572     AliRICHSegmentationV0*  segmentation;
 573     AliRICHChamber*       iChamber;
 574     AliRICHGeometry*  geometry;
 575
 576     iChamber = &(pRICH->Chamber(0));
 577     segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel();
 578     geometry=iChamber->GetGeometryModel();
 579
 580     new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
 581
 582     Float_t padplane_width = segmentation->GetPadPlaneWidth();
 583     Float_t padplane_length = segmentation->GetPadPlaneLength();
 584
 585     //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());
 586
 587     new TBRIK("PHOTO","PHOTO","void", padplane_width/2,.1,padplane_length/2);
 588
 589     //printf("\n\n\n\n\n Padplane   w: %f l: %f \n\n\n\n\n", padplane_width/2,padplane_length/2);
 590     //printf("\n\n\n\n\n Padplane   w: %f l: %f \n\n\n\n\n", segmentation->GetPadPlaneWidth(), segmentation->GetPadPlaneLength());
 591
 592     Float_t offset       = 490 + 1.276 - geometry->GetGapThickness()/2;        //distance from center of mother volume to methane
 593     Float_t deltaphi     = 19.5;                                               //phi angle between center of chambers - z direction
 594     Float_t deltatheta   = 20;                                                 //theta angle between center of chambers - x direction
 595     Float_t cosphi       = TMath::Cos(deltaphi*TMath::Pi()/180);
 596     Float_t sinphi       = TMath::Sin(deltaphi*TMath::Pi()/180);
 597     Float_t costheta     = TMath::Cos(deltatheta*TMath::Pi()/180);
 598     Float_t sintheta     = TMath::Sin(deltatheta*TMath::Pi()/180);
 599
 600     //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
 601
 602     new TRotMatrix("rot993","rot993",90., 0.               , 90. - deltaphi, 90.             , deltaphi, -90.           );
 603     new TRotMatrix("rot994","rot994",90., -deltatheta      , 90.           , 90.- deltatheta , 0.      , 0.             );
 604     new TRotMatrix("rot995","rot995",90., 0.               , 90.           , 90.             , 0.      , 0.             );
 605     new TRotMatrix("rot996","rot996",90.,  deltatheta      , 90.           , 90 + deltatheta , 0.      , 0.             );
 606     new TRotMatrix("rot997","rot997",90., 360. - deltatheta, 108.2         , 90.- deltatheta ,18.2     , 90 - deltatheta);
 607     new TRotMatrix("rot998","rot998",90., 0.               , 90 + deltaphi , 90.             , deltaphi, 90.            );
 608     new TRotMatrix("rot999","rot999",90., deltatheta       , 108.2         , 90.+ deltatheta ,18.2     , 90 + deltatheta);
 609
 610     Float_t pos1[3]={0.                , offset*cosphi         , offset*sinphi};
 611     Float_t pos2[3]={offset*sintheta   , offset*costheta       , 0. };
 612     Float_t pos3[3]={0.                , offset                , 0.};
 613     Float_t pos4[3]={-offset*sintheta  , offset*costheta       , 0.};
 614     Float_t pos5[3]={offset*sinphi     , offset*costheta*cosphi, -offset*sinphi};
 615     Float_t pos6[3]={0.                , offset*cosphi         , -offset*sinphi};
 616     Float_t pos7[3]={ -offset*sinphi   , offset*costheta*cosphi, -offset*sinphi};
 617
 618
 619     top->cd();
 620     //Float_t pos1[3]={0,471.8999,165.2599};
 621     //Chamber(0).SetChamberTransform(pos1[0],pos1[1],pos1[2],
 622     //new TRotMatrix("rot993","rot993",90,0,70.69,90,19.30999,-90);
 623     node = new TNode("RICH1","RICH1","S_RICH",pos1[0],pos1[1],pos1[2],"rot993");
 624     node->SetLineColor(kColorRICH);
 625     node->cd();
 626     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 627     subnode->SetLineColor(kGreen);
 628     fNodes->Add(subnode);
 629     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 630     subnode->SetLineColor(kGreen);
 631     fNodes->Add(subnode);
 632     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 633     subnode->SetLineColor(kGreen);
 634     fNodes->Add(subnode);
 635     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 636     subnode->SetLineColor(kGreen);
 637     fNodes->Add(subnode);
 638     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 639     subnode->SetLineColor(kGreen);
 640     fNodes->Add(subnode);
 641     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 642     subnode->SetLineColor(kGreen);
 643     fNodes->Add(subnode);
 644     fNodes->Add(node);
 645
 646
 647     top->cd();
 648     //Float_t pos2[3]={171,470,0};
 649     //Chamber(1).SetChamberTransform(pos2[0],pos2[1],pos2[2],
 650     //new TRotMatrix("rot994","rot994",90,-20,90,70,0,0);
 651     node = new TNode("RICH2","RICH2","S_RICH",pos2[0],pos2[1],pos2[2],"rot994");
 652     node->SetLineColor(kColorRICH);
 653     node->cd();
 654     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 655     subnode->SetLineColor(kGreen);
 656     fNodes->Add(subnode);
 657     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 658     subnode->SetLineColor(kGreen);
 659     fNodes->Add(subnode);
 660     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 661     subnode->SetLineColor(kGreen);
 662     fNodes->Add(subnode);
 663     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 664     subnode->SetLineColor(kGreen);
 665     fNodes->Add(subnode);
 666     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 667     subnode->SetLineColor(kGreen);
 668     fNodes->Add(subnode);
 669     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 670     subnode->SetLineColor(kGreen);
 671     fNodes->Add(subnode);
 672     fNodes->Add(node);
 673
 674
 675     top->cd();
 676     //Float_t pos3[3]={0,500,0};
 677     //Chamber(2).SetChamberTransform(pos3[0],pos3[1],pos3[2],
 678     //new TRotMatrix("rot995","rot995",90,0,90,90,0,0);
 679     node = new TNode("RICH3","RICH3","S_RICH",pos3[0],pos3[1],pos3[2],"rot995");
 680     node->SetLineColor(kColorRICH);
 681     node->cd();
 682     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 683     subnode->SetLineColor(kGreen);
 684     fNodes->Add(subnode);
 685     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 686     subnode->SetLineColor(kGreen);
 687     fNodes->Add(subnode);
 688     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 689     subnode->SetLineColor(kGreen);
 690     fNodes->Add(subnode);
 691     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 692     subnode->SetLineColor(kGreen);
 693     fNodes->Add(subnode);
 694     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 695     subnode->SetLineColor(kGreen);
 696     fNodes->Add(subnode);
 697     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 698     subnode->SetLineColor(kGreen);
 699     fNodes->Add(subnode);
 700     fNodes->Add(node);
 701
 702     top->cd();
 703     //Float_t pos4[3]={-171,470,0};
 704     //Chamber(3).SetChamberTransform(pos4[0],pos4[1],pos4[2],
 705     //new TRotMatrix("rot996","rot996",90,20,90,110,0,0);
 706     node = new TNode("RICH4","RICH4","S_RICH",pos4[0],pos4[1],pos4[2],"rot996");
 707     node->SetLineColor(kColorRICH);
 708     node->cd();
 709     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 710     subnode->SetLineColor(kGreen);
 711     fNodes->Add(subnode);
 712     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 713     subnode->SetLineColor(kGreen);
 714     fNodes->Add(subnode);
 715     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 716     subnode->SetLineColor(kGreen);
 717     fNodes->Add(subnode);
 718     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 719     subnode->SetLineColor(kGreen);
 720     fNodes->Add(subnode);
 721     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 722     subnode->SetLineColor(kGreen);
 723     fNodes->Add(subnode);
 724     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 725     subnode->SetLineColor(kGreen);
 726     fNodes->Add(subnode);
 727     fNodes->Add(node);
 728
 729
 730     top->cd();
 731     //Float_t pos5[3]={161.3999,443.3999,-165.3};
 732     //Chamber(4).SetChamberTransform(pos5[0],pos5[1],pos5[2],
 733     //new TRotMatrix("rot997","rot997",90,340,108.1999,70,18.2,70);
 734     node = new TNode("RICH5","RICH5","S_RICH",pos5[0],pos5[1],pos5[2],"rot997");
 735     node->SetLineColor(kColorRICH);
 736     node->cd();
 737     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 738     subnode->SetLineColor(kGreen);
 739     fNodes->Add(subnode);
 740     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 741     subnode->SetLineColor(kGreen);
 742     fNodes->Add(subnode);
 743     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 744     subnode->SetLineColor(kGreen);
 745     fNodes->Add(subnode);
 746     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 747     subnode->SetLineColor(kGreen);
 748     fNodes->Add(subnode);
 749     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 750     subnode->SetLineColor(kGreen);
 751     fNodes->Add(subnode);
 752     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 753     subnode->SetLineColor(kGreen);
 754     fNodes->Add(subnode);
 755     fNodes->Add(node);
 756
 757
 758     top->cd();
 759     //Float_t pos6[3]={0., 471.9, -165.3,};
 760     //Chamber(5).SetChamberTransform(pos6[0],pos6[1],pos6[2],
 761     //new TRotMatrix("rot998","rot998",90,0,109.3099,90,19.30999,90);
 762     node = new TNode("RICH6","RICH6","S_RICH",pos6[0],pos6[1],pos6[2],"rot998");
 763     node->SetLineColor(kColorRICH);
 764     fNodes->Add(node);node->cd();
 765     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 766     subnode->SetLineColor(kGreen);
 767     fNodes->Add(subnode);
 768     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 769     subnode->SetLineColor(kGreen);
 770     fNodes->Add(subnode);
 771     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 772     subnode->SetLineColor(kGreen);
 773     fNodes->Add(subnode);
 774     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 775     subnode->SetLineColor(kGreen);
 776     fNodes->Add(subnode);
 777     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 778     subnode->SetLineColor(kGreen);
 779     fNodes->Add(subnode);
 780     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 781     subnode->SetLineColor(kGreen);
 782     fNodes->Add(subnode);
 783
 784
 785     top->cd();
 786     //Float_t pos7[3]={-161.399,443.3999,-165.3};
 787     //Chamber(6).SetChamberTransform(pos7[0],pos7[1],pos7[2],
 788     //new TRotMatrix("rot999","rot999",90,20,108.1999,110,18.2,110);
 789     node = new TNode("RICH7","RICH7","S_RICH",pos7[0],pos7[1],pos7[2],"rot999");
 790     node->SetLineColor(kColorRICH);
 791     node->cd();
 792     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 793     subnode->SetLineColor(kGreen);
 794     fNodes->Add(subnode);
 795     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,padplane_length/2 + segmentation->DeadZone()/2,"");
 796     subnode->SetLineColor(kGreen);
 797     fNodes->Add(subnode);
 798     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,padplane_length/2 + segmentation->DeadZone()/2,"");
 799     subnode->SetLineColor(kGreen);
 800     fNodes->Add(subnode);
 801     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",padplane_width + segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 802     subnode->SetLineColor(kGreen);
 803     fNodes->Add(subnode);
 804     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 805     subnode->SetLineColor(kGreen);
 806     fNodes->Add(subnode);
 807     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-padplane_width - segmentation->DeadZone(),5,-padplane_length/2 - segmentation->DeadZone()/2,"");
 808     subnode->SetLineColor(kGreen);
 809     fNodes->Add(subnode);
 810     fNodes->Add(node);
 811
 812 }
 813
 814 //___________________________________________
 815 void AliRICH::CreateGeometry()
 816 {
 817     //
 818     // Create the geometry for RICH version 1
 819     //
 820     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
 821     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
 822     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
 823     //
 824     //Begin_Html
 825     /*
 826       <img src="picts/AliRICHv1.gif">
 827     */
 828     //End_Html
 829     //Begin_Html
 830     /*
 831       <img src="picts/AliRICHv1Tree.gif">
 832     */
 833     //End_Html
 834
 835   AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
 836   AliRICHSegmentationV0*  segmentation;
 837   AliRICHGeometry*  geometry;
 838   AliRICHChamber*       iChamber;
 839
 840   iChamber = &(pRICH->Chamber(0));
 841   segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel();
 842   geometry=iChamber->GetGeometryModel();
 843
 844   Float_t distance;
 845   distance = geometry->GetFreonThickness()/2 + geometry->GetQuartzThickness() + geometry->GetGapThickness();
 846   geometry->SetRadiatorToPads(distance);
 847
 848   //Opaque quartz thickness
 849   Float_t oqua_thickness = .5;
 850   //CsI dimensions
 851
 852   //Float_t csi_length = 160*.8 + 2.6;
 853   //Float_t csi_width = 144*.84 + 2*2.6;
 854
 855   Float_t csi_width = segmentation->Npx()*segmentation->Dpx() + segmentation->DeadZone();
 856   Float_t csi_length = segmentation->Npy()*segmentation->Dpy() + 2*segmentation->DeadZone();
 857
 858   //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());
 859
 860   Int_t *idtmed = fIdtmed->GetArray()-999;
 861
 862     Int_t i;
 863     Float_t zs;
 864     Int_t idrotm[1099];
 865     Float_t par[3];
 866
 867     // --- Define the RICH detector
 868     //     External aluminium box
 869     par[0] = 68.8;
 870     par[1] = 13;                 //Original Settings
 871     par[2] = 70.86;
 872     /*par[0] = 73.15;
 873     par[1] = 11.5;
 874     par[2] = 71.1;*/
 875     gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 876
 877     //     Air
 878     par[0] = 66.3;
 879     par[1] = 13;                 //Original Settings
 880     par[2] = 68.35;
 881     /*par[0] = 66.55;
 882     par[1] = 11.5;
 883     par[2] = 64.8;*/
 884     gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], par, 3);
 885
 886     //    Air 2 (cutting the lower part of the box)
 887
 888     par[0] = 1.25;
 889     par[1] = 3;                 //Original Settings
 890     par[2] = 70.86;
 891     gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
 892
 893     //    Air 3 (cutting the lower part of the box)
 894
 895     par[0] = 66.3;
 896     par[1] = 3;                 //Original Settings
 897     par[2] = 1.2505;
 898     gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
 899
 900     //     Honeycomb
 901     par[0] = 66.3;
 902     par[1] = .188;                 //Original Settings
 903     par[2] = 68.35;
 904     /*par[0] = 66.55;
 905     par[1] = .188;
 906     par[2] = 63.1;*/
 907     gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
 908
 909     //     Aluminium sheet
 910     par[0] = 66.3;
 911     par[1] = .025;                 //Original Settings
 912     par[2] = 68.35;
 913     /*par[0] = 66.5;
 914     par[1] = .025;
 915     par[2] = 63.1;*/
 916     gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], par, 3);
 917
 918     //     Quartz
 919     par[0] = geometry->GetQuartzWidth()/2;
 920     par[1] = geometry->GetQuartzThickness()/2;
 921     par[2] = geometry->GetQuartzLength()/2;
 922     /*par[0] = 63.1;
 923     par[1] = .25;                  //Original Settings
 924     par[2] = 65.5;*/
 925     /*par[0] = geometry->GetQuartzWidth()/2;
 926     par[1] = geometry->GetQuartzThickness()/2;
 927     par[2] = geometry->GetQuartzLength()/2;*/
 928     //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]);
 929     gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
 930
 931     //     Spacers (cylinders)
 932     par[0] = 0.;
 933     par[1] = .5;
 934     par[2] = geometry->GetFreonThickness()/2;
 935     gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);
 936
 937     //     Feet (freon slabs supports)
 938
 939     par[0] = .7;
 940     par[1] = .3;
 941     par[2] = 1.9;
 942     gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
 943
 944     //     Opaque quartz
 945     par[0] = geometry->GetQuartzWidth()/2;
 946     par[1] = .2;
 947     par[2] = geometry->GetQuartzLength()/2;
 948     /*par[0] = 61.95;
 949     par[1] = .2;                   //Original Settings
 950     par[2] = 66.5;*/
 951     /*par[0] = 66.5;
 952     par[1] = .2;
 953     par[2] = 61.95;*/
 954     gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
 955
 956     //     Frame of opaque quartz
 957     par[0] = geometry->GetOuterFreonWidth()/2;
 958     //+ oqua_thickness;
 959     par[1] = geometry->GetFreonThickness()/2;
 960     par[2] = geometry->GetOuterFreonLength()/2;
 961     //+ oqua_thickness;
 962     /*par[0] = 20.65;
 963     par[1] = .5;                   //Original Settings
 964     par[2] = 66.5;*/
 965     /*par[0] = 66.5;
 966     par[1] = .5;
 967     par[2] = 20.65;*/
 968     gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
 969
 970     par[0] = geometry->GetInnerFreonWidth()/2;
 971     par[1] = geometry->GetFreonThickness()/2;
 972     par[2] = geometry->GetInnerFreonLength()/2;
 973     gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
 974
 975     //     Little bar of opaque quartz
 976     //par[0] = .275;
 977     //par[1] = geometry->GetQuartzThickness()/2;
 978     //par[2] = geometry->GetInnerFreonLength()/2 - 2.4;
 979     //par[2] = geometry->GetInnerFreonLength()/2;
 980     //+ oqua_thickness;
 981     /*par[0] = .275;
 982     par[1] = .25;                   //Original Settings
 983     par[2] = 63.1;*/
 984     /*par[0] = 63.1;
 985     par[1] = .25;
 986     par[2] = .275;*/
 987     //gMC->Gsvolu("BARR", "BOX ", idtmed[1007], par, 3);
 988
 989     //     Freon
 990     par[0] = geometry->GetOuterFreonWidth()/2 - oqua_thickness;
 991     par[1] = geometry->GetFreonThickness()/2;
 992     par[2] = geometry->GetOuterFreonLength()/2 - 2*oqua_thickness;
 993     /*par[0] = 20.15;
 994     par[1] = .5;                   //Original Settings
 995     par[2] = 65.5;*/
 996     /*par[0] = 65.5;
 997     par[1] = .5;
 998     par[2] = 20.15;*/
 999     gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
1000
1001     par[0] = geometry->GetInnerFreonWidth()/2 - oqua_thickness;
1002     par[1] = geometry->GetFreonThickness()/2;
1003     par[2] = geometry->GetInnerFreonLength()/2 - 2*oqua_thickness;
1004     gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);
1005
1006     //     Methane
1007     //par[0] = 64.8;
1008     par[0] = csi_width/2;
1009     par[1] = geometry->GetGapThickness()/2;
1010     //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]);
1011     //par[2] = 64.8;
1012     par[2] = csi_length/2;
1013     gMC->Gsvolu("META", "BOX ", idtmed[1004], par, 3);
1014
1015     //     Methane gap
1016     //par[0] = 64.8;
1017     par[0] = csi_width/2;
1018     par[1] = geometry->GetProximityGapThickness()/2;
1019     //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]);
1020     //par[2] = 64.8;
1021     par[2] = csi_length/2;
1022     gMC->Gsvolu("GAP ", "BOX ", idtmed[1008], par, 3);
1023
1024     //     CsI photocathode
1025     //par[0] = 64.8;
1026     par[0] = csi_width/2;
1027     par[1] = .25;
1028     //par[2] = 64.8;
1029     par[2] = csi_length/2;
1030     gMC->Gsvolu("CSI ", "BOX ", idtmed[1005], par, 3);
1031
1032     //     Anode grid
1033     par[0] = 0.;
1034     par[1] = .001;
1035     par[2] = 20.;
1036     gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
1037
1038     // Wire supports
1039     // Bar of metal
1040
1041     par[0] = csi_width/2;
1042     par[1] = 1.05;
1043     par[2] = 1.05;
1044     gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
1045
1046     // Ceramic pick up (base)
1047
1048     par[0] =  csi_width/2;
1049     par[1] = .25;
1050     par[2] = 1.05;
1051     gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
1052
1053     // Ceramic pick up (head)
1054
1055     par[0] = csi_width/2;
1056     par[1] = .1;
1057     par[2] = .1;
1058     gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
1059
1060     // Aluminium supports for methane and CsI
1061     // Short bar
1062
1063     par[0] = csi_width/2;
1064     par[1] = geometry->GetGapThickness()/2 + .25;
1065     par[2] = (68.35 - csi_length/2)/2;
1066     gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
1067
1068     // Long bar
1069
1070     par[0] = (66.3 - csi_width/2)/2;
1071     par[1] = geometry->GetGapThickness()/2 + .25;
1072     par[2] = csi_length/2 + 68.35 - csi_length/2;
1073     gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
1074
1075     // Aluminium supports for freon
1076     // Short bar
1077
1078     par[0] = geometry->GetQuartzWidth()/2;
1079     par[1] = .3;
1080     par[2] = (68.35 - geometry->GetQuartzLength()/2)/2;
1081     gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);
1082
1083     // Long bar
1084
1085     par[0] = (66.3 - geometry->GetQuartzWidth()/2)/2;
1086     par[1] = .3;
1087     par[2] = geometry->GetQuartzLength()/2 + 68.35 - geometry->GetQuartzLength()/2;
1088     gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);
1089
1090     // PCB backplane
1091
1092     par[0] = csi_width/2;
1093     par[1] = .25;
1094     par[2] = csi_length/4 -.5025;
1095     gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
1096
1097
1098     // Backplane supports
1099
1100     // Aluminium slab
1101
1102     par[0] = 33.15;
1103     par[1] = 2;
1104     par[2] = 21.65;
1105     gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);
1106
1107     // Big hole
1108
1109     par[0] = 9.05;
1110     par[1] = 2;
1111     par[2] = 4.4625;
1112     gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
1113
1114     // Small hole
1115
1116     par[0] = 5.7;
1117     par[1] = 2;
1118     par[2] = 4.4625;
1119     gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
1120
1121     // Place holes inside backplane support
1122
1123     gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
1124     gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
1125     gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
1126     gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
1127     gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1128     gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1129     gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1130     gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1131     gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
1132     gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
1133     gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1134     gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
1135     gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1136     gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
1137     gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1138     gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
1139
1140
1141
1142     // --- Places the detectors defined with GSVOLU
1143     //     Place material inside RICH
1144     gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
1145     gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1146     gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
1147     gMC->Gspos("AIR3", 1, "RICH", 0.,  1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
1148     gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35,  68.35 + 1.25, 0, "ONLY");
1149
1150
1151     gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
1152     gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- geometry->GetGapThickness()/2  - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
1153     gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
1154     gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1155     gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
1156     gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1157     gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, 36.9, 0, "ONLY");
1158     gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1159     gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1160     gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1161     gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .4 - .3, -36.9, 0, "ONLY");
1162     gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()- .2, 0., 0, "ONLY");
1163
1164     // Supports placing
1165
1166     // Methane supports
1167     gMC->Gspos("SMLG", 1, "SRIC", csi_width/2 + (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1168     gMC->Gspos("SMLG", 2, "SRIC", - csi_width/2 - (66.3 - csi_width/2)/2, 1.276 + .25, 0., 0, "ONLY");
1169     gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, csi_length/2 + (68.35 - csi_length/2)/2, 0, "ONLY");
1170     gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - csi_length/2 - (68.35 - csi_length/2)/2, 0, "ONLY");
1171
1172     //Freon supports
1173
1174     Float_t supp_y = 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness() - .2 + .3; //y position of freon supports
1175
1176     gMC->Gspos("SFLG", 1, "SRIC", geometry->GetQuartzWidth()/2 + (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1177     gMC->Gspos("SFLG", 2, "SRIC", - geometry->GetQuartzWidth()/2 - (66.3 - geometry->GetQuartzWidth()/2)/2, supp_y, 0., 0, "ONLY");
1178     gMC->Gspos("SFSH", 1, "SRIC", 0., supp_y, geometry->GetQuartzLength()/2 + (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1179     gMC->Gspos("SFSH", 2, "SRIC", 0., supp_y, - geometry->GetQuartzLength()/2 - (68.35 - geometry->GetQuartzLength()/2)/2, 0, "ONLY");
1180
1181     AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
1182
1183      //Placing of the spacers inside the freon slabs
1184
1185     Int_t nspacers = 30;
1186     //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);
1187
1188     //printf("Nspacers: %d", nspacers);
1189
1190     for (i = 0; i < nspacers/3; i++) {
1191         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1192         gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
1193     }
1194
1195     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1196         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1197         gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
1198     }
1199
1200     for (i = (nspacers*2)/3; i < nspacers; ++i) {
1201         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1202         gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1203     }
1204
1205     for (i = 0; i < nspacers/3; i++) {
1206         zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
1207         gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings
1208     }
1209
1210     for (i = nspacers/3; i < (nspacers*2)/3; i++) {
1211         zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1212         gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings
1213     }
1214
1215     for (i = (nspacers*2)/3; i < nspacers; ++i) {
1216         zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
1217         gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings
1218     }
1219
1220
1221     gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
1222     gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
1223     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)
1224 //    printf("Opaque quartz in SRIC %f\n", 1.276 - geometry->GetGapThickness()/2- geometry->GetQuartzThickness() -geometry->GetFreonThickness()/2);
1225     gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness() - geometry->GetFreonThickness()/2, 0., 0, "ONLY");          //Original settings
1226     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)
1227     //gMC->Gspos("BARR", 1, "QUAR", - geometry->GetInnerFreonWidth()/2 - oqua_thickness, 0., 0., 0, "ONLY");           //Original settings (-21.65)
1228     //gMC->Gspos("BARR", 2, "QUAR",  geometry->GetInnerFreonWidth()/2 + oqua_thickness, 0., 0., 0, "ONLY");            //Original settings (21.65)
1229     gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - geometry->GetGapThickness()/2 - geometry->GetQuartzThickness()/2, 0., 0, "ONLY");
1230     gMC->Gspos("GAP ", 1, "META", 0., geometry->GetGapThickness()/2 - geometry->GetProximityGapThickness()/2 - 0.0001, 0., 0, "ONLY");
1231     gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
1232     gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + geometry->GetGapThickness()/2 + .25, 0., 0, "ONLY");
1233     printf("CSI pos: %f\n",1.276 + geometry->GetGapThickness()/2 + .25);
1234
1235     // Wire support placing
1236
1237     gMC->Gspos("WSG2", 1, "GAP ", 0., geometry->GetProximityGapThickness()/2 - .1, 0., 0, "ONLY");
1238     gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
1239     gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
1240
1241     // Backplane placing
1242
1243     gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
1244     gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
1245     gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1246     gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
1247     gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1248     gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + geometry->GetGapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
1249
1250     // PCB placing
1251
1252     gMC->Gspos("PCB ", 1, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, csi_width/4 + .5025 + 2.5, 0, "ONLY");
1253     gMC->Gspos("PCB ", 2, "SRIC ", 0.,  1.276 + geometry->GetGapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
1254
1255
1256
1257     //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);
1258
1259     //     Place RICH inside ALICE apparatus
1260
1261     /* old values
1262
1263       AliMatrix(idrotm[1000], 90., 0., 70.69, 90., 19.31, -90.);
1264       AliMatrix(idrotm[1001], 90., -20., 90., 70., 0., 0.);
1265       AliMatrix(idrotm[1002], 90., 0., 90., 90., 0., 0.);
1266       AliMatrix(idrotm[1003], 90., 20., 90., 110., 0., 0.);
1267       AliMatrix(idrotm[1004], 90., 340., 108.2, 70., 18.2, 70.);
1268       AliMatrix(idrotm[1005], 90., 0., 109.31, 90., 19.31, 90.);
1269       AliMatrix(idrotm[1006], 90., 20., 108.2, 110., 18.2, 110.);
1270
1271       gMC->Gspos("RICH", 1, "ALIC", 0., 471.9, 165.26,     idrotm[1000], "ONLY");
1272       gMC->Gspos("RICH", 2, "ALIC", 171., 470., 0.,        idrotm[1001], "ONLY");
1273       gMC->Gspos("RICH", 3, "ALIC", 0., 500., 0.,          idrotm[1002], "ONLY");
1274       gMC->Gspos("RICH", 4, "ALIC", -171., 470., 0.,       idrotm[1003], "ONLY");
1275       gMC->Gspos("RICH", 5, "ALIC", 161.4, 443.4, -165.3,  idrotm[1004], "ONLY");
1276       gMC->Gspos("RICH", 6, "ALIC", 0., 471.9, -165.3,     idrotm[1005], "ONLY");
1277       gMC->Gspos("RICH", 7, "ALIC", -161.4, 443.4, -165.3, idrotm[1006], "ONLY");*/
1278
1279      // The placing of the chambers is measured from the vertex to the base of the methane vessel (490 cm)
1280
1281     Float_t offset       = 490 + 1.276 - geometry->GetGapThickness()/2;        //distance from center of mother volume to methane
1282     Float_t deltaphi     = 19.5;                                               //phi angle between center of chambers - z direction
1283     Float_t deltatheta   = 20;                                                 //theta angle between center of chambers - x direction
1284     Float_t cosphi       = TMath::Cos(deltaphi*TMath::Pi()/180);
1285     Float_t sinphi       = TMath::Sin(deltaphi*TMath::Pi()/180);
1286     Float_t costheta     = TMath::Cos(deltatheta*TMath::Pi()/180);
1287     Float_t sintheta     = TMath::Sin(deltatheta*TMath::Pi()/180);
1288
1289     //printf("\n\n%f %f %f %f %f %f %f\n\n",offset,deltatheta,deltaphi,cosphi,costheta,sinphi,sintheta);
1290
1291     AliMatrix(idrotm[1000], 90., 0.               , 90. - deltaphi, 90.             , deltaphi, -90.           );
1292     AliMatrix(idrotm[1001], 90., -deltatheta      , 90.           , 90.- deltatheta , 0.      , 0.             );
1293     AliMatrix(idrotm[1002], 90., 0.               , 90.           , 90.             , 0.      , 0.             );
1294     AliMatrix(idrotm[1003], 90.,  deltatheta      , 90.           , 90 + deltatheta , 0.      , 0.             );
1295     AliMatrix(idrotm[1004], 90., 360. - deltatheta, 108.2         , 90.- deltatheta ,18.2     , 90 - deltatheta);
1296     AliMatrix(idrotm[1005], 90., 0.               , 90 + deltaphi , 90.             , deltaphi, 90.            );
1297     AliMatrix(idrotm[1006], 90., deltatheta       , 108.2         , 90.+ deltatheta ,18.2     , 90 + deltatheta);
1298
1299     gMC->Gspos("RICH", 1, "ALIC", 0.                , offset*cosphi         , offset*sinphi ,idrotm[1000], "ONLY");
1300     gMC->Gspos("RICH", 2, "ALIC", (offset)*sintheta , offset*costheta       , 0.            ,idrotm[1001], "ONLY");
1301     gMC->Gspos("RICH", 3, "ALIC", 0.                , offset                , 0.            ,idrotm[1002], "ONLY");
1302     gMC->Gspos("RICH", 4, "ALIC", -(offset)*sintheta, offset*costheta       , 0.            ,idrotm[1003], "ONLY");
1303     gMC->Gspos("RICH", 5, "ALIC", (offset)*sinphi   , offset*costheta*cosphi, -offset*sinphi,idrotm[1004], "ONLY");
1304     gMC->Gspos("RICH", 6, "ALIC", 0.                , offset*cosphi         , -offset*sinphi,idrotm[1005], "ONLY");
1305     gMC->Gspos("RICH", 7, "ALIC", -(offset)*sinphi  , offset*costheta*cosphi, -offset*sinphi,idrotm[1006], "ONLY");
1306
1307 }
1308
1309
1310 //___________________________________________
1311 void AliRICH::CreateMaterials()
1312 {
1313     //
1314     // *** DEFINITION OF AVAILABLE RICH MATERIALS ***
1315     // ORIGIN    : NICK VAN EIJNDHOVEN
1316     // Modified by:  N. Colonna (INFN - BARI, Nicola.Colonna@ba.infn.it)
1317     //               R.A. Fini  (INFN - BARI, Rosanna.Fini@ba.infn.it)
1318     //               R.A. Loconsole (Bari University, loco@riscom.ba.infn.it)
1319     //
1320     Int_t   isxfld = gAlice->Field()->Integ();
1321     Float_t sxmgmx = gAlice->Field()->Max();
1322     Int_t i;
1323
1324     /************************************Antonnelo's Values (14-vectors)*****************************************/
1325     /*
1326     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,
1327                            6.7e-9,6.88e-9,7.08e-9,7.3e-9,7.51e-9,7.74e-9,8e-9 };
1328     Float_t rIndexQuarz[14] = { 1.528309,1.533333,
1329                                  1.538243,1.544223,1.550568,1.55777,
1330                                  1.565463,1.574765,1.584831,1.597027,
1331                                1.611858,1.6277,1.6472,1.6724 };
1332     Float_t rIndexOpaqueQuarz[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1333     Float_t rIndexMethane[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1334     Float_t rIndexGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1335     Float_t abscoFreon[14] = { 179.0987,179.0987,
1336                                 179.0987,179.0987,179.0987,142.92,56.65,13.95,10.43,7.07,2.03,.5773,.33496,0. };
1337     //Float_t abscoFreon[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1338         //                       1e-5,1e-5,1e-5,1e-5,1e-5 };
1339     Float_t abscoQuarz[14] = { 64.035,39.98,35.665,31.262,27.527,22.815,21.04,17.52,
1340                                 14.177,9.282,4.0925,1.149,.3627,.10857 };
1341     Float_t abscoOpaqueQuarz[14] = { 1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,
1342                                  1e-5,1e-5,1e-5,1e-5,1e-5 };
1343     Float_t abscoCsI[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1344                               1e-4,1e-4,1e-4,1e-4 };
1345     Float_t abscoMethane[14] = { 1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,1e6,
1346                                   1e6,1e6,1e6 };
1347     Float_t abscoGrid[14] = { 1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,1e-4,
1348                               1e-4,1e-4,1e-4,1e-4 };
1349     Float_t efficAll[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1350     Float_t efficCsI[14] = { 6e-4,.005,.0075,.01125,.045,.117,.135,.16575,
1351                               .17425,.1785,.1836,.1904,.1938,.221 };
1352     Float_t efficGrid[14] = { 1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1. };
1353     */
1354
1355
1356     /**********************************End of Antonnelo's Values**********************************/
1357
1358     /**********************************Values from rich_media.f (31-vectors)**********************************/
1359
1360
1361     //Photons energy intervals
1362     Float_t ppckov[26];
1363     for (i=0;i<26;i++)
1364     {
1365         ppckov[i] = (Float_t(i)*0.1+5.5)*1e-9;
1366         //printf ("Energy intervals: %e\n",ppckov[i]);
1367     }
1368
1369
1370     //Refraction index for quarz
1371     Float_t rIndexQuarz[26];
1372     Float_t  e1= 10.666;
1373     Float_t  e2= 18.125;
1374     Float_t  f1= 46.411;
1375     Float_t  f2= 228.71;
1376     for (i=0;i<26;i++)
1377     {
1378         Float_t ene=ppckov[i]*1e9;
1379         Float_t a=f1/(e1*e1 - ene*ene);
1380         Float_t b=f2/(e2*e2 - ene*ene);
1381         rIndexQuarz[i] = TMath::Sqrt(1. + a + b );
1382         //printf ("rIndexQuarz: %e\n",rIndexQuarz[i]);
1383     }
1384
1385     //Refraction index for opaque quarz, methane and grid
1386     Float_t rIndexOpaqueQuarz[26];
1387     Float_t rIndexMethane[26];
1388     Float_t rIndexGrid[26];
1389     for (i=0;i<26;i++)
1390     {
1391         rIndexOpaqueQuarz[i]=1;
1392         rIndexMethane[i]=1.000444;
1393         rIndexGrid[i]=1;
1394         //printf ("rIndexOpaqueQuarz , etc: %e, %e, %e\n",rIndexOpaqueQuarz[i], rIndexMethane[i], rIndexGrid[i]=1);
1395     }
1396
1397     //Absorption index for freon
1398     Float_t abscoFreon[26] = {179.0987, 179.0987, 179.0987, 179.0987, 179.0987,  179.0987, 179.0987, 179.0987,
1399                                179.0987, 142.9206, 56.64957, 25.58622, 13.95293, 12.03905, 10.42953, 8.804196,
1400                                7.069031, 4.461292, 2.028366, 1.293013, .577267,   .40746,  .334964, 0., 0., 0.};
1401
1402     //Absorption index for quarz
1403     /*Float_t Qzt [21] = {.0,.0,.005,.04,.35,.647,.769,.808,.829,.844,.853,.858,.869,.887,.903,.902,.902,
1404                         .906,.907,.907,.907};
1405     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,
1406                        215.0,220.0,225.0,230.0,235.0,240.0,245.0,250.0};
1407     Float_t abscoQuarz[31];
1408     for (Int_t i=0;i<31;i++)
1409     {
1410         Float_t Xlam = 1237.79 / (ppckov[i]*1e9);
1411         if (Xlam <= 160) abscoQuarz[i] = 0;
1412         if (Xlam > 250) abscoQuarz[i] = 1;
1413         else
1414         {
1415             for (Int_t j=0;j<21;j++)
1416             {
1417                 //printf ("Passed\n");
1418                 if (Xlam > Wavl2[j] && Xlam < Wavl2[j+1])
1419                 {
1420                     Float_t Dabs = (Qzt[j+1] - Qzt[j])/(Wavl2[j+1] - Wavl2[j]);
1421                     Float_t Abso = Qzt[j] + Dabs*(Xlam - Wavl2[j]);
1422                     abscoQuarz[i] = -5.0/(TMath::Log(Abso));
1423                 }
1424             }
1425         }
1426         printf ("abscoQuarz: %e abscoFreon: %e for energy: %e\n",abscoQuarz[i],abscoFreon[i],ppckov[i]);
1427     }*/
1428
1429     /*Float_t abscoQuarz[31] = {49.64211, 48.41296, 47.46989, 46.50492, 45.13682, 44.47883, 43.1929 , 41.30922, 40.5943 ,
1430                                39.82956, 38.98623, 38.6247 , 38.43448, 37.41084, 36.22575, 33.74852, 30.73901, 24.25086,
1431                                17.94531, 11.88753, 5.99128,  3.83503,  2.36661,  1.53155, 1.30582, 1.08574, .8779708,
1432                                .675275, 0., 0., 0.};
1433
1434     for (Int_t i=0;i<31;i++)
1435     {
1436         abscoQuarz[i] = abscoQuarz[i]/10;
1437     }*/
1438
1439     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,
1440                                 19.266, 17.525, 15.878, 14.177, 11.719, 9.282, 6.62, 4.0925, 2.601, 1.149, .667, .3627,
1441                                 .192, .1497, .10857};
1442
1443     //Absorption index for methane
1444     Float_t abscoMethane[26];
1445     for (i=0;i<26;i++)
1446     {
1447         abscoMethane[i]=AbsoCH4(ppckov[i]*1e9);
1448         //printf("abscoMethane: %e for energy: %e\n", abscoMethane[i],ppckov[i]*1e9);
1449     }
1450
1451     //Absorption index for opaque quarz, csi and grid, efficiency for all and grid
1452     Float_t abscoOpaqueQuarz[26];
1453     Float_t abscoCsI[26];
1454     Float_t abscoGrid[26];
1455     Float_t efficAll[26];
1456     Float_t efficGrid[26];
1457     for (i=0;i<26;i++)
1458     {
1459         abscoOpaqueQuarz[i]=1e-5;
1460         abscoCsI[i]=1e-4;
1461         abscoGrid[i]=1e-4;
1462         efficAll[i]=1;
1463         efficGrid[i]=1;
1464         //printf ("All must be 1: %e,  %e,  %e,  %e,  %e\n",abscoOpaqueQuarz[i],abscoCsI[i],abscoGrid[i],efficAll[i],efficGrid[i]);
1465     }
1466
1467     //Efficiency for csi
1468
1469     Float_t efficCsI[26] = {0.000199999995, 0.000600000028, 0.000699999975, 0.00499999989, 0.00749999983, 0.010125,
1470                              0.0242999997, 0.0405000001, 0.0688500032, 0.105299994, 0.121500008, 0.141749993, 0.157949999,
1471                              0.162, 0.166050002, 0.167669997, 0.174299985, 0.176789999, 0.179279998, 0.182599992, 0.18592,
1472                              0.187579989, 0.189239994, 0.190899998, 0.207499996, 0.215799987};
1473
1474
1475
1476     //FRESNEL LOSS CORRECTION FOR PERPENDICULAR INCIDENCE AND
1477     //UNPOLARIZED PHOTONS
1478
1479     for (i=0;i<26;i++)
1480     {
1481         efficCsI[i] = efficCsI[i]/(1.-Fresnel(ppckov[i]*1e9,1.,0));
1482         //printf ("Fresnel result: %e for energy: %e\n",Fresnel(ppckov[i]*1e9,1.,0),ppckov[i]*1e9);
1483     }
1484
1485     /*******************************************End of rich_media.f***************************************/
1486
1487
1488
1489
1490
1491
1492     Float_t afre[2], agri, amet[2], aqua[2], ahon, zfre[2], zgri, zhon,
1493     zmet[2], zqua[2];
1494     Int_t nlmatfre;
1495     Float_t densquao;
1496     Int_t nlmatmet, nlmatqua;
1497     Float_t wmatquao[2], rIndexFreon[26];
1498     Float_t aquao[2], epsil, stmin, zquao[2];
1499     Int_t nlmatquao;
1500     Float_t radlal, densal, tmaxfd, deemax, stemax;
1501     Float_t aal, zal, radlgri, densfre, radlhon, densgri, denshon,densqua, densmet, wmatfre[2], wmatmet[2], wmatqua[2];
1502
1503     Int_t *idtmed = fIdtmed->GetArray()-999;
1504
1505     // --- Photon energy (GeV)
1506     // --- Refraction indexes
1507     for (i = 0; i < 26; ++i) {
1508       rIndexFreon[i] = ppckov[i] * .0172 * 1e9 + 1.177;
1509       //rIndexFreon[i] = 1;
1510         //printf ("rIndexFreon: %e \n efficCsI: %e for energy: %e\n",rIndexFreon[i], efficCsI[i], ppckov[i]);
1511     }
1512
1513     // --- Detection efficiencies (quantum efficiency for CsI)
1514     // --- Define parameters for honeycomb.
1515     //     Used carbon of equivalent rad. lenght
1516
1517     ahon    = 12.01;
1518     zhon    = 6.;
1519     denshon = 0.1;
1520     radlhon = 18.8;
1521
1522     // --- Parameters to include in GSMIXT, relative to Quarz (SiO2)
1523
1524     aqua[0]    = 28.09;
1525     aqua[1]    = 16.;
1526     zqua[0]    = 14.;
1527     zqua[1]    = 8.;
1528     densqua    = 2.64;
1529     nlmatqua   = -2;
1530     wmatqua[0] = 1.;
1531     wmatqua[1] = 2.;
1532
1533     // --- Parameters to include in GSMIXT, relative to opaque Quarz (SiO2)
1534
1535     aquao[0]    = 28.09;
1536     aquao[1]    = 16.;
1537     zquao[0]    = 14.;
1538     zquao[1]    = 8.;
1539     densquao    = 2.64;
1540     nlmatquao   = -2;
1541     wmatquao[0] = 1.;
1542     wmatquao[1] = 2.;
1543
1544     // --- Parameters to include in GSMIXT, relative to Freon (C6F14)
1545
1546     afre[0]    = 12.;
1547     afre[1]    = 19.;
1548     zfre[0]    = 6.;
1549     zfre[1]    = 9.;
1550     densfre    = 1.7;
1551     nlmatfre   = -2;
1552     wmatfre[0] = 6.;
1553     wmatfre[1] = 14.;
1554
1555     // --- Parameters to include in GSMIXT, relative to methane (CH4)
1556
1557     amet[0]    = 12.01;
1558     amet[1]    = 1.;
1559     zmet[0]    = 6.;
1560     zmet[1]    = 1.;
1561     densmet    = 7.17e-4;
1562     nlmatmet   = -2;
1563     wmatmet[0] = 1.;
1564     wmatmet[1] = 4.;
1565
1566     // --- Parameters to include in GSMIXT, relative to anode grid (Cu)
1567
1568     agri    = 63.54;
1569     zgri    = 29.;
1570     densgri = 8.96;
1571     radlgri = 1.43;
1572
1573     // --- Parameters to include in GSMATE related to aluminium sheet
1574
1575     aal    = 26.98;
1576     zal    = 13.;
1577     densal = 2.7;
1578     radlal = 8.9;
1579
1580     // --- Glass parameters
1581
1582     Float_t aglass[5]={12.01, 28.09, 16.,   10.8,  23.};
1583     Float_t zglass[5]={ 6.,   14.,    8.,    5.,   11.};
1584     Float_t wglass[5]={ 0.5,  0.105, 0.355, 0.03,  0.01};
1585     Float_t dglass=1.74;
1586
1587
1588     AliMaterial(1, "Air     $", 14.61, 7.3, .001205, 30420., 67500);
1589     AliMaterial(6, "HON", ahon, zhon, denshon, radlhon, 0);
1590     AliMaterial(16, "CSI", ahon, zhon, denshon, radlhon, 0);
1591     AliMixture(20, "QUA", aqua, zqua, densqua, nlmatqua, wmatqua);
1592     AliMixture(21, "QUAO", aquao, zquao, densquao, nlmatquao, wmatquao);
1593     AliMixture(30, "FRE", afre, zfre, densfre, nlmatfre, wmatfre);
1594     AliMixture(40, "MET", amet, zmet, densmet, nlmatmet, wmatmet);
1595     AliMixture(41, "METG", amet, zmet, densmet, nlmatmet, wmatmet);
1596     AliMaterial(11, "GRI", agri, zgri, densgri, radlgri, 0);
1597     AliMaterial(50, "ALUM", aal, zal, densal, radlal, 0);
1598     AliMixture(32, "GLASS",aglass, zglass, dglass, 5, wglass);
1599     AliMaterial(31, "COPPER$",   63.54,    29.,   8.96,  1.4, 0.);
1600
1601     tmaxfd = -10.;
1602     stemax = -.1;
1603     deemax = -.2;
1604     epsil  = .001;
1605     stmin  = -.001;
1606
1607     AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1608     AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1609     AliMedium(3, "QUARZO$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1610     AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1611     AliMedium(5, "METANO$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1612     AliMedium(6, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
1613     AliMedium(7, "GRIGLIA$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1614     AliMedium(8, "QUARZOO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1615     AliMedium(9, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, .1, -deemax, epsil, -stmin);
1616     AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1617     AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1618     AliMedium(12, "PCB_COPPER", 31, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
1619
1620
1621     gMC->SetCerenkov(idtmed[1000], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1622     gMC->SetCerenkov(idtmed[1001], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1623     gMC->SetCerenkov(idtmed[1002], 26, ppckov, abscoQuarz, efficAll,rIndexQuarz);
1624     gMC->SetCerenkov(idtmed[1003], 26, ppckov, abscoFreon, efficAll,rIndexFreon);
1625     gMC->SetCerenkov(idtmed[1004], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1626     gMC->SetCerenkov(idtmed[1005], 26, ppckov, abscoCsI, efficCsI, rIndexMethane);
1627     gMC->SetCerenkov(idtmed[1006], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1628     gMC->SetCerenkov(idtmed[1007], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1629     gMC->SetCerenkov(idtmed[1008], 26, ppckov, abscoMethane, efficAll, rIndexMethane);
1630     gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
1631     gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
1632 }
1633
1634 //___________________________________________
1635
1636 Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
1637 {
1638
1639     //ENE(EV), PDOTI=COS(INC.ANG.), PDOTR=COS(POL.PLANE ROT.ANG.)
1640
1641     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,
1642                       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,
1643                       7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
1644
1645
1646     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
1647                         2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
1648                         2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
1649                         1.72,1.53};
1650
1651     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
1652                         0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
1653                         0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
1654                         1.714,1.498};
1655     Float_t xe=ene;
1656     Int_t  j=Int_t(xe*10)-49;
1657     Float_t cn=csin[j]+((csin[j+1]-csin[j])/0.1)*(xe-en[j]);
1658     Float_t ck=csik[j]+((csik[j+1]-csik[j])/0.1)*(xe-en[j]);
1659
1660     //FORMULAE FROM HANDBOOK OF OPTICS, 33.23 OR
1661     //W.R. HUNTER, J.O.S.A. 54 (1964),15 , J.O.S.A. 55(1965),1197
1662
1663     Float_t sinin=TMath::Sqrt(1-pdoti*pdoti);
1664     Float_t tanin=sinin/pdoti;
1665
1666     Float_t c1=cn*cn-ck*ck-sinin*sinin;
1667     Float_t c2=4*cn*cn*ck*ck;
1668     Float_t aO=TMath::Sqrt(0.5*(TMath::Sqrt(c1*c1+c2)+c1));
1669     Float_t b2=0.5*(TMath::Sqrt(c1*c1+c2)-c1);
1670
1671     Float_t rs=((aO-pdoti)*(aO-pdoti)+b2)/((aO+pdoti)*(aO+pdoti)+b2);
1672     Float_t rp=rs*((aO-sinin*tanin)*(aO-sinin*tanin)+b2)/((aO+sinin*tanin)*(aO+sinin*tanin)+b2);
1673
1674
1675     //CORRECTION FACTOR FOR SURFACE ROUGHNESS
1676     //B.J. STAGG  APPLIED OPTICS, 30(1991),4113
1677
1678     Float_t sigraf=18.;
1679     Float_t lamb=1240/ene;
1680     Float_t fresn;
1681
1682     Float_t  rO=TMath::Exp(-(4*TMath::Pi()*pdoti*sigraf/lamb)*(4*TMath::Pi()*pdoti*sigraf/lamb));
1683
1684     if(pola)
1685     {
1686         Float_t pdotr=0.8;                                 //DEGREE OF POLARIZATION : 1->P , -1->S
1687         fresn=0.5*(rp*(1+pdotr)+rs*(1-pdotr));
1688     }
1689     else
1690         fresn=0.5*(rp+rs);
1691
1692     fresn = fresn*rO;
1693     return(fresn);
1694 }
1695
1696 //__________________________________________
1697 Float_t AliRICH::AbsoCH4(Float_t x)
1698 {
1699
1700     //KLOSCH,SCH4(9),WL(9),EM(9),ALENGTH(31)
1701     Float_t sch4[9] = {.12,.16,.23,.38,.86,2.8,7.9,28.,80.};              //MB X 10^22
1702     //Float_t wl[9] = {153.,152.,151.,150.,149.,148.,147.,146.,145};
1703     Float_t em[9] = {8.1,8.158,8.212,8.267,8.322,8.378,8.435,8.493,8.55};
1704     const Float_t kLosch=2.686763E19;                                      // LOSCHMIDT NUMBER IN CM-3
1705     const Float_t kIgas1=100, kIgas2=0, kOxy=10., kWater=5., kPressure=750.,kTemperature=283.;
1706     Float_t pn=kPressure/760.;
1707     Float_t tn=kTemperature/273.16;
1708
1709
1710 // ------- METHANE CROSS SECTION -----------------
1711 // ASTROPH. J. 214, L47 (1978)
1712
1713     Float_t sm=0;
1714     if (x<7.75)
1715         sm=.06e-22;
1716
1717     if(x>=7.75 && x<=8.1)
1718     {
1719         Float_t c0=-1.655279e-1;
1720         Float_t c1=6.307392e-2;
1721         Float_t c2=-8.011441e-3;
1722         Float_t c3=3.392126e-4;
1723         sm=(c0+c1*x+c2*x*x+c3*x*x*x)*1.e-18;
1724     }
1725
1726     if (x> 8.1)
1727     {
1728         Int_t j=0;
1729         while (x<=em[j] && x>=em[j+1])
1730         {
1731             j++;
1732             Float_t a=(sch4[j+1]-sch4[j])/(em[j+1]-em[j]);
1733             sm=(sch4[j]+a*(x-em[j]))*1e-22;
1734         }
1735     }
1736
1737     Float_t dm=(kIgas1/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1738     Float_t abslm=1./sm/dm;
1739
1740 //    ------- ISOBUTHANE CROSS SECTION --------------
1741 //     i-C4H10 (ai) abs. length from curves in
1742 //     Lu-McDonald paper for BARI RICH workshop .
1743 //     -----------------------------------------------------------
1744
1745     Float_t ai;
1746     Float_t absli;
1747     if (kIgas2 != 0)
1748     {
1749         if (x<7.25)
1750             ai=100000000.;
1751
1752         if(x>=7.25 && x<7.375)
1753             ai=24.3;
1754
1755         if(x>=7.375)
1756             ai=.0000000001;
1757
1758         Float_t si = 1./(ai*kLosch*273.16/293.);                    // ISOB. CRO.SEC.IN CM2
1759         Float_t di=(kIgas2/100.)*(1.-((kOxy+kWater)/1.e6))*kLosch*pn/tn;
1760         absli =1./si/di;
1761     }
1762     else
1763         absli=1.e18;
1764 //    ---------------------------------------------------------
1765 //
1766 //       transmission of O2
1767 //
1768 //       y= path in cm, x=energy in eV
1769 //       so= cross section for UV absorption in cm2
1770 //       do= O2 molecular density in cm-3
1771 //    ---------------------------------------------------------
1772
1773     Float_t abslo;
1774     Float_t so=0;
1775     if(x>=6.0)
1776     {
1777         if(x>=6.0 && x<6.5)
1778         {
1779             so=3.392709e-13 * TMath::Exp(2.864104 *x);
1780             so=so*1e-18;
1781         }
1782
1783         if(x>=6.5 && x<7.0)
1784         {
1785             so=2.910039e-34 * TMath::Exp(10.3337*x);
1786             so=so*1e-18;
1787         }
1788
1789
1790         if (x>=7.0)
1791         {
1792             Float_t a0=-73770.76;
1793             Float_t a1=46190.69;
1794             Float_t a2=-11475.44;
1795             Float_t a3=1412.611;
1796             Float_t a4=-86.07027;
1797             Float_t a5=2.074234;
1798             so= a0+(a1*x)+(a2*x*x)+(a3*x*x*x)+(a4*x*x*x*x)+(a5*x*x*x*x*x);
1799             so=so*1e-18;
1800         }
1801
1802         Float_t dox=(kOxy/1e6)*kLosch*pn/tn;
1803         abslo=1./so/dox;
1804     }
1805     else
1806         abslo=1.e18;
1807 //     ---------------------------------------------------------
1808 //
1809 //       transmission of H2O
1810 //
1811 //       y= path in cm, x=energy in eV
1812 //       sw= cross section for UV absorption in cm2
1813 //       dw= H2O molecular density in cm-3
1814 //     ---------------------------------------------------------
1815
1816     Float_t abslw;
1817
1818     Float_t b0=29231.65;
1819     Float_t b1=-15807.74;
1820     Float_t b2=3192.926;
1821     Float_t b3=-285.4809;
1822     Float_t b4=9.533944;
1823
1824     if(x>6.75)
1825     {
1826         Float_t sw= b0+(b1*x)+(b2*x*x)+(b3*x*x*x)+(b4*x*x*x*x);
1827         sw=sw*1e-18;
1828         Float_t dw=(kWater/1e6)*kLosch*pn/tn;
1829         abslw=1./sw/dw;
1830     }
1831     else
1832         abslw=1.e18;
1833
1834 //    ---------------------------------------------------------
1835
1836     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
1837     return (alength);
1838 }
1839
1840
1841
1842 //___________________________________________
1843 Int_t AliRICH::DistancetoPrimitive(Int_t , Int_t )
1844 {
1845
1846 // Default value
1847
1848     return 9999;
1849 }
1850
1851 //___________________________________________
1852 void AliRICH::MakeBranch(Option_t* option, const char *file)
1853 {
1854   // Create Tree branches for the RICH.
1855
1856     const Int_t kBufferSize = 4000;
1857     char branchname[20];
1858
1859     AliDetector::MakeBranch(option,file);
1860
1861     const char *cH = strstr(option,"H");
1862     const char *cD = strstr(option,"D");
1863     const char *cR = strstr(option,"R");
1864     const char *cS = strstr(option,"S");
1865
1866
1867     if (cH) {
1868       sprintf(branchname,"%sCerenkov",GetName());
1869       if (fCerenkovs   && gAlice->TreeH()) {
1870         //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1871         MakeBranchInTree(gAlice->TreeH(),branchname, &fCerenkovs, kBufferSize, file) ;
1872         //branch->SetAutoDelete(kFALSE);
1873       }
1874       sprintf(branchname,"%sSDigits",GetName());
1875       if (fSDigits   && gAlice->TreeH()) {
1876         //TBranch* branch = MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1877         MakeBranchInTree(gAlice->TreeH(),branchname, &fSDigits, kBufferSize, file) ;
1878         //branch->SetAutoDelete(kFALSE);
1879         //printf("Making branch %sSDigits in TreeH\n",GetName());
1880       }
1881     }
1882
1883     if (cS) {
1884       sprintf(branchname,"%sSDigits",GetName());
1885       if (fSDigits   && gAlice->TreeS()) {
1886         //TBranch* branch = MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1887         MakeBranchInTree(gAlice->TreeS(),branchname, &fSDigits, kBufferSize, file) ;
1888         //branch->SetAutoDelete(kFALSE);
1889         //printf("Making branch %sSDigits in TreeS\n",GetName());
1890       }
1891     }
1892
1893     if (cD) {
1894     //
1895     // one branch for digits per chamber
1896     //
1897       Int_t i;
1898
1899       for (i=0; i<kNCH ;i++) {
1900         sprintf(branchname,"%sDigits%d",GetName(),i+1);
1901         if (fDchambers   && gAlice->TreeD()) {
1902            //TBranch* branch = MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1903           MakeBranchInTree(gAlice->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, file) ;
1904            //branch->SetAutoDelete(kFALSE);
1905           //printf("Making Branch %sDigits%d\n",GetName(),i+1);
1906         }
1907       }
1908     }
1909
1910     if (cR) {
1911     //
1912     // one branch for raw clusters per chamber
1913     //
1914
1915       //printf("Called MakeBranch for TreeR\n");
1916
1917       Int_t i;
1918
1919       for (i=0; i<kNCH ;i++) {
1920         sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
1921         if (fRawClusters && gAlice->TreeR()) {
1922            //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1923           MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, file) ;
1924            //branch->SetAutoDelete(kFALSE);
1925         }
1926       }
1927      //
1928      // one branch for rec hits per chamber
1929      //
1930      for (i=0; i<kNCH ;i++) {
1931        sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
1932        if (fRecHits1D   && gAlice->TreeR()) {
1933          //TBranch* branch = MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1934          MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, file) ;
1935          //branch->SetAutoDelete(kFALSE);
1936        }
1937      }
1938      for (i=0; i<kNCH ;i++) {
1939        sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
1940        if (fRecHits3D   && gAlice->TreeR()) {
1941          MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, file) ;
1942          //branch->SetAutoDelete(kFALSE);
1943       }
1944     }
1945   }
1946 }
1947
1948 //___________________________________________
1949 void AliRICH::SetTreeAddress()
1950 {
1951   // Set branch address for the Hits and Digits Tree.
1952   char branchname[20];
1953   Int_t i;
1954
1955     AliDetector::SetTreeAddress();
1956
1957     TBranch *branch;
1958     TTree *treeH = gAlice->TreeH();
1959     TTree *treeD = gAlice->TreeD();
1960     TTree *treeR = gAlice->TreeR();
1961     TTree *treeS = gAlice->TreeS();
1962
1963     if (treeH) {
1964       if (fCerenkovs) {
1965             branch = treeH->GetBranch("RICHCerenkov");
1966             if (branch) branch->SetAddress(&fCerenkovs);
1967         }
1968     if (fSDigits) {
1969         branch = treeH->GetBranch("RICHSDigits");
1970         if (branch)
1971           {
1972             branch->SetAddress(&fSDigits);
1973             //printf("Setting sdigits branch address at %p in TreeH\n",&fSDigits);
1974           }
1975       }
1976     }
1977
1978     if (treeS) {
1979       if (fSDigits) {
1980         branch = treeS->GetBranch("RICHSDigits");
1981         if (branch)
1982           {
1983             branch->SetAddress(&fSDigits);
1984             //printf("Setting sdigits branch address at %p in TreeS\n",&fSDigits);
1985           }
1986       }
1987     }
1988
1989
1990     if (treeD) {
1991         for (int i=0; i<kNCH; i++) {
1992             sprintf(branchname,"%sDigits%d",GetName(),i+1);
1993             if (fDchambers) {
1994               branch = treeD->GetBranch(branchname);
1995               if (branch) branch->SetAddress(&((*fDchambers)[i]));
1996             }
1997         }
1998     }
1999   if (treeR) {
2000       for (i=0; i<kNCH; i++) {
2001           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
2002           if (fRawClusters) {
2003               branch = treeR->GetBranch(branchname);
2004               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
2005           }
2006       }
2007
2008       for (i=0; i<kNCH; i++) {
2009         sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
2010         if (fRecHits1D) {
2011           branch = treeR->GetBranch(branchname);
2012           if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
2013           }
2014       }
2015
2016      for (i=0; i<kNCH; i++) {
2017         sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
2018         if (fRecHits3D) {
2019           branch = treeR->GetBranch(branchname);
2020           if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
2021           }
2022       }
2023
2024   }
2025 }
2026 //___________________________________________
2027 void AliRICH::ResetHits()
2028 {
2029   // Reset number of clusters and the cluster array for this detector
2030     AliDetector::ResetHits();
2031     fNSDigits   = 0;
2032     fNcerenkovs = 0;
2033     if (fSDigits)  fSDigits->Clear();
2034     if (fCerenkovs) fCerenkovs->Clear();
2035 }
2036
2037
2038 //____________________________________________
2039 void AliRICH::ResetDigits()
2040 {
2041   //
2042   // Reset number of digits and the digits array for this detector
2043   //
2044     for ( int i=0;i<kNCH;i++ ) {
2045       //PH      if (fDchambers && (*fDchambers)[i])   (*fDchambers)[i]->Clear();
2046         if (fDchambers && fDchambers->At(i))   fDchambers->At(i)->Clear();
2047         if (fNdch)  fNdch[i]=0;
2048     }
2049 }
2050
2051 //____________________________________________
2052 void AliRICH::ResetRawClusters()
2053 {
2054   //
2055   // Reset number of raw clusters and the raw clust array for this detector
2056   //
2057     for ( int i=0;i<kNCH;i++ ) {
2058       //PH      if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
2059         if (fRawClusters->At(i))    ((TClonesArray*)fRawClusters->At(i))->Clear();
2060         if (fNrawch)  fNrawch[i]=0;
2061     }
2062 }
2063
2064 //____________________________________________
2065 void AliRICH::ResetRecHits1D()
2066 {
2067   //
2068   // Reset number of raw clusters and the raw clust array for this detector
2069   //
2070
2071   for ( int i=0;i<kNCH;i++ ) {
2072     //PH        if ((*fRecHits1D)[i])    ((TClonesArray*)(*fRecHits1D)[i])->Clear();
2073         if (fRecHits1D->At(i))    ((TClonesArray*)fRecHits1D->At(i))->Clear();
2074         if (fNrechits1D)  fNrechits1D[i]=0;
2075     }
2076 }
2077
2078 //____________________________________________
2079 void AliRICH::ResetRecHits3D()
2080 {
2081   //
2082   // Reset number of raw clusters and the raw clust array for this detector
2083   //
2084
2085   for ( int i=0;i<kNCH;i++ ) {
2086     //PH        if ((*fRecHits3D)[i])    ((TClonesArray*)(*fRecHits3D)[i])->Clear();
2087         if (fRecHits3D->At(i))    ((TClonesArray*)fRecHits3D->At(i))->Clear();
2088         if (fNrechits3D)  fNrechits3D[i]=0;
2089     }
2090 }
2091
2092
2093 //___________________________________________
2094 void AliRICH::StepManager()
2095 {
2096 // Full Step Manager
2097
2098     Int_t          copy, id;
2099     static Int_t   idvol;
2100     static Int_t   vol[2];
2101     Int_t          ipart;
2102     static Float_t hits[22];
2103     static Float_t ckovData[19];
2104     TLorentzVector position;
2105     TLorentzVector momentum;
2106     Float_t        pos[3];
2107     Float_t        mom[4];
2108     Float_t        localPos[3];
2109     Float_t        localMom[4];
2110     Float_t        localTheta,localPhi;
2111     Float_t        theta,phi;
2112     Float_t        destep, step;
2113     Double_t        ranf[2];
2114     Int_t          nPads;
2115     Float_t        coscerenkov;
2116     static Float_t eloss, xhit, yhit, tlength;
2117     const  Float_t kBig=1.e10;
2118
2119     TClonesArray &lhits = *fHits;
2120     TParticle *current = (TParticle*)(*gAlice->Particles())[gAlice->CurrentTrack()];
2121
2122  //if (current->Energy()>1)
2123    //{
2124
2125     // Only gas gap inside chamber
2126     // Tag chambers and record hits when track enters
2127
2128
2129     id=gMC->CurrentVolID(copy);
2130     idvol = copy-1;
2131     Float_t cherenkovLoss=0;
2132     //gAlice->KeepTrack(gAlice->CurrentTrack());
2133
2134     gMC->TrackPosition(position);
2135     pos[0]=position(0);
2136     pos[1]=position(1);
2137     pos[2]=position(2);
2138     //bzero((char *)ckovData,sizeof(ckovData)*19);
2139     ckovData[1] = pos[0];                 // X-position for hit
2140     ckovData[2] = pos[1];                 // Y-position for hit
2141     ckovData[3] = pos[2];                 // Z-position for hit
2142     ckovData[6] = 0;                      // dummy track length
2143     //ckovData[11] = gAlice->CurrentTrack();
2144
2145     //printf("\n+++++++++++\nTrack: %d\n++++++++++++\n",gAlice->CurrentTrack());
2146
2147     //AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH");
2148
2149     /********************Store production parameters for Cerenkov photons************************/
2150 //is it a Cerenkov photon?
2151     if (gMC->TrackPid() == 50000050) {
2152
2153       //if (gMC->VolId("GAP ")==gMC->CurrentVolID(copy))
2154         //{
2155           Float_t ckovEnergy = current->Energy();
2156           //energy interval for tracking
2157           if  (ckovEnergy > 5.6e-09 && ckovEnergy < 7.8e-09 )
2158             //if (ckovEnergy > 0)
2159             {
2160               if (gMC->IsTrackEntering()){        //is track entering?
2161                 //printf("Track entered (1)\n");
2162                 if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2163                   {                                                          //is it in freo?
2164                     if (gMC->IsNewTrack()){                          //is it the first step?
2165                       //printf("I'm in!\n");
2166                       Int_t mother = current->GetFirstMother();
2167
2168                       //printf("Second Mother:%d\n",current->GetSecondMother());
2169
2170                       ckovData[10] = mother;
2171                       ckovData[11] = gAlice->CurrentTrack();
2172                       ckovData[12] = 1;             //Media where photon was produced 1->Freon, 2->Quarz
2173                       //printf("Produced in FREO\n");
2174                       fCkovNumber++;
2175                       fFreonProd=1;
2176                       //printf("Index: %d\n",fCkovNumber);
2177                     }    //first step question
2178                   }        //freo question
2179
2180                 if (gMC->IsNewTrack()){                                  //is it first step?
2181                   if (gMC->VolId("QUAR")==gMC->CurrentVolID(copy))             //is it in quarz?
2182                     {
2183                       ckovData[12] = 2;
2184                       //printf("Produced in QUAR\n");
2185                     }    //quarz question
2186                 }        //first step question
2187
2188                 //printf("Before %d\n",fFreonProd);
2189               }   //track entering question
2190
2191               if (ckovData[12] == 1)                                        //was it produced in Freon?
2192                 //if (fFreonProd == 1)
2193                 {
2194                   if (gMC->IsTrackEntering()){                                     //is track entering?
2195                     //printf("Track entered (2)\n");
2196                     //printf("Current volume (should be META): %s\n",gMC->CurrentVolName());
2197                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("META"),gMC->CurrentVolID(copy));
2198                     if (gMC->VolId("META")==gMC->CurrentVolID(copy))                //is it in gap?
2199                       {
2200                         //printf("Got in META\n");
2201                         gMC->TrackMomentum(momentum);
2202                         mom[0]=momentum(0);
2203                         mom[1]=momentum(1);
2204                         mom[2]=momentum(2);
2205                         mom[3]=momentum(3);
2206
2207                         gMC->Gmtod(mom,localMom,2);
2208                         Float_t cophi = TMath::Cos(TMath::ATan2(localMom[0], localMom[1]));
2209                         Float_t t = (1. - .025 / cophi) * (1. - .05 /  cophi);
2210                         /**************** Photons lost in second grid have to be calculated by hand************/
2211                         gMC->GetRandom()->RndmArray(1,ranf);
2212                         if (ranf[0] > t) {
2213                           gMC->StopTrack();
2214                           ckovData[13] = 5;
2215                           AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2216                           //printf("Added One (1)!\n");
2217                           //printf("Lost one in grid\n");
2218                         }
2219                         /**********************************************************************************/
2220                       }    //gap
2221
2222                     //printf("Current volume (should be CSI) (1): %s\n",gMC->CurrentVolName());
2223                     //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2224                     if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))             //is it in csi?
2225                       {
2226                         //printf("Got in CSI\n");
2227                         gMC->TrackMomentum(momentum);
2228                         mom[0]=momentum(0);
2229                         mom[1]=momentum(1);
2230                         mom[2]=momentum(2);
2231                         mom[3]=momentum(3);
2232
2233                         gMC->Gmtod(mom,localMom,2);
2234                         /********* Photons lost by Fresnel reflection have to be calculated by hand********/
2235                         /***********************Cerenkov phtons (always polarised)*************************/
2236                         Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2237                         Double_t localRt = TMath::Sqrt(localTc);
2238                         localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])));
2239                         Double_t cotheta = TMath::Abs(cos(localTheta));
2240                         Float_t t = Fresnel(ckovEnergy*1e9,cotheta,1);
2241                             gMC->GetRandom()->RndmArray(1,ranf);
2242                             if (ranf[0] < t) {
2243                               gMC->StopTrack();
2244                               ckovData[13] = 6;
2245                               AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2246
2247                               //printf("Added One (2)!\n");
2248                               //printf("Lost by Fresnel\n");
2249                             }
2250                             /**********************************************************************************/
2251                       }
2252                   } //track entering?
2253
2254
2255                   /********************Evaluation of losses************************/
2256                   /******************still in the old fashion**********************/
2257
2258                   TArrayI procs;
2259                   Int_t i1 = gMC->StepProcesses(procs);            //number of physics mechanisms acting on the particle
2260                   for (Int_t i = 0; i < i1; ++i) {
2261                     //        Reflection loss
2262                     if (procs[i] == kPLightReflection) {        //was it reflected
2263                       ckovData[13]=10;
2264                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2265                         ckovData[13]=1;
2266                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2267                         ckovData[13]=2;
2268                       //gMC->StopTrack();
2269                       //AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2270                     } //reflection question
2271
2272                     //        Absorption loss
2273                     else if (procs[i] == kPLightAbsorption) {              //was it absorbed?
2274                       //printf("Got in absorption\n");
2275                       ckovData[13]=20;
2276                       if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2277                         ckovData[13]=11;
2278                       if (gMC->CurrentVolID(copy) == gMC->VolId("QUAR"))
2279                         ckovData[13]=12;
2280                       if (gMC->CurrentVolID(copy) == gMC->VolId("META"))
2281                         ckovData[13]=13;
2282                       if (gMC->CurrentVolID(copy) == gMC->VolId("GAP "))
2283                         ckovData[13]=13;
2284
2285                       if (gMC->CurrentVolID(copy) == gMC->VolId("SRIC"))
2286                         ckovData[13]=15;
2287
2288                       //        CsI inefficiency
2289                       if (gMC->CurrentVolID(copy) == gMC->VolId("CSI ")) {
2290                         ckovData[13]=16;
2291                       }
2292                       gMC->StopTrack();
2293                       AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2294                       //printf("Added One (3)!\n");
2295                       //printf("Added cerenkov %d\n",fCkovNumber);
2296                     } //absorption question
2297
2298
2299                     //        Photon goes out of tracking scope
2300                     else if (procs[i] == kPStop) {                 //is it below energy treshold?
2301                       ckovData[13]=21;
2302                       gMC->StopTrack();
2303                       AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2304                       //printf("Added One (4)!\n");
2305                     }   // energy treshold question
2306                   }  //number of mechanisms cycle
2307                   /**********************End of evaluation************************/
2308                 } //freon production question
2309             } //energy interval question
2310         //}//inside the proximity gap question
2311     } //cerenkov photon question
2312
2313     /**************************************End of Production Parameters Storing*********************/
2314
2315
2316     /*******************************Treat photons that hit the CsI (Ckovs and Feedbacks)************/
2317
2318     if (gMC->TrackPid() == 50000050 || gMC->TrackPid() == 50000051) {
2319       //printf("Cerenkov\n");
2320
2321       //if (gMC->TrackPid() == 50000051)
2322         //printf("Tracking a feedback\n");
2323
2324       if (gMC->VolId("CSI ")==gMC->CurrentVolID(copy))
2325         {
2326           //printf("Current volume (should be CSI) (2): %s\n",gMC->CurrentVolName());
2327           //printf("VolId: %d, CurrentVolID: %d\n",gMC->VolId("CSI "),gMC->CurrentVolID(copy));
2328           //printf("Got in CSI\n");
2329           //printf("Tracking a %d\n",gMC->TrackPid());
2330           if (gMC->Edep() > 0.){
2331                 gMC->TrackPosition(position);
2332                 gMC->TrackMomentum(momentum);
2333                 pos[0]=position(0);
2334                 pos[1]=position(1);
2335                 pos[2]=position(2);
2336                 mom[0]=momentum(0);
2337                 mom[1]=momentum(1);
2338                 mom[2]=momentum(2);
2339                 mom[3]=momentum(3);
2340                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2341                 Double_t rt = TMath::Sqrt(tc);
2342                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2343                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2344
2345                 gMC->CurrentVolOffID(2,copy);
2346                 vol[0]=copy;
2347                 idvol=vol[0]-1;
2348
2349
2350                 gMC->Gmtod(pos,localPos,1);
2351
2352                 //Chamber(idvol).GlobaltoLocal(pos,localPos);
2353
2354                 gMC->Gmtod(mom,localMom,2);
2355
2356                 //Chamber(idvol).GlobaltoLocal(mom,localMom);
2357
2358                 gMC->CurrentVolOffID(2,copy);
2359                 vol[0]=copy;
2360                 idvol=vol[0]-1;
2361
2362                 //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2363                         //->Sector(localPos[0], localPos[2]);
2364                 //printf("Sector:%d\n",sector);
2365
2366                 /*if (gMC->TrackPid() == 50000051){
2367                   fFeedbacks++;
2368                   printf("Feedbacks:%d\n",fFeedbacks);
2369                 }*/
2370
2371         //PH            ((AliRICHChamber*) (*fChambers)[idvol])
2372                 ((AliRICHChamber*)fChambers->At(idvol))
2373                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2374                 if(idvol<kNCH) {
2375                     ckovData[0] = gMC->TrackPid();        // particle type
2376                     ckovData[1] = pos[0];                 // X-position for hit
2377                     ckovData[2] = pos[1];                 // Y-position for hit
2378                     ckovData[3] = pos[2];                 // Z-position for hit
2379                     ckovData[4] = theta;                      // theta angle of incidence
2380                     ckovData[5] = phi;                      // phi angle of incidence
2381                     ckovData[8] = (Float_t) fNSDigits;      // first sdigit
2382                     ckovData[9] = -1;                       // last pad hit
2383                     ckovData[13] = 4;                       // photon was detected
2384                     ckovData[14] = mom[0];
2385                     ckovData[15] = mom[1];
2386                     ckovData[16] = mom[2];
2387
2388                     destep = gMC->Edep();
2389                     gMC->SetMaxStep(kBig);
2390                     cherenkovLoss  += destep;
2391                     ckovData[7]=cherenkovLoss;
2392
2393                     nPads = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);
2394
2395                     if (fNSDigits > (Int_t)ckovData[8]) {
2396                         ckovData[8]= ckovData[8]+1;
2397                         ckovData[9]= (Float_t) fNSDigits;
2398                     }
2399
2400                     //printf("Cerenkov loss: %f\n", cherenkovLoss);
2401
2402                     ckovData[17] = nPads;
2403                     //printf("nPads:%d",nPads);
2404
2405                     //TClonesArray *Hits = RICH->Hits();
2406                     AliRICHHit *mipHit =  (AliRICHHit*) (fHits->UncheckedAt(0));
2407                     if (mipHit)
2408                       {
2409                         mom[0] = current->Px();
2410                         mom[1] = current->Py();
2411                         mom[2] = current->Pz();
2412                         Float_t mipPx = mipHit->MomX();
2413                         Float_t mipPy = mipHit->MomY();
2414                         Float_t mipPz = mipHit->MomZ();
2415
2416                         Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
2417                         Float_t rt = TMath::Sqrt(r);
2418                         Float_t mipR = mipPx*mipPx + mipPy*mipPy + mipPz*mipPz;
2419                         Float_t mipRt = TMath::Sqrt(mipR);
2420                         if ((rt*mipRt) > 0)
2421                           {
2422                             coscerenkov = (mom[0]*mipPx + mom[1]*mipPy + mom[2]*mipPz)/(rt*mipRt);
2423                           }
2424                         else
2425                           {
2426                             coscerenkov = 0;
2427                           }
2428                         Float_t cherenkov = TMath::ACos(coscerenkov);
2429                         ckovData[18]=cherenkov;
2430                       }
2431                     //if (sector != -1)
2432                     //{
2433                     AddHit(gAlice->CurrentTrack(),vol,ckovData);
2434                     AddCerenkov(gAlice->CurrentTrack(),vol,ckovData);
2435                     //printf("Added One (5)!\n");
2436                     //}
2437                 }
2438             }
2439         }
2440     }
2441
2442     /***********************************************End of photon hits*********************************************/
2443
2444
2445     /**********************************************Charged particles treatment*************************************/
2446
2447     else if (gMC->TrackCharge())
2448     //else if (1 == 1)
2449       {
2450 //If MIP
2451         /*if (gMC->IsTrackEntering())
2452           {
2453             hits[13]=20;//is track entering?
2454           }*/
2455         if (gMC->VolId("FRE1")==gMC->CurrentVolID(copy) || gMC->VolId("FRE2")==gMC->CurrentVolID(copy))
2456           {
2457             gMC->TrackMomentum(momentum);
2458             mom[0]=momentum(0);
2459             mom[1]=momentum(1);
2460             mom[2]=momentum(2);
2461             mom[3]=momentum(3);
2462             hits [19] = mom[0];
2463             hits [20] = mom[1];
2464             hits [21] = mom[2];
2465             fFreonProd=1;
2466           }
2467
2468         if (gMC->VolId("GAP ")== gMC->CurrentVolID(copy)) {
2469 // Get current particle id (ipart), track position (pos)  and momentum (mom)
2470
2471             gMC->CurrentVolOffID(3,copy);
2472             vol[0]=copy;
2473             idvol=vol[0]-1;
2474
2475             //Int_t sector=((AliRICHChamber*) (*fChambers)[idvol])
2476                         //->Sector(localPos[0], localPos[2]);
2477             //printf("Sector:%d\n",sector);
2478
2479             gMC->TrackPosition(position);
2480             gMC->TrackMomentum(momentum);
2481             pos[0]=position(0);
2482             pos[1]=position(1);
2483             pos[2]=position(2);
2484             mom[0]=momentum(0);
2485             mom[1]=momentum(1);
2486             mom[2]=momentum(2);
2487             mom[3]=momentum(3);
2488
2489             gMC->Gmtod(pos,localPos,1);
2490
2491             //Chamber(idvol).GlobaltoLocal(pos,localPos);
2492
2493             gMC->Gmtod(mom,localMom,2);
2494
2495             //Chamber(idvol).GlobaltoLocal(mom,localMom);
2496
2497             ipart  = gMC->TrackPid();
2498             //
2499             // momentum loss and steplength in last step
2500             destep = gMC->Edep();
2501             step   = gMC->TrackStep();
2502
2503             //
2504             // record hits when track enters ...
2505             if( gMC->IsTrackEntering()) {
2506 //              gMC->SetMaxStep(fMaxStepGas);
2507                 Double_t tc = mom[0]*mom[0]+mom[1]*mom[1];
2508                 Double_t rt = TMath::Sqrt(tc);
2509                 theta   = Float_t(TMath::ATan2(rt,Double_t(mom[2])))*kRaddeg;
2510                 phi     = Float_t(TMath::ATan2(Double_t(mom[1]),Double_t(mom[0])))*kRaddeg;
2511
2512
2513                 Double_t localTc = localMom[0]*localMom[0]+localMom[2]*localMom[2];
2514                 Double_t localRt = TMath::Sqrt(localTc);
2515                 localTheta   = Float_t(TMath::ATan2(localRt,Double_t(localMom[1])))*kRaddeg;
2516                 localPhi     = Float_t(TMath::ATan2(Double_t(localMom[2]),Double_t(localMom[0])))*kRaddeg;
2517
2518                 hits[0] = Float_t(ipart);         // particle type
2519                 hits[1] = localPos[0];                 // X-position for hit
2520                 hits[2] = localPos[1];                 // Y-position for hit
2521                 hits[3] = localPos[2];                 // Z-position for hit
2522                 hits[4] = localTheta;                  // theta angle of incidence
2523                 hits[5] = localPhi;                    // phi angle of incidence
2524                 hits[8] = (Float_t) fNSDigits;    // first sdigit
2525                 hits[9] = -1;                     // last pad hit
2526                 hits[13] = fFreonProd;           // did id hit the freon?
2527                 hits[14] = mom[0];
2528                 hits[15] = mom[1];
2529                 hits[16] = mom[2];
2530                 hits[18] = 0;               // dummy cerenkov angle
2531
2532                 tlength = 0;
2533                 eloss   = 0;
2534                 fFreonProd = 0;
2535
2536                 Chamber(idvol).LocaltoGlobal(localPos,hits+1);
2537
2538
2539                 //To make chamber coordinates x-y had to pass localPos[0], localPos[2]
2540                 xhit    = localPos[0];
2541                 yhit    = localPos[2];
2542                 // Only if not trigger chamber
2543                 if(idvol<kNCH) {
2544                     //
2545                     //  Initialize hit position (cursor) in the segmentation model
2546           //PH              ((AliRICHChamber*) (*fChambers)[idvol])
2547                     ((AliRICHChamber*)fChambers->At(idvol))
2548                         ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2549                 }
2550             }
2551
2552             //
2553             // Calculate the charge induced on a pad (disintegration) in case
2554             //
2555             // Mip left chamber ...
2556             if( gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()){
2557                 gMC->SetMaxStep(kBig);
2558                 eloss   += destep;
2559                 tlength += step;
2560
2561
2562                 // Only if not trigger chamber
2563                 if(idvol<kNCH) {
2564                   if (eloss > 0)
2565                     {
2566                       if(gMC->TrackPid() == kNeutron)
2567                         printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2568                       nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2569                       hits[17] = nPads;
2570                       //printf("nPads:%d",nPads);
2571                     }
2572                 }
2573
2574                 hits[6]=tlength;
2575                 hits[7]=eloss;
2576                 if (fNSDigits > (Int_t)hits[8]) {
2577                     hits[8]= hits[8]+1;
2578                     hits[9]= (Float_t) fNSDigits;
2579                 }
2580
2581                 //if(sector !=-1)
2582                 new(lhits[fNhits++]) AliRICHHit(fIshunt,gAlice->CurrentTrack(),vol,hits);
2583                 eloss = 0;
2584                 //
2585                 // Check additional signal generation conditions
2586                 // defined by the segmentation
2587                 // model (boundary crossing conditions)
2588             } else if
2589           //PH          (((AliRICHChamber*) (*fChambers)[idvol])
2590                 (((AliRICHChamber*)fChambers->At(idvol))
2591                  ->SigGenCond(localPos[0], localPos[2], localPos[1]))
2592             {
2593           //PH          ((AliRICHChamber*) (*fChambers)[idvol])
2594                 ((AliRICHChamber*)fChambers->At(idvol))
2595                     ->SigGenInit(localPos[0], localPos[2], localPos[1]);
2596                 if (eloss > 0)
2597                   {
2598                     if(gMC->TrackPid() == kNeutron)
2599                       printf("\n\n\n\n\n Neutron Making Pad Hit!!! \n\n\n\n");
2600                     nPads = Hits2SDigits(xhit,yhit,eloss,idvol,kMip);
2601                     hits[17] = nPads;
2602                     //printf("Npads:%d",NPads);
2603                   }
2604                 xhit     = localPos[0];
2605                 yhit     = localPos[2];
2606                 eloss    = destep;
2607                 tlength += step ;
2608                 //
2609                 // nothing special  happened, add up energy loss
2610             } else {
2611                 eloss   += destep;
2612                 tlength += step ;
2613             }
2614         }
2615       }
2616     /*************************************************End of MIP treatment**************************************/
2617    //}
2618 }//void AliRICH::StepManager()
2619
2620 void AliRICH::FindClusters(Int_t nev,Int_t lastEntry)
2621 {
2622
2623 //
2624 // Loop on chambers and on cathode planes
2625 //
2626     for (Int_t icat=1;icat<2;icat++) {
2627         gAlice->ResetDigits();
2628         gAlice->TreeD()->GetEvent(0);
2629         for (Int_t ich=0;ich<kNCH;ich++) {
2630       //PH        AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
2631           AliRICHChamber* iChamber=(AliRICHChamber*)fChambers->At(ich);
2632           TClonesArray *pRICHdigits  = this->DigitsAddress(ich);
2633           if (pRICHdigits == 0)
2634               continue;
2635           //
2636           // Get ready the current chamber stuff
2637           //
2638           AliRICHResponse* response = iChamber->GetResponseModel();
2639           AliSegmentation*  seg = iChamber->GetSegmentationModel();
2640           AliRICHClusterFinder* rec = iChamber->GetReconstructionModel();
2641           if (seg) {
2642               rec->SetSegmentation(seg);
2643               rec->SetResponse(response);
2644               rec->SetDigits(pRICHdigits);
2645               rec->SetChamber(ich);
2646               if (nev==0) rec->CalibrateCOG();
2647               rec->FindRawClusters();
2648           }
2649           TClonesArray *fRch;
2650           fRch=RawClustAddress(ich);
2651           fRch->Sort();
2652         } // for ich
2653
2654         gAlice->TreeR()->Fill();
2655         TClonesArray *fRch;
2656         for (int i=0;i<kNCH;i++) {
2657             fRch=RawClustAddress(i);
2658             int nraw=fRch->GetEntriesFast();
2659             printf ("Chamber %d, raw clusters %d\n",i,nraw);
2660         }
2661
2662         ResetRawClusters();
2663
2664     } // for icat
2665
2666     char hname[30];
2667     sprintf(hname,"TreeR%d",nev);
2668     gAlice->TreeR()->Write(hname,kOverwrite,0);
2669     gAlice->TreeR()->Reset();
2670
2671     //gObjectTable->Print();
2672 }
2673
2674 AliRICHSDigit* AliRICH::FirstPad(AliRICHHit*  hit,TClonesArray *clusters )
2675 {
2676 //
2677     // Initialise the pad iterator
2678     // Return the address of the first sdigit for hit
2679     TClonesArray *theClusters = clusters;
2680     Int_t nclust = theClusters->GetEntriesFast();
2681     if (nclust && hit->PHlast() > 0) {
2682         sMaxIterPad=Int_t(hit->PHlast());
2683         sCurIterPad=Int_t(hit->PHfirst());
2684         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2685     } else {
2686         return 0;
2687     }
2688
2689 }
2690
2691 AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters)
2692 {
2693
2694   // Iterates over pads
2695
2696     sCurIterPad++;
2697     if (sCurIterPad <= sMaxIterPad) {
2698         return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
2699     } else {
2700         return 0;
2701     }
2702 }
2703
2704 AliRICH& AliRICH::operator=(const AliRICH& rhs)
2705 {
2706 // Assignment operator
2707     return *this;
2708
2709 }
2710
2711 void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
2712 {
2713
2714   Int_t NpadX = 162;                 // number of pads on X
2715   Int_t NpadY = 162;                 // number of pads on Y
2716
2717   Int_t Pad[162][162];
2718   for (Int_t i=0;i<NpadX;i++) {
2719     for (Int_t j=0;j<NpadY;j++) {
2720       Pad[i][j]=0;
2721     }
2722   }
2723
2724   //  Create some histograms
2725
2726   TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
2727   TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
2728   TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
2729   TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
2730   TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
2731   TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
2732   TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
2733   TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
2734   TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
2735   TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
2736   TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
2737   TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
2738   TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
2739   TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
2740   TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
2741   TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
2742   TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2743   TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
2744   TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
2745   TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2746   TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
2747   TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
2748   TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
2749   TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
2750   TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
2751   //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
2752   TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
2753   TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
2754   TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
2755   TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
2756   TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
2757   TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
2758
2759
2760
2761
2762 //   Start loop over events
2763
2764   Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
2765   Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
2766   Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
2767   TRandom* random=0;
2768
2769    for (int nev=0; nev<= evNumber2; nev++) {
2770        Int_t nparticles = gAlice->GetEvent(nev);
2771
2772
2773        printf ("Event number       : %d\n",nev);
2774        printf ("Number of particles: %d\n",nparticles);
2775        if (nev < evNumber1) continue;
2776        if (nparticles <= 0) return;
2777
2778 // Get pointers to RICH detector and Hits containers
2779
2780        AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
2781
2782        TTree *treeH = gAlice->TreeH();
2783        Int_t ntracks =(Int_t) treeH->GetEntries();
2784
2785 // Start loop on tracks in the hits containers
2786
2787        for (Int_t track=0; track<ntracks;track++) {
2788            printf ("Processing Track: %d\n",track);
2789            gAlice->ResetHits();
2790            treeH->GetEvent(track);
2791
2792            for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
2793                mHit;
2794                mHit=(AliRICHHit*)pRICH->NextHit())
2795              {
2796                //Int_t nch  = mHit->fChamber;              // chamber number
2797                //Float_t x  = mHit->X();                    // x-pos of hit
2798                //Float_t y  = mHit->Z();                    // y-pos
2799                //Float_t z  = mHit->Y();
2800                //Float_t phi = mHit->Phi();                 //Phi angle of incidence
2801                Float_t theta = mHit->Theta();             //Theta angle of incidence
2802                Float_t px = mHit->MomX();
2803                Float_t py = mHit->MomY();
2804                Int_t index = mHit->Track();
2805                Int_t particle = (Int_t)(mHit->Particle());
2806                Float_t R;
2807                Float_t PTfinal;
2808                Float_t PTvertex;
2809
2810               TParticle *current = gAlice->Particle(index);
2811
2812               //Float_t energy=current->Energy();
2813
2814               R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
2815               PTfinal=TMath::Sqrt(px*px + py*py);
2816               PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
2817
2818
2819
2820               if (TMath::Abs(particle) < 10000000)
2821                 {
2822                   hitsTheta->Fill(theta,(float) 1);
2823                   if (R<5)
2824                     {
2825                       if (PTvertex>.5 && PTvertex<=1)
2826                         {
2827                           hitsTheta500MeV->Fill(theta,(float) 1);
2828                         }
2829                       if (PTvertex>1 && PTvertex<=2)
2830                         {
2831                           hitsTheta1GeV->Fill(theta,(float) 1);
2832                         }
2833                       if (PTvertex>2 && PTvertex<=3)
2834                         {
2835                           hitsTheta2GeV->Fill(theta,(float) 1);
2836                         }
2837                       if (PTvertex>3)
2838                         {
2839                           hitsTheta3GeV->Fill(theta,(float) 1);
2840                         }
2841                     }
2842
2843                 }
2844
2845               //if (nch == 3)
2846                 //{
2847
2848               //printf("Particle type: %d\n",current->GetPdgCode());
2849               if (TMath::Abs(particle) < 50000051)
2850                 {
2851                   //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
2852                   if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
2853                     {
2854                       //gMC->Rndm(&random, 1);
2855                       if (random->Rndm() < .1)
2856                         production->Fill(current->Vz(),R,(float) 1);
2857                       if (TMath::Abs(particle) == 50000050)
2858                         //if (TMath::Abs(particle) > 50000000)
2859                         {
2860                           photons +=1;
2861                           if (R<5)
2862                             {
2863                               primaryphotons +=1;
2864                               if (current->Energy()>0.001)
2865                                 highprimaryphotons +=1;
2866                             }
2867                         }
2868                       if (TMath::Abs(particle) == 2112)
2869                         {
2870                           neutron +=1;
2871                           if (current->Energy()>0.0001)
2872                             highneutrons +=1;
2873                         }
2874                     }
2875                   if (TMath::Abs(particle) < 50000000)
2876                     {
2877                       production->Fill(current->Vz(),R,(float) 1);
2878                       //printf("Adding %d at %f\n",particle,R);
2879                     }
2880                   //mip->Fill(x,y,(float) 1);
2881                 }
2882
2883               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2884                 {
2885                   if (R<5)
2886                     {
2887                       pionptspectravertex->Fill(PTvertex,(float) 1);
2888                       pionptspectrafinal->Fill(PTfinal,(float) 1);
2889                     }
2890                 }
2891
2892               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2893                   || TMath::Abs(particle)==311)
2894                 {
2895                   if (R<5)
2896                     {
2897                       kaonptspectravertex->Fill(PTvertex,(float) 1);
2898                       kaonptspectrafinal->Fill(PTfinal,(float) 1);
2899                     }
2900                 }
2901
2902
2903               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
2904                 {
2905                   pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2906                   //printf ("fParticle: %d, PDG code:%d\n",particle,current->GetPdgCode());
2907                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2908                     pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2909                   if (R>250 && R<450)
2910                     {
2911                       pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2912                       //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2913                     }
2914                   //printf("Pion mass: %e\n",current->GetCalcMass());
2915                   pion +=1;
2916                   if (TMath::Abs(particle)==211)
2917                     {
2918                       chargedpions +=1;
2919                       if (R<5)
2920                         {
2921                           primarypions +=1;
2922                           if (current->Energy()>1)
2923                             highprimarypions +=1;
2924                         }
2925                     }
2926                 }
2927               if (TMath::Abs(particle)==2212)
2928                 {
2929                   protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2930                   //ptspectra->Fill(Pt,(float) 1);
2931                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2932                     protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2933                   if (R>250 && R<450)
2934                     protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2935                   //printf("\n\n\n\n\n\n\nProton mass: %e\n\n\n\n\n\n\n\n\n",current->GetCalcMass());
2936                   proton +=1;
2937                 }
2938               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
2939                   || TMath::Abs(particle)==311)
2940                 {
2941                   kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2942                   //ptspectra->Fill(Pt,(float) 1);
2943                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2944                     kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2945                   if (R>250 && R<450)
2946                     kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2947                   //printf("Kaon mass: %e\n",current->GetCalcMass());
2948                   kaon +=1;
2949                   if (TMath::Abs(particle)==321)
2950                     {
2951                       chargedkaons +=1;
2952                       if (R<5)
2953                         {
2954                           primarykaons +=1;
2955                           if (current->Energy()>1)
2956                             highprimarykaons +=1;
2957                         }
2958                     }
2959                 }
2960               if (TMath::Abs(particle)==11)
2961                 {
2962                   electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2963                   //ptspectra->Fill(Pt,(float) 1);
2964                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2965                     electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2966                   if (R>250 && R<450)
2967                     electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2968                   //printf("Electron mass: %e\n",current->GetCalcMass());
2969                   if (particle == 11)
2970                     electron +=1;
2971                   if (particle == -11)
2972                     positron +=1;
2973                 }
2974               if (TMath::Abs(particle)==13)
2975                 {
2976                   muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2977                   //ptspectra->Fill(Pt,(float) 1);
2978                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2979                     muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2980                   if (R>250 && R<450)
2981                     muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2982                   //printf("Muon mass: %e\n",current->GetCalcMass());
2983                   muon +=1;
2984                 }
2985               if (TMath::Abs(particle)==2112)
2986                 {
2987                   neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2988                   //ptspectra->Fill(Pt,(float) 1);
2989                   if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
2990                     neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2991                   if (R>250 && R<450)
2992                     {
2993                       neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
2994                       //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\R:%f\n\n\n\n\n\n\n\n\n",R);
2995                     }
2996                   //printf("Neutron mass: %e\n",current->GetCalcMass());
2997                   neutron +=1;
2998                 }
2999               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3000                 {
3001                   if (current->Energy()-current->GetCalcMass()>1)
3002                     {
3003                       chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3004                       if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
3005                         chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3006                       if (R>250 && R<450)
3007                         chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
3008                     }
3009                 }
3010               //printf("Hits:%d\n",hit);
3011               //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3012               // Fill the histograms
3013               //Nh1+=nhits;
3014               //h->Fill(x,y,(float) 1);
3015               //}
3016               //}
3017            }
3018
3019        }
3020
3021    }
3022    //   }
3023
3024    TStyle *mystyle=new TStyle("Plain","mystyle");
3025    mystyle->SetPalette(1,0);
3026    mystyle->cd();
3027
3028    //Create canvases, set the view range, show histograms
3029
3030     TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
3031     c2->Divide(2,2);
3032     //c2->SetFillColor(42);
3033
3034     c2->cd(1);
3035     hitsTheta500MeV->SetFillColor(5);
3036     hitsTheta500MeV->Draw();
3037     c2->cd(2);
3038     hitsTheta1GeV->SetFillColor(5);
3039     hitsTheta1GeV->Draw();
3040     c2->cd(3);
3041     hitsTheta2GeV->SetFillColor(5);
3042     hitsTheta2GeV->Draw();
3043     c2->cd(4);
3044     hitsTheta3GeV->SetFillColor(5);
3045     hitsTheta3GeV->Draw();
3046
3047
3048
3049     TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
3050     c15->cd();
3051     production->SetFillColor(42);
3052     production->SetXTitle("z (m)");
3053     production->SetYTitle("R (m)");
3054     production->Draw();
3055
3056     TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
3057     c10->Divide(2,2);
3058     c10->cd(1);
3059     pionptspectravertex->SetFillColor(5);
3060     pionptspectravertex->SetXTitle("Pt (GeV)");
3061     pionptspectravertex->Draw();
3062     c10->cd(2);
3063     pionptspectrafinal->SetFillColor(5);
3064     pionptspectrafinal->SetXTitle("Pt (GeV)");
3065     pionptspectrafinal->Draw();
3066     c10->cd(3);
3067     kaonptspectravertex->SetFillColor(5);
3068     kaonptspectravertex->SetXTitle("Pt (GeV)");
3069     kaonptspectravertex->Draw();
3070     c10->cd(4);
3071     kaonptspectrafinal->SetFillColor(5);
3072     kaonptspectrafinal->SetXTitle("Pt (GeV)");
3073     kaonptspectrafinal->Draw();
3074
3075
3076    TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
3077    c16->Divide(2,1);
3078
3079    c16->cd(1);
3080    //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
3081    electronspectra1->SetFillColor(5);
3082    electronspectra1->SetXTitle("log(GeV)");
3083    electronspectra2->SetFillColor(46);
3084    electronspectra2->SetXTitle("log(GeV)");
3085    electronspectra3->SetFillColor(10);
3086    electronspectra3->SetXTitle("log(GeV)");
3087    //c13->SetLogx();
3088    electronspectra1->Draw();
3089    electronspectra2->Draw("same");
3090    electronspectra3->Draw("same");
3091
3092    c16->cd(2);
3093    //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
3094    muonspectra1->SetFillColor(5);
3095    muonspectra1->SetXTitle("log(GeV)");
3096    muonspectra2->SetFillColor(46);
3097    muonspectra2->SetXTitle("log(GeV)");
3098    muonspectra3->SetFillColor(10);
3099    muonspectra3->SetXTitle("log(GeV)");
3100    //c14->SetLogx();
3101    muonspectra1->Draw();
3102    muonspectra2->Draw("same");
3103    muonspectra3->Draw("same");
3104
3105    //c16->cd(3);
3106    //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
3107    //neutronspectra1->SetFillColor(42);
3108    //neutronspectra1->SetXTitle("log(GeV)");
3109    //neutronspectra2->SetFillColor(46);
3110    //neutronspectra2->SetXTitle("log(GeV)");
3111    //neutronspectra3->SetFillColor(10);
3112    //neutronspectra3->SetXTitle("log(GeV)");
3113    //c16->SetLogx();
3114    //neutronspectra1->Draw();
3115    //neutronspectra2->Draw("same");
3116    //neutronspectra3->Draw("same");
3117
3118    TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
3119    //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
3120    c9->Divide(2,2);
3121
3122    c9->cd(1);
3123    pionspectra1->SetFillColor(5);
3124    pionspectra1->SetXTitle("log(GeV)");
3125    pionspectra2->SetFillColor(46);
3126    pionspectra2->SetXTitle("log(GeV)");
3127    pionspectra3->SetFillColor(10);
3128    pionspectra3->SetXTitle("log(GeV)");
3129    //c9->SetLogx();
3130    pionspectra1->Draw();
3131    pionspectra2->Draw("same");
3132    pionspectra3->Draw("same");
3133
3134    c9->cd(2);
3135    //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
3136    protonspectra1->SetFillColor(5);
3137    protonspectra1->SetXTitle("log(GeV)");
3138    protonspectra2->SetFillColor(46);
3139    protonspectra2->SetXTitle("log(GeV)");
3140    protonspectra3->SetFillColor(10);
3141    protonspectra3->SetXTitle("log(GeV)");
3142    //c10->SetLogx();
3143    protonspectra1->Draw();
3144    protonspectra2->Draw("same");
3145    protonspectra3->Draw("same");
3146
3147    c9->cd(3);
3148    //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700);
3149    kaonspectra1->SetFillColor(5);
3150    kaonspectra1->SetXTitle("log(GeV)");
3151    kaonspectra2->SetFillColor(46);
3152    kaonspectra2->SetXTitle("log(GeV)");
3153    kaonspectra3->SetFillColor(10);
3154    kaonspectra3->SetXTitle("log(GeV)");
3155    //c11->SetLogx();
3156    kaonspectra1->Draw();
3157    kaonspectra2->Draw("same");
3158    kaonspectra3->Draw("same");
3159
3160    c9->cd(4);
3161    //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
3162    chargedspectra1->SetFillColor(5);
3163    chargedspectra1->SetXTitle("log(GeV)");
3164    chargedspectra2->SetFillColor(46);
3165    chargedspectra2->SetXTitle("log(GeV)");
3166    chargedspectra3->SetFillColor(10);
3167    chargedspectra3->SetXTitle("log(GeV)");
3168    //c12->SetLogx();
3169    chargedspectra1->Draw();
3170    chargedspectra2->Draw("same");
3171    chargedspectra3->Draw("same");
3172
3173
3174
3175    printf("*****************************************\n");
3176    printf("* Particle                   *  Counts  *\n");
3177    printf("*****************************************\n");
3178
3179    printf("* Pions:                     *   %4d   *\n",pion);
3180    printf("* Charged Pions:             *   %4d   *\n",chargedpions);
3181    printf("* Primary Pions:             *   %4d   *\n",primarypions);
3182    printf("* Primary Pions (p>1GeV/c):  *   %4d   *\n",highprimarypions);
3183    printf("* Kaons:                     *   %4d   *\n",kaon);
3184    printf("* Charged Kaons:             *   %4d   *\n",chargedkaons);
3185    printf("* Primary Kaons:             *   %4d   *\n",primarykaons);
3186    printf("* Primary Kaons (p>1GeV/c):  *   %4d   *\n",highprimarykaons);
3187    printf("* Muons:                     *   %4d   *\n",muon);
3188    printf("* Electrons:                 *   %4d   *\n",electron);
3189    printf("* Positrons:                 *   %4d   *\n",positron);
3190    printf("* Protons:                   *   %4d   *\n",proton);
3191    printf("* All Charged:               *   %4d   *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
3192    printf("*****************************************\n");
3193    //printf("* Photons:                   *   %3.1f   *\n",photons);
3194    //printf("* Primary Photons:           *   %3.1f   *\n",primaryphotons);
3195    //printf("* Primary Photons (p>1MeV/c):*   %3.1f   *\n",highprimaryphotons);
3196    //printf("*****************************************\n");
3197    //printf("* Neutrons:                  *   %3.1f   *\n",neutron);
3198    //printf("* Neutrons (p>100keV/c):     *   %3.1f   *\n",highneutrons);
3199    //printf("*****************************************\n");
3200
3201    if (gAlice->TreeD())
3202      {
3203        gAlice->TreeD()->GetEvent(0);
3204
3205        Float_t occ[7];
3206        Float_t sum=0;
3207        Float_t mean=0;
3208        printf("\n*****************************************\n");
3209        printf("* Chamber   * Digits      * Occupancy   *\n");
3210        printf("*****************************************\n");
3211
3212        for (Int_t ich=0;ich<7;ich++)
3213          {
3214            TClonesArray *Digits = DigitsAddress(ich);    //  Raw clusters branch
3215            Int_t ndigits = Digits->GetEntriesFast();
3216            occ[ich] = Float_t(ndigits)/(160*144);
3217            sum += Float_t(ndigits)/(160*144);
3218            printf("*   %d      *    %d      *   %3.1f%%     *\n",ich,ndigits,occ[ich]*100);
3219          }
3220        mean = sum/7;
3221        printf("*****************************************\n");
3222        printf("* Mean occupancy          *   %3.1f%%     *\n",mean*100);
3223        printf("*****************************************\n");
3224      }
3225
3226   printf("\nEnd of analysis\n");
3227
3228 }
3229
3230 //_________________________________________________________________________________________________
3231
3232
3233 void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
3234 {
3235
3236 AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
3237    AliRICHSegmentationV0*  segmentation;
3238    AliRICHChamber*       chamber;
3239
3240    chamber = &(pRICH->Chamber(0));
3241    segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel();
3242
3243    Int_t NpadX = segmentation->Npx();                 // number of pads on X
3244    Int_t NpadY = segmentation->Npy();                 // number of pads on Y
3245
3246    //Int_t Pad[144][160];
3247    /*for (Int_t i=0;i<NpadX;i++) {
3248      for (Int_t j=0;j<NpadY;j++) {
3249        Pad[i][j]=0;
3250      }
3251    } */
3252
3253
3254    Int_t xmin= -NpadX/2;
3255    Int_t xmax=  NpadX/2;
3256    Int_t ymin= -NpadY/2;
3257    Int_t ymax=  NpadY/2;
3258
3259    Float_t PTfinal = 0;
3260    Int_t pionCount = 0;
3261    Int_t kaonCount = 0;
3262    Int_t protonCount = 0;
3263
3264    TH2F *feedback = 0;
3265    TH2F *mip = 0;
3266    TH2F *cerenkov = 0;
3267    TH2F *h = 0;
3268    TH1F *hitsX = 0;
3269    TH1F *hitsY = 0;
3270
3271    TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-25,25,150,-45,5);
3272
3273    if (diaglevel == 1)
3274      {
3275        printf("Single Ring Hits\n");
3276        feedback = new TH2F("feedback","Feedback hit distribution",150,-20,20,150,-35,5);
3277        mip = new TH2F("mip","Mip hit distribution",150,-20,20,150,-35,5);
3278        cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-20,20,150,-35,5);
3279        h = new TH2F("h","Detector hit distribution",150,-20,20,150,-35,5);
3280        hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-50,50);
3281        hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,50);
3282      }
3283    else
3284      {
3285        printf("Full Event Hits\n");
3286
3287        feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
3288        mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
3289        cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
3290        h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300);
3291        hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
3292        hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
3293      }
3294
3295
3296
3297    TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3298    TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3299    TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3300    TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3301    TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3302    TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3303    TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
3304
3305    TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
3306    TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",100,.35,.8);
3307    TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
3308    TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
3309    TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
3310    TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
3311    TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
3312    TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
3313    TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
3314    //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
3315    TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
3316    TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
3317    TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
3318    TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
3319    TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
3320    TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
3321    TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
3322    TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
3323    TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
3324    TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
3325    TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
3326    TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",50,0,360);
3327    TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of theta angle of incidence",50,0,15);
3328    TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",50,.5,1);
3329    TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",100,0,15);
3330    TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",100,0,360);
3331    TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",100,.35,.8);
3332    TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",100,.35,.8);
3333    TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
3334    TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
3335    TH2F *identification = new TH2F("identification","Particle Identification",100,1,5,100,0,.8);
3336    TH1F *OriginalOmega = new TH1F("Original Omega","Cerenkov angle per track",100,.35,.8);
3337    TH1F *OriginalPhi = new TH1F("Original Phi","Distribution of phi angle of incidence per track",100,0,360);
3338    TH1F *OriginalTheta = new TH1F("Original Theta","Distribution of theta angle per track",100,0,15);
3339    TH1F *OmegaError = new TH1F("Omega Error","Difference between original an reconstructed cerenkov angle",100,0,.2);
3340    TH1F *PhiError = new TH1F("Phi Error","Difference between original an reconstructed phi angle",100,0,360);
3341    TH1F *ThetaError = new TH1F("Theta Error","Difference between original an reconstructed phi angle",100,0,15);
3342
3343
3344 //   Start loop over events
3345
3346    Int_t Nh=0;
3347    Int_t pads=0;
3348    Int_t Nh1=0;
3349    Int_t mothers[80000];
3350    Int_t mothers2[80000];
3351    Float_t mom[3];
3352    Int_t nraw=0;
3353    Int_t phot=0;
3354    Int_t feed=0;
3355    Int_t padmip=0;
3356    Float_t x=0,y=0;
3357
3358    Float_t chiSquareOmega = 0;
3359    Float_t chiSquareTheta = 0;
3360    Float_t chiSquarePhi = 0;
3361
3362    Float_t recEffEvent = 0;
3363    Float_t recEffTotal = 0;
3364
3365    Float_t trackglob[3];
3366    Float_t trackloc[3];
3367
3368
3369    for (Int_t i=0;i<100;i++) mothers[i]=0;
3370
3371    for (int nev=0; nev<= evNumber2; nev++) {
3372        Int_t nparticles = gAlice->GetEvent(nev);
3373
3374
3375        //cout<<"nev  "<<nev<<endl;
3376        printf ("\n**********************************\nProcessing Event: %d\n",nev);
3377        //cout<<"nparticles  "<<nparticles<<endl;
3378        printf ("Particles       : %d\n\n",nparticles);
3379        if (nev < evNumber1) continue;
3380        if (nparticles <= 0) return;
3381
3382 // Get pointers to RICH detector and Hits containers
3383
3384
3385        TTree *TH = gAlice->TreeH();
3386        Stat_t ntracks = TH->GetEntries();
3387
3388        // Start loop on tracks in the hits containers
3389        //Int_t Nc=0;
3390        for (Int_t track=0; track<ntracks;track++) {
3391
3392          printf ("\nProcessing Track: %d\n",track);
3393          gAlice->ResetHits();
3394          TH->GetEvent(track);
3395          Int_t nhits = pRICH->Hits()->GetEntriesFast();
3396          if (nhits) Nh+=nhits;
3397          printf("Hits            : %d\n",nhits);
3398          for(AliRICHHit* mHit=(AliRICHHit*)pRICH->FirstHit(-1);
3399              mHit;
3400              mHit=(AliRICHHit*)pRICH->NextHit())
3401            {
3402              Int_t nch  = mHit->Chamber();              // chamber number
3403              trackglob[0] = mHit->X();                 // x-pos of hit
3404              trackglob[1] = mHit->Y();
3405              trackglob[2] = mHit->Z();                 // y-pos of hit
3406              //x  = mHit->X();                           // x-pos of hit
3407              //y  = mHit->Z();                           // y-pos
3408              Float_t phi = mHit->Phi();                 //Phi angle of incidence
3409              Float_t theta = mHit->Theta();             //Theta angle of incidence
3410              Int_t index = mHit->Track();
3411              Int_t particle = (Int_t)(mHit->Particle());
3412              //Int_t freon = (Int_t)(mHit->fLoss);
3413              Float_t px = mHit->MomX();
3414              Float_t py = mHit->MomY();
3415
3416              if (TMath::Abs(particle) < 10000000)
3417                {
3418                  PTfinal=TMath::Sqrt(px*px + py*py);
3419                  //printf("PTfinal 0: %f\n",PTfinal);
3420                }
3421
3422              chamber = &(pRICH->Chamber(nch-1));
3423
3424              //printf("Nch:%d\n",nch);
3425
3426              chamber->GlobaltoLocal(trackglob,trackloc);
3427
3428              chamber->LocaltoGlobal(trackloc,trackglob);
3429
3430
3431              x=trackloc[0];
3432              y=trackloc[2];
3433
3434              hitsX->Fill(x,(float) 1);
3435              hitsY->Fill(y,(float) 1);
3436
3437               //printf("Particle:%9d\n",particle);
3438
3439               TParticle *current = (TParticle*)gAlice->Particle(index);
3440               //printf("Particle type: %d\n",sizeoff(Particles));
3441
3442               hitsTheta->Fill(theta,(float) 1);
3443               //hitsPhi->Fill(phi,(float) 1);
3444               //if (pRICH->GetDebugLevel() == -1)
3445               //printf("Theta:%f, Phi:%f\n",theta,phi);
3446
3447               //printf("Debug Level:%d\n",pRICH->GetDebugLevel());
3448
3449               if (current->GetPdgCode() < 10000000)
3450                 {
3451                   mip->Fill(x,y,(float) 1);
3452                   //printf("adding mip\n");
3453                   //if (current->Energy() - current->GetCalcMass()>1 && freon==1)
3454                   //{
3455                   hitsPhi->Fill(TMath::Abs(phi),(float) 1);
3456                   //hitsTheta->Fill(theta,(float) 1);
3457                   //printf("Theta:%f, Phi:%f\n",theta,phi);
3458                   //}
3459                 }
3460
3461               if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
3462                 {
3463                   pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3464                 }
3465               if (TMath::Abs(particle)==2212)
3466                 {
3467                   protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3468                 }
3469               if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310
3470                   || TMath::Abs(particle)==311)
3471                 {
3472                   kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3473                 }
3474               if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
3475                 {
3476                   if (current->Energy() - current->GetCalcMass()>1)
3477                     chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
3478                 }
3479               //printf("Hits:%d\n",hit);
3480               //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
3481               // Fill the histograms
3482               Nh1+=nhits;
3483               h->Fill(x,y,(float) 1);
3484                   //}
3485               //}
3486            }
3487
3488            Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
3489            //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
3490            //totalphotonsevent->Fill(ncerenkovs,(float) 1);
3491
3492            if (ncerenkovs) {
3493              printf("Cerenkovs       : %d\n",ncerenkovs);
3494              totalphotonsevent->Fill(ncerenkovs,(float) 1);
3495              for (Int_t hit=0;hit<ncerenkovs;hit++) {
3496                AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
3497                Int_t nchamber = cHit->fChamber;     // chamber number
3498                Int_t index =    cHit->Track();
3499                //Int_t pindex =   (Int_t)(cHit->fIndex);
3500                trackglob[0] = cHit->X();                 // x-pos of hit
3501                trackglob[1] = cHit->Y();
3502                trackglob[2] = cHit->Z();                 // y-pos of hit
3503                //Float_t cx  =      cHit->X();                // x-position
3504                //Float_t cy  =      cHit->Z();                // y-position
3505                Int_t cmother =  cHit->fCMother;      // Index of mother particle
3506                Int_t closs =    (Int_t)(cHit->fLoss);           // How did the particle get lost?
3507                Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
3508
3509                chamber = &(pRICH->Chamber(nchamber-1));
3510
3511                //printf("Nch:%d\n",nch);
3512
3513                chamber->GlobaltoLocal(trackglob,trackloc);
3514
3515                chamber->LocaltoGlobal(trackloc,trackglob);
3516
3517
3518                Float_t cx=trackloc[0];
3519                Float_t cy=trackloc[2];
3520
3521                //printf ("Cerenkov hit number %d/%d, X:%f, Y:%f\n",hit,ncerenkovs,cx,cy);
3522
3523
3524                //printf("Particle:%9d\n",index);
3525
3526                TParticle *current = (TParticle*)gAlice->Particle(index);
3527                Float_t energyckov = current->Energy();
3528
3529                if (current->GetPdgCode() == 50000051)
3530                  {
3531                    if (closs==4)
3532                      {
3533                        feedback->Fill(cx,cy,(float) 1);
3534                        feed++;
3535                      }
3536                  }
3537                if (current->GetPdgCode() == 50000050)
3538                  {
3539
3540                    if (closs !=4)
3541                      {
3542                        phspectra2->Fill(energyckov*1e9,(float) 1);
3543                      }
3544
3545                    if (closs==4)
3546                      {
3547                        cerenkov->Fill(cx,cy,(float) 1);
3548
3549                        //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy);
3550
3551                        //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
3552                        AliRICHHit* mipHit = (AliRICHHit*) pRICH->Hits()->UncheckedAt(0);
3553                        mom[0] = current->Px();
3554                        mom[1] = current->Py();
3555                        mom[2] = current->Pz();
3556                        //mom[0] = cHit->fMomX;
3557                        // mom[1] = cHit->fMomZ;
3558                        //mom[2] = cHit->fMomY;
3559                        //Float_t energymip = MIP->Energy();
3560                        //Float_t Mip_px = mipHit->fMomFreoX;
3561                        //Float_t Mip_py = mipHit->fMomFreoY;
3562                        //Float_t Mip_pz = mipHit->fMomFreoZ;
3563                        //Float_t Mip_px = MIP->Px();
3564                        //Float_t Mip_py = MIP->Py();
3565                        //Float_t Mip_pz = MIP->Pz();
3566
3567
3568
3569                        //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
3570                        //Float_t rt = TMath::Sqrt(r);
3571                        //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz;
3572                        //Float_t Mip_rt = TMath::Sqrt(Mip_r);
3573                        //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
3574                        //Float_t cherenkov = TMath::ACos(coscerenkov);
3575                        ckovangle->Fill(cherenkov,(float) 1);                           //Cerenkov angle calculus
3576                        //printf("Cherenkov: %f\n",cherenkov);
3577                        Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
3578                        hckphi->Fill(ckphi,(float) 1);
3579
3580
3581                        //Float_t mix = MIP->Vx();
3582                        //Float_t miy = MIP->Vy();
3583                        Float_t mx = mipHit->X();
3584                        Float_t my = mipHit->Z();
3585                        //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
3586                        Float_t dx = trackglob[0] - mx;
3587                        Float_t dy = trackglob[2] - my;
3588                        //printf("Dx:%f, Dy:%f\n",dx,dy);
3589                        Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
3590                        //printf("Final radius:%f\n",final_radius);
3591                        radius->Fill(final_radius,(float) 1);
3592
3593                        phspectra1->Fill(energyckov*1e9,(float) 1);
3594                        phot++;
3595                      }
3596                    for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
3597                      if (cmother == nmothers){
3598                        if (closs == 4)
3599                          mothers2[cmother]++;
3600                        mothers[cmother]++;
3601                      }
3602                    }
3603                  }
3604              }
3605            }
3606
3607
3608            if(gAlice->TreeR())
3609              {
3610                Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
3611                gAlice->TreeR()->GetEvent(nent-1);
3612                TClonesArray *Rawclusters = pRICH->RawClustAddress(2);    //  Raw clusters branch
3613                //printf ("Rawclusters:%p",Rawclusters);
3614                Int_t nrawclusters = Rawclusters->GetEntriesFast();
3615
3616                if (nrawclusters) {
3617                  printf("Raw Clusters    : %d\n",nrawclusters);
3618                  for (Int_t hit=0;hit<nrawclusters;hit++) {
3619                    AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
3620                    //Int_t nchamber = rcHit->fChamber;     // chamber number
3621                    //Int_t nhit = cHit->fHitNumber;        // hit number
3622                    Int_t qtot = rcHit->fQ;                 // charge
3623                    Float_t fx  =  rcHit->fX;                 // x-position
3624                    Float_t fy  =  rcHit->fY;                 // y-position
3625                    //Int_t type = rcHit->fCtype;             // cluster type ?
3626                    Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
3627                    pads += mult;
3628                    if (qtot > 0) {
3629                      //printf ("fx: %d, fy: %d\n",fx,fy);
3630                      if (fx>(x-4) && fx<(x+4)  && fy>(y-4) && fy<(y+4)) {
3631                        //printf("There %d \n",mult);
3632                        padmip+=mult;
3633                      } else {
3634                        padnumber->Fill(mult,(float) 1);
3635                        nraw++;
3636                        if (mult<4) Clcharge->Fill(qtot,(float) 1);
3637                      }
3638
3639                    }
3640                  }
3641                }
3642
3643
3644                TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3645                Int_t nrechits1D = RecHits1D->GetEntriesFast();
3646                //printf (" nrechits:%d\n",nrechits);
3647
3648                if(nrechits1D)
3649                  {
3650                    for (Int_t hit=0;hit<nrechits1D;hit++) {
3651                      AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
3652                      Float_t r_omega = recHit1D->fOmega;                  // Cerenkov angle
3653                      Float_t *cer_pho = recHit1D->fCerPerPhoton;        // Cerenkov angle per photon
3654                      Int_t *padsx = recHit1D->fPadsUsedX;           // Pads Used fo reconstruction (x)
3655                      Int_t *padsy = recHit1D->fPadsUsedY;           // Pads Used fo reconstruction (y)
3656                      Int_t goodPhotons = recHit1D->fGoodPhotons;    // Number of pads used for reconstruction
3657
3658                      Omega1D->Fill(r_omega,(float) 1);
3659
3660                      for (Int_t i=0; i<goodPhotons; i++)
3661                        {
3662                          PhotonCer->Fill(cer_pho[i],(float) 1);
3663                          PadsUsed->Fill(padsx[i],padsy[i],1);
3664                          //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
3665                        }
3666
3667                      //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
3668                    }
3669                  }
3670
3671
3672                TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3673                Int_t nrechits3D = RecHits3D->GetEntriesFast();
3674                //printf (" nrechits:%d\n",nrechits);
3675
3676                if(nrechits3D)
3677                  {
3678                    recEffEvent = 0;
3679
3680                    //for (Int_t hit=0;hit<nrechits3D;hit++) {
3681                    AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(track);
3682                    Float_t r_omega    = recHit3D->fOmega;                  // Cerenkov angle
3683                    Float_t r_theta    = recHit3D->fTheta;                  // Theta angle of incidence
3684                    Float_t r_phi      = recHit3D->fPhi;                    // Phi angle if incidence
3685                    Float_t meanradius = recHit3D->fMeanRadius;              // Mean radius for reconstructed point
3686                    Float_t originalOmega = recHit3D->fOriginalOmega;       // Real Cerenkov angle
3687                    Float_t originalTheta = recHit3D->fOriginalTheta;       // Real incidence angle
3688                    Float_t originalPhi = recHit3D->fOriginalPhi;           // Real azimuthal angle
3689
3690
3691                    //correction to track cerenkov angle
3692                    originalOmega = (Float_t) ckovangle->GetMean();
3693
3694                    if(diaglevel == 4)
3695                      {
3696                        printf("\nMean cerenkov angle: %f\n", originalOmega);
3697                        printf("Reconstructed cerenkov angle: %f\n",r_omega);
3698                      }
3699
3700                    Float_t omegaError = TMath::Abs(originalOmega - r_omega);
3701                    Float_t thetaError = TMath::Abs(originalTheta - r_theta);
3702                    Float_t phiError   = TMath::Abs(originalPhi - r_phi);
3703
3704                    //chiSquareOmega += (omegaError/originalOmega)*(omegaError/originalOmega);
3705                    //chiSquareTheta += (thetaError/originalTheta)*(thetaError/originalTheta);
3706                    //chiSquarePhi += (phiError/originalPhi)*(phiError/originalPhi);
3707
3708                    if(TMath::Abs(omegaError) < 0.015)
3709                      recEffEvent += 1;
3710
3711
3712
3713                    //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY);
3714
3715                    Omega3D->Fill(r_omega,(float) 1);
3716                    Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
3717                    Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
3718                    MeanRadius->Fill(meanradius,(float) 1);
3719                    identification->Fill(PTfinal, r_omega,1);
3720                    OriginalOmega->Fill(originalOmega, (float) 1);
3721                    OriginalTheta->Fill(originalTheta, (float) 1);
3722                    OriginalPhi->Fill(TMath::Abs(originalPhi), (float) 1);
3723                    OmegaError->Fill(omegaError, (float) 1);
3724                    ThetaError->Fill(thetaError, (float) 1);
3725                    PhiError->Fill(phiError, (float) 1);
3726
3727                    recEffEvent = recEffEvent;
3728                    recEffTotal += recEffEvent;
3729
3730                    Float_t pioncer = acos(sqrt((.139*.139+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
3731                    Float_t kaoncer = acos(sqrt((.439*.439+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
3732                    Float_t protoncer = acos(sqrt((.938*.938+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
3733
3734                    Float_t piondist = TMath::Abs(r_omega - pioncer);
3735                    Float_t kaondist = TMath::Abs(r_omega - kaoncer);
3736                    Float_t protondist = TMath::Abs(r_omega - protoncer);
3737
3738                    if(diaglevel == 4)
3739                      {
3740                        if(pioncer<r_omega)
3741                          {
3742                            printf("Identified as a PION!\n");
3743                            pionCount += 1;
3744                          }
3745                        if(kaoncer<r_omega && pioncer>r_omega)
3746                          {
3747                            if(kaondist>piondist)
3748                              {
3749                                printf("Identified as a PION!\n");
3750                                pionCount += 1;
3751                              }
3752                            else
3753                              {
3754                                printf("Identified as a KAON!\n");
3755                                kaonCount += 1;
3756                              }
3757                          }                       }
3758                        if(protoncer<r_omega && kaoncer>r_omega)
3759                          {
3760                            if(kaondist>protondist)
3761                              {
3762                                printf("Identified as a PROTON!\n");
3763                                protonCount += 1;
3764                              }
3765                            else
3766                              {
3767                                printf("Identified as a KAON!\n");
3768                                pionCount += 1;
3769                              }
3770                          }
3771                        if(protoncer>r_omega)
3772                          {
3773                            printf("Identified as a PROTON!\n");
3774                            protonCount += 1;
3775                          }
3776
3777                        printf("\nReconstruction efficiency: %5.2f%%\n", recEffEvent*100);
3778                  }
3779              }
3780        }
3781
3782
3783        for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
3784          totalphotonstrack->Fill(mothers[nmothers],(float) 1);
3785          mother->Fill(mothers2[nmothers],(float) 1);
3786          //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
3787        }
3788
3789        clusev->Fill(nraw,(float) 1);
3790        photev->Fill(phot,(float) 1);
3791        feedev->Fill(feed,(float) 1);
3792        padsmip->Fill(padmip,(float) 1);
3793        padscl->Fill(pads,(float) 1);
3794        //printf("Photons:%d\n",phot);
3795        phot = 0;
3796        feed = 0;
3797        pads = 0;
3798        nraw=0;
3799        padmip=0;
3800
3801
3802
3803        gAlice->ResetDigits();
3804        //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
3805        gAlice->TreeD()->GetEvent(0);
3806
3807        if (diaglevel < 4)
3808          {
3809
3810
3811            TClonesArray *Digits  = pRICH->DigitsAddress(2);
3812            Int_t ndigits = Digits->GetEntriesFast();
3813            printf("Digits          : %d\n",ndigits);
3814            padsev->Fill(ndigits,(float) 1);
3815            for (Int_t hit=0;hit<ndigits;hit++) {
3816              AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3817              Int_t qtot = dHit->Signal();                // charge
3818              Int_t ipx  = dHit->PadX();               // pad number on X
3819              Int_t ipy  = dHit->PadY();               // pad number on Y
3820              //printf("%d, %d\n",ipx,ipy);
3821              if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
3822            }
3823          }
3824
3825        if (diaglevel == 5)
3826          {
3827            for (Int_t ich=0;ich<7;ich++)
3828              {
3829                TClonesArray *Digits = pRICH->DigitsAddress(ich);    //  Raw clusters branch
3830                Int_t ndigits = Digits->GetEntriesFast();
3831                //printf("Digits:%d\n",ndigits);
3832                padsev->Fill(ndigits,(float) 1);
3833                if (ndigits) {
3834                  for (Int_t hit=0;hit<ndigits;hit++) {
3835                    AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
3836                    //Int_t nchamber = dHit->GetChamber();     // chamber number
3837                    //Int_t nhit = dHit->fHitNumber;          // hit number
3838                    Int_t qtot = dHit->Signal();                // charge
3839                    Int_t ipx  = dHit->PadX();               // pad number on X
3840                    Int_t ipy  = dHit->PadY();               // pad number on Y
3841                    //Int_t iqpad  = dHit->fQpad;           // charge per pad
3842                    //Int_t rpad  = dHit->fRSec;            // R-position of pad
3843                    //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
3844                    if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
3845                    if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
3846                    if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
3847                    if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
3848                    if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
3849                    if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
3850                    if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
3851                    if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
3852                  }
3853                }
3854              }
3855          }
3856    }
3857
3858    if(diaglevel == 4)
3859      {
3860
3861        Stat_t omegaE;
3862        Stat_t thetaE;
3863        Stat_t phiE;
3864
3865        Stat_t omegaO;
3866        Stat_t thetaO;
3867        Stat_t phiO;
3868
3869        for(Int_t i=0;i<99;i++)
3870          {
3871            omegaE = OriginalOmega->GetBinContent(i);
3872            if(omegaE != 0)
3873              {
3874                omegaO = Omega3D->GetBinContent(i);
3875                chiSquareOmega += (TMath::Power(omegaE,2) - TMath::Power(omegaO,2))/omegaO;
3876              }
3877
3878            thetaE = OriginalTheta->GetBinContent(i);
3879            if(thetaE != 0)
3880              {
3881                thetaO = Theta->GetBinContent(i);
3882                chiSquareTheta += (TMath::Power(thetaE,2) - TMath::Power(thetaO,2))/thetaO;
3883              }
3884
3885            phiE = OriginalPhi->GetBinContent(i);
3886            if(phiE != 0)
3887              {
3888                phiO = Phi->GetBinContent(i);
3889                chiSquarePhi += (TMath::Power(phiE,2) - TMath::Power(phiO,2))/phiO;
3890              }
3891
3892            //printf(" o: %f  t: %f  p: %f\n", OriginalOmega->GetBinContent(i), OriginalTheta->GetBinContent(i),OriginalPhi->GetBinContent(i));
3893
3894          }
3895
3896
3897
3898        printf("\nChi square test values:   Omega - %f\n", chiSquareOmega);
3899        printf("                          Theta - %f\n", chiSquareTheta);
3900        printf("                          Phi   - %f\n", chiSquarePhi);
3901
3902        printf("\nKolmogorov test values:   Omega - %5.4f\n", Omega3D->KolmogorovTest(OriginalOmega));
3903        printf("                          Theta - %5.4f\n", Theta->KolmogorovTest(OriginalTheta));
3904        printf("                          Phi   - %5.4f\n", Phi->KolmogorovTest(OriginalPhi));
3905
3906        recEffTotal = recEffTotal/evNumber2;
3907        printf("\nTotal reconstruction efficiency: %5.2f%%\n", recEffTotal*100);
3908        printf("\n Pions: %d\n Kaons: %d\n Protons:%d\n",pionCount, kaonCount, protonCount);
3909
3910      }
3911
3912
3913    //Create canvases, set the view range, show histograms
3914
3915    TCanvas *c1 = 0;
3916    TCanvas *c2 = 0;
3917    TCanvas *c3 = 0;
3918    TCanvas *c4 = 0;
3919    TCanvas *c5 = 0;
3920    TCanvas *c6 = 0;
3921    TCanvas *c7 = 0;
3922    TCanvas *c8 = 0;
3923    TCanvas *c9 = 0;
3924    TCanvas *c10 = 0;
3925    TCanvas *c11 = 0;
3926    TCanvas *c12 = 0;
3927    TCanvas *c13 = 0;
3928
3929    //TF1* expo = 0;
3930    //TF1* gaus = 0;
3931
3932    TStyle *mystyle=new TStyle("Plain","mystyle");
3933    mystyle->SetPalette(1,0);
3934    //mystyle->SetTitleYSize(0.2);
3935    //mystyle->SetStatW(0.19);
3936    //mystyle->SetStatH(0.1);
3937    //mystyle->SetStatFontSize(0.01);
3938    //mystyle->SetTitleYSize(0.3);
3939    mystyle->SetFuncColor(2);
3940    //mystyle->SetOptStat(0111);
3941    mystyle->SetDrawBorder(0);
3942    mystyle->SetTitleBorderSize(0);
3943    mystyle->SetOptFit(1111);
3944    mystyle->cd();
3945
3946
3947    TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
3948    Int_t nrechits3D = RecHits3D->GetEntriesFast();
3949    TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
3950    Int_t nrechits1D = RecHits1D->GetEntriesFast();
3951
3952   switch(diaglevel)
3953      {
3954      case 1:
3955
3956        c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
3957        hc0->SetXTitle("ix (npads)");
3958        hc0->Draw("colz");
3959
3960 //
3961        c2 = new TCanvas("c2","Hits per type",100,100,600,700);
3962        c2->Divide(2,2);
3963        //c4->SetFillColor(42);
3964
3965        c2->cd(1);
3966        feedback->SetXTitle("x (cm)");
3967        feedback->SetYTitle("y (cm)");
3968        feedback->Draw("colz");
3969
3970        c2->cd(2);
3971        //mip->SetFillColor(5);
3972        mip->SetXTitle("x (cm)");
3973        mip->SetYTitle("y (cm)");
3974        mip->Draw("colz");
3975
3976        c2->cd(3);
3977        //cerenkov->SetFillColor(5);
3978        cerenkov->SetXTitle("x (cm)");
3979        cerenkov->SetYTitle("y (cm)");
3980        cerenkov->Draw("colz");
3981
3982        c2->cd(4);
3983        //h->SetFillColor(5);
3984        h->SetXTitle("x (cm)");
3985        h->SetYTitle("y (cm)");
3986        h->Draw("colz");
3987
3988        c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
3989        c3->Divide(2,1);
3990        //c10->SetFillColor(42);
3991
3992        c3->cd(1);
3993        hitsX->SetFillColor(5);
3994        hitsX->SetXTitle("(cm)");
3995        hitsX->Draw();
3996
3997        c3->cd(2);
3998        hitsY->SetFillColor(5);
3999        hitsY->SetXTitle("(cm)");
4000        hitsY->Draw();
4001
4002
4003        break;
4004 //
4005      case 2:
4006
4007        c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
4008        c4->Divide(2,1);
4009        //c6->SetFillColor(42);
4010
4011        c4->cd(1);
4012        phspectra2->SetFillColor(5);
4013        phspectra2->SetXTitle("energy (eV)");
4014        phspectra2->Draw();
4015        c4->cd(2);
4016        phspectra1->SetFillColor(5);
4017        phspectra1->SetXTitle("energy (eV)");
4018        phspectra1->Draw();
4019
4020        c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
4021        c5->Divide(2,2);
4022        //c9->SetFillColor(42);
4023
4024        c5->cd(1);
4025        pionspectra->SetFillColor(5);
4026        pionspectra->SetXTitle("(GeV)");
4027        pionspectra->Draw();
4028
4029        c5->cd(2);
4030        protonspectra->SetFillColor(5);
4031        protonspectra->SetXTitle("(GeV)");
4032        protonspectra->Draw();
4033
4034        c5->cd(3);
4035        kaonspectra->SetFillColor(5);
4036        kaonspectra->SetXTitle("(GeV)");
4037        kaonspectra->Draw();
4038
4039        c5->cd(4);
4040        chargedspectra->SetFillColor(5);
4041        chargedspectra->SetXTitle("(GeV)");
4042        chargedspectra->Draw();
4043
4044        break;
4045
4046      case 3:
4047
4048
4049        if(gAlice->TreeR())
4050          {
4051            c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
4052            c6->Divide(2,2);
4053            //c3->SetFillColor(42);
4054
4055            c6->cd(1);
4056            //TPad* c6_1;
4057            //c6_1->SetLogy();
4058            Clcharge->SetFillColor(5);
4059            Clcharge->SetXTitle("ADC counts");
4060            if (evNumber2>10)
4061              {
4062                Clcharge->Fit("expo");
4063                //expo->SetLineColor(2);
4064                //expo->SetLineWidth(3);
4065              }
4066            Clcharge->Draw();
4067
4068            c6->cd(2);
4069            padnumber->SetFillColor(5);
4070            padnumber->SetXTitle("(counts)");
4071            padnumber->Draw();
4072
4073            c6->cd(3);
4074            clusev->SetFillColor(5);
4075            clusev->SetXTitle("(counts)");
4076            if (evNumber2>10)
4077              {
4078                clusev->Fit("gaus");
4079                //gaus->SetLineColor(2);
4080                //gaus->SetLineWidth(3);
4081              }
4082            clusev->Draw();
4083
4084            c6->cd(4);
4085            padsmip->SetFillColor(5);
4086            padsmip->SetXTitle("(counts)");
4087            padsmip->Draw();
4088          }
4089
4090        if(evNumber2<1)
4091          {
4092            c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
4093            mother->SetFillColor(5);
4094            mother->SetXTitle("counts");
4095            mother->Draw();
4096          }
4097
4098        c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
4099        c7->Divide(2,2);
4100        //c7->SetFillColor(42);
4101
4102        c7->cd(1);
4103        totalphotonsevent->SetFillColor(5);
4104        totalphotonsevent->SetXTitle("Photons (counts)");
4105        if (evNumber2>10)
4106            {
4107              totalphotonsevent->Fit("gaus");
4108              //gaus->SetLineColor(2);
4109              //gaus->SetLineWidth(3);
4110            }
4111        totalphotonsevent->Draw();
4112
4113        c7->cd(2);
4114        photev->SetFillColor(5);
4115        photev->SetXTitle("(counts)");
4116        if (evNumber2>10)
4117          {
4118            photev->Fit("gaus");
4119            //gaus->SetLineColor(2);
4120            //gaus->SetLineWidth(3);
4121          }
4122        photev->Draw();
4123
4124        c7->cd(3);
4125        feedev->SetFillColor(5);
4126        feedev->SetXTitle("(counts)");
4127        if (evNumber2>10)
4128          {
4129            feedev->Fit("gaus");
4130            //gaus->SetLineColor(2);
4131            //gaus->SetLineWidth(3);
4132          }
4133        feedev->Draw();
4134
4135        c7->cd(4);
4136        padsev->SetFillColor(5);
4137        padsev->SetXTitle("(counts)");
4138        if (evNumber2>10)
4139          {
4140            padsev->Fit("gaus");
4141            //gaus->SetLineColor(2);
4142            //gaus->SetLineWidth(3);
4143          }
4144        padsev->Draw();
4145
4146        break;
4147
4148      case 4:
4149
4150
4151        if(nrechits3D)
4152          {
4153            c8 = new TCanvas("c8","3D reconstruction of Phi angle",50,50,300,1050);
4154            c8->Divide(1,3);
4155            //c2->SetFillColor(42);
4156
4157
4158            // data per hit
4159            c8->cd(1);
4160            hitsPhi->SetFillColor(5);
4161            if (evNumber2>10)
4162              hitsPhi->Fit("gaus");
4163            hitsPhi->Draw();
4164
4165             //data per track
4166            c8->cd(2);
4167            OriginalPhi->SetFillColor(5);
4168            if (evNumber2>10)
4169              OriginalPhi->Fit("gaus");
4170            OriginalPhi->Draw();
4171
4172            //recontructed data
4173            c8->cd(3);
4174            Phi->SetFillColor(5);
4175            if (evNumber2>10)
4176              Phi->Fit("gaus");
4177            Phi->Draw();
4178
4179            c9 = new TCanvas("c9","3D reconstruction of theta angle",75,75,300,1050);
4180            c9->Divide(1,3);
4181
4182            // data per hit
4183            c9->cd(1);
4184            hitsTheta->SetFillColor(5);
4185            if (evNumber2>10)
4186              hitsTheta->Fit("gaus");
4187            hitsTheta->Draw();
4188
4189            //data per track
4190            c9->cd(2);
4191            OriginalTheta->SetFillColor(5);
4192            if (evNumber2>10)
4193              OriginalTheta->Fit("gaus");
4194            OriginalTheta->Draw();
4195
4196            //recontructed data
4197            c9->cd(3);
4198            Theta->SetFillColor(5);
4199            if (evNumber2>10)
4200              Theta->Fit("gaus");
4201            Theta->Draw();
4202
4203            c10 = new TCanvas("c10","3D reconstruction of cherenkov angle",100,100,300,1050);
4204            c10->Divide(1,3);
4205
4206            // data per hit
4207            c10->cd(1);
4208            ckovangle->SetFillColor(5);
4209            ckovangle->SetXTitle("angle (radians)");
4210            if (evNumber2>10)
4211              ckovangle->Fit("gaus");
4212            ckovangle->Draw();
4213
4214            //data per track
4215            c10->cd(2);
4216            OriginalOmega->SetFillColor(5);
4217            OriginalOmega->SetXTitle("angle (radians)");
4218            if (evNumber2>10)
4219              OriginalOmega->Fit("gaus");
4220            OriginalOmega->Draw();
4221
4222            //recontructed data
4223            c10->cd(3);
4224            Omega3D->SetFillColor(5);
4225            Omega3D->SetXTitle("angle (radians)");
4226            if (evNumber2>10)
4227              Omega3D->Fit("gaus");
4228            Omega3D->Draw();
4229
4230
4231            c11 = new TCanvas("c11","3D reconstruction of mean radius",125,125,300,700);
4232            c11->Divide(1,2);
4233
4234            // data per hit
4235            c11->cd(1);
4236            radius->SetFillColor(5);
4237            radius->SetXTitle("radius (cm)");
4238            radius->Draw();
4239
4240            //recontructed data
4241            c11->cd(2);
4242            MeanRadius->SetFillColor(5);
4243            MeanRadius->SetXTitle("radius (cm)");
4244            MeanRadius->Draw();
4245
4246
4247            c12 = new TCanvas("c12","Cerenkov angle vs. Momentum",150,150,550,350);
4248
4249            c12->cd(1);
4250            identification->SetFillColor(5);
4251            identification->SetXTitle("Momentum (GeV/c)");
4252            identification->SetYTitle("Cherenkov angle (radians)");
4253
4254            //Float_t pionmass=.139;
4255            //Float_t kaonmass=.493;
4256            //Float_t protonmass=.938;
4257            //Float_t n=1.295;
4258
4259            TF1 *pionplot = new TF1("pion","acos(sqrt((.139*.139+x*x)/(x*x*1.285*1.285)))",1,5);
4260            TF1 *kaonplot = new TF1("kaon","acos(sqrt((.439*.439+x*x)/(x*x*1.285*1.285)))",1,5);
4261            TF1 *protonplot = new TF1("proton","acos(sqrt((.938*.938+x*x)/(x*x*1.285*1.285)))",1,5);
4262
4263            identification->Draw();
4264
4265            pionplot->SetLineColor(5);
4266            pionplot->Draw("same");
4267
4268            kaonplot->SetLineColor(4);
4269            kaonplot->Draw("same");
4270
4271            protonplot->SetLineColor(3);
4272            protonplot->Draw("same");
4273            //identification->Draw("same");
4274
4275
4276
4277            c13 = new TCanvas("c13","Reconstruction Errors",200,200,900,350);
4278            c13->Divide(3,1);
4279
4280            c13->cd(1);
4281            PhiError->SetFillColor(5);
4282            if (evNumber2>10)
4283              PhiError->Fit("gaus");
4284            PhiError->Draw();
4285            c13->cd(2);
4286            ThetaError->SetFillColor(5);
4287            if (evNumber2>10)
4288              ThetaError->Fit("gaus");
4289            ThetaError->Draw();
4290            c13->cd(3);
4291            OmegaError->SetFillColor(5);
4292            OmegaError->SetXTitle("angle (radians)");
4293            if (evNumber2>10)
4294              OmegaError->Fit("gaus");
4295            OmegaError->Draw();
4296
4297          }
4298
4299        if(nrechits1D)
4300          {
4301            c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
4302            c9->Divide(3,2);
4303            //c5->SetFillColor(42);
4304
4305            c9->cd(1);
4306            ckovangle->SetFillColor(5);
4307            ckovangle->SetXTitle("angle (radians)");
4308            ckovangle->Draw();
4309
4310            c9->cd(2);
4311            radius->SetFillColor(5);
4312            radius->SetXTitle("radius (cm)");
4313            radius->Draw();
4314
4315            c9->cd(3);
4316            hc0->SetXTitle("pads");
4317            hc0->Draw("box");
4318
4319            c9->cd(5);
4320            Omega1D->SetFillColor(5);
4321            Omega1D->SetXTitle("angle (radians)");
4322            Omega1D->Draw();
4323
4324            c9->cd(4);
4325            PhotonCer->SetFillColor(5);
4326            PhotonCer->SetXTitle("angle (radians)");
4327            PhotonCer->Draw();
4328
4329            c9->cd(6);
4330            PadsUsed->SetXTitle("pads");
4331            PadsUsed->Draw("box");
4332          }
4333
4334        break;
4335
4336      case 5:
4337
4338        printf("Drawing histograms.../n");
4339
4340        //if (ndigits)
4341          //{
4342        c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
4343        c1->Divide(4,2);
4344        //c1->SetFillColor(42);
4345
4346        c10->cd(1);
4347        hc1->SetXTitle("ix (npads)");
4348        hc1->Draw("box");
4349        c10->cd(2);
4350        hc2->SetXTitle("ix (npads)");
4351        hc2->Draw("box");
4352        c10->cd(3);
4353        hc3->SetXTitle("ix (npads)");
4354        hc3->Draw("box");
4355        c10->cd(4);
4356        hc4->SetXTitle("ix (npads)");
4357        hc4->Draw("box");
4358        c10->cd(5);
4359        hc5->SetXTitle("ix (npads)");
4360        hc5->Draw("box");
4361        c10->cd(6);
4362        hc6->SetXTitle("ix (npads)");
4363        hc6->Draw("box");
4364        c10->cd(7);
4365        hc7->SetXTitle("ix (npads)");
4366        hc7->Draw("box");
4367        c10->cd(8);
4368        hc0->SetXTitle("ix (npads)");
4369        hc0->Draw("box");
4370          //}
4371 //
4372        c11 = new TCanvas("c11","Hits per type",100,100,600,700);
4373        c11->Divide(2,2);
4374        //c4->SetFillColor(42);
4375
4376        c11->cd(1);
4377        feedback->SetXTitle("x (cm)");
4378        feedback->SetYTitle("y (cm)");
4379        feedback->Draw();
4380
4381        c11->cd(2);
4382        //mip->SetFillColor(5);
4383        mip->SetXTitle("x (cm)");
4384        mip->SetYTitle("y (cm)");
4385        mip->Draw();
4386
4387        c11->cd(3);
4388        //cerenkov->SetFillColor(5);
4389        cerenkov->SetXTitle("x (cm)");
4390        cerenkov->SetYTitle("y (cm)");
4391        cerenkov->Draw();
4392
4393        c11->cd(4);
4394        //h->SetFillColor(5);
4395        h->SetXTitle("x (cm)");
4396        h->SetYTitle("y (cm)");
4397        h->Draw();
4398
4399        c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
4400        c12->Divide(2,1);
4401        //c10->SetFillColor(42);
4402
4403        c12->cd(1);
4404        hitsX->SetFillColor(5);
4405        hitsX->SetXTitle("(cm)");
4406        hitsX->Draw();
4407
4408        c12->cd(2);
4409        hitsY->SetFillColor(5);
4410        hitsY->SetXTitle("(cm)");
4411        hitsY->Draw();
4412
4413        break;
4414
4415      }
4416
4417
4418    // calculate the number of pads which give a signal
4419
4420
4421    //Int_t Np=0;
4422    /*for (Int_t i=0;i< NpadX;i++) {
4423        for (Int_t j=0;j< NpadY;j++) {
4424            if (Pad[i][j]>=6){
4425                Np+=1;
4426            }
4427        }
4428    }*/
4429    //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
4430    printf("\nEnd of analysis\n");
4431    printf("**********************************\n");
4432 }
4433
4434 ////////////////////////////////////////////////////////////////////////
4435 void AliRICH::MakeBranchInTreeD(TTree *treeD, const char *file)
4436 {
4437     //
4438     // Create TreeD branches for the RICH.
4439     //
4440
4441   const Int_t kBufferSize = 4000;
4442   char branchname[30];
4443
4444   //
4445   // one branch for digits per chamber
4446   //
4447   for (Int_t i=0; i<kNCH ;i++) {
4448     sprintf(branchname,"%sDigits%d",GetName(),i+1);
4449     if (fDchambers && treeD) {
4450       MakeBranchInTree(treeD,
4451                        branchname, &((*fDchambers)[i]), kBufferSize, file);
4452 //      printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
4453     }
4454   }
4455 }
4456 ////////////////////////////////////////////////////////////////////////