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