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