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