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