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