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