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