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