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