all old staff moved to v3
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Sep 2003 13:48:18 +0000 (13:48 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 17 Sep 2003 13:48:18 +0000 (13:48 +0000)
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHv0.cxx
RICH/AliRICHv3.cxx
RICH/AliRICHv3.h

index f0cff20..d9a31a0 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id$ */
 
-////////////////////////////////////////////////
-//  Manager and hits classes for set:RICH     //
-////////////////////////////////////////////////
-
-#include <Riostream.h>
 #include <strings.h>
 
 #include <TArrayF.h>
 #include <TBRIK.h>
-#include <TCanvas.h>
-#include <TF1.h>
 #include <TFile.h>
 #include <TGeometry.h>
-#include <TH1.h>
-#include <TH2.h>
 #include <TNode.h> 
 #include <TObjArray.h>
 #include <TObject.h>
 #include <TParticle.h>
 #include <TPDGCode.h>
 #include <TRandom.h> 
-#include <TStyle.h>
 #include <TTUBE.h>
 #include <TTree.h>
 #include <TVector.h>
-#include "AliConst.h"
 #include "AliMagF.h"
 #include "AliPoints.h"
 #include "AliRICH.h"
+#include "AliRICHParam.h"
 #include "AliRICHClusterFinder.h"
-#include "AliRICHDigit.h"
 #include "AliRICHDigitizer.h"
 #include "AliRICHHitMapA1.h"
 #include "AliRICHMerger.h"
 #include "AliRICHTransientDigit.h"
 #include "AliRun.h"
 #include "AliRunDigitizer.h"
-#include "AliSegmentation.h"
-#include "AliRICHParam.h"
-
-static Int_t sMaxIterPad=0;    // Static variables for the pad-hit iterator routines
-static Int_t sCurIterPad=0;
+#include "AliRICHSegmentationV1.h"
+#include "AliRICHResponseV0.h"
  
 ClassImp(AliRICHhit)
 ClassImp(AliRICHdigit)
@@ -95,42 +80,37 @@ AliRICH::AliRICH()
       fNrechits3D[i] = 0;
   }
   fpParam=0;
-//kir  fFileName = 0;
-//kir  fMerger = 0;
 }//AliRICH::AliRICH()
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRICH::AliRICH(const char *name, const char *title)
         :AliDetector(name,title)
 {//Named ctor
   if(GetDebug())Info("named ctor","Start.");
   fpParam     =new AliRICHParam;
-  fHits       =new TClonesArray("AliRICHhit",1000  );
-  fCerenkovs  =new TClonesArray("AliRICHCerenkov",1000);
-  fSDigits    =new TClonesArray("AliRICHdigit",100000);
-  gAlice->AddHitList(fHits);
-  gAlice->AddHitList(fCerenkovs);
-  fNsdigits   =0;
-  fNcerenkovs =0;
+  fNcerenkovs =fNsdigits   =0;//fNhits and fNdigits reset in AliDetector ctor
+  fHits       =new TClonesArray("AliRICHhit",1000  ); gAlice->AddHitList(fHits);//hits
+  fCerenkovs  =new TClonesArray("AliRICHCerenkov",1000);gAlice->AddHitList(fCerenkovs);//cerenkovs ??? to be removed    
+  fSDigits    =new TClonesArray("AliRICHdigit",100000);//sdigits
+  fDigits     =new TClonesArray("AliRICHdigit",100000);//digits
+  
   fIshunt     =0;
-  fDchambers  =new TObjArray(kNCH);
-  fRawClusters=new TObjArray(kNCH);
-  fRecHits1D  =new TObjArray(kNCH);
-  fRecHits3D  =new TObjArray(kNCH);
+  fDchambers  =new TObjArray(kNCH);//digits             ??? to be removed
+  fRawClusters=new TObjArray(kNCH);//clusters
+  fRecHits1D  =new TObjArray(kNCH);//recos Bari
+  fRecHits3D  =new TObjArray(kNCH);//recos Lisbon
   for(int i=0;i<kNCH;i++) {
-    fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i); 
+    fNdch[i]=fNrawch[i]=0;
+    fDchambers  ->AddAt(new TClonesArray("AliRICHDigit",10000), i); //??? to be removed
     fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i); 
-    fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
-    fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
-    fNdch[i]=0;
-    fNrawch[i]=0;
+    fRecHits1D  ->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
+    fRecHits3D  ->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
   }
   SetMarkerColor(kRed);
   fCkovNumber=fFreonProd=0;
-//kir  fFileName = 0;
-//kir  fMerger = 0;
+  CreateChambers();
   if(GetDebug())Info("named ctor","Stop.");
 }//AliRICH::AliRICH(const char *name, const char *title)
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRICH::~AliRICH()
 {//dtor
   if(GetDebug()) Info("dtor","Start.");
@@ -163,92 +143,65 @@ AliRICH::~AliRICH()
     }                     
   if(GetDebug()) Info("dtor","Stop.");    
 }//AliRICH::~AliRICH()
-//______________________________________________________________________________
-Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
-{//calls the charge disintegration method of the current chamber and adds all generated sdigits to the list of digits
-   
-   Int_t iChamber,iPadX,iPadY,iAdc,iTrack;
-   Float_t list[4][500];
-   Int_t iNdigits;
-        
-
-  ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, iNdigits, list, res);
-    Int_t ic=0;
-    
-  for(Int_t i=0; i<iNdigits; i++) {
-    if(Int_t(list[0][i]) > 0) {
-           ic++;
-           iAdc = Int_t(list[0][i]);
-           iPadX = Int_t(list[1][i]);
-           iPadY = Int_t(list[2][i]);
-           iChamber = Int_t(list[3][i]);
-
-           
-           AddSDigit(iChamber,iPadX,iPadY,iAdc,iTrack);
-       }
-    }
-    
-   if(gAlice->TreeS()){
-       gAlice->TreeS()->Fill();
-       gAlice->TreeS()->Write(0,TObject::kOverwrite);
-   }
-   return iNdigits;
-}//Int_t AliRICH::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::Hits2SDigits()
 {//Create a list of sdigits corresponding to list of hits. Every hit generates sdigit.
   if(GetDebug()) Info("Hit2SDigits","Start.");
   
-  for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//loop on events
+  for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
     fLoader->GetRunLoader()->GetEvent(iEventN);
   
     if(!fLoader->TreeH()) fLoader->LoadHits();
     if(!fLoader->TreeS()) fLoader->MakeTree("S");
     MakeBranch("S");
-  
-    for(int iPrimN=0;iPrimN<TreeH()->GetEntries();iPrimN++){//loop on primary tracks
+    
+    AliRICHSegmentationV1 *pSeg=new AliRICHSegmentationV1;
+    AliRICHResponseV0     *pRes=new AliRICHResponseV0;
+    
+    Float_t dx=Param()->SigmaIntegration()*Param()->ChargeSpreadX();
+    Float_t dy=Param()->SigmaIntegration()*Param()->ChargeSpreadY();
+    Float_t charge;
+    
+    for(int iPrimN=0;iPrimN<TreeH()->GetEntries();iPrimN++){//prims loop
       fLoader->TreeH()->GetEntry(iPrimN); 
-      for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//loop on hits for given primary track  
-        AddSDigit(4,13,24,55,4);//chamber-xpad-ypad-qdc-track1-2-3
-      }//loop on hits for given primary track
-    }//loop on primary tracks
-  
+      for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop
+        AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit
+        
+        if(pHit->Y()>0)
+          charge=Param()->TotalCharge(pHit->Particle(),pHit->Eloss(),pHit->Y()-Param()->SectorSizeY());
+        else
+          charge=Param()->TotalCharge(pHit->Particle(),pHit->Eloss(),pHit->Y()+Param()->SectorSizeY());
+                  
+        for(pSeg->FirstPad(pHit->X(),pHit->Y(),0,dx,dy);pSeg->MorePads();pSeg->NextPad()){//pads loop
+            
+          AddSDigit(pHit->Chamber(),pSeg->Ix(),pSeg->Iy(),
+                                    Int_t(charge*TMath::Abs(pRes->IntXY(pSeg))),
+                iPrimN);//chamber-xpad-ypad-qdc-track1-2-3
+        }//pads loop
+      }//hits loop
+    }//prims loop
+    
+    delete pSeg;
+    
     fLoader->TreeS()->Fill();
     fLoader->WriteSDigits("OVERWRITE");
-  }//loop on events
+  }//events loop
   
+  fLoader->UnloadHits();
+  fLoader->UnloadSDigits();  
   if(GetDebug()) Info("Hit2SDigits","Stop.");
-}
-//______________________________________________________________________________
+}//void AliRICH::Hits2SDigits()
+//__________________________________________________________________________________________________
 void AliRICH::SDigits2Digits()
 {//Generate digits from sdigits.
   if(GetDebug()) Info("SDigits2Digits","Start.");
-   //AliRICHChamber*       iChamber;
-  
-  
-   //for(Int_t i=0;i<7;i++) {
-   //iChamber = &(Chamber(i));
-   //iChamber->GenerateTresholds();
-   //}
-  
-   //int nparticles = gAlice->GetNtrack();
-   //cout << "Particles (RICH):" <<nparticles<<endl;
-   //if (nparticles <= 0) return;
-   //if (!fMerger) {
-   //fMerger = new AliRICHMerger();
-   //}
-
-
-   //fMerger->Init();
-   //fMerger->Digitise(nev,flag);
-
-   AliRunDigitizer * manager = new AliRunDigitizer(1,1);
-   manager->SetInputStream(0,"galice.root");
-   //AliRICHDigitizer *dRICH  = new AliRICHDigitizer(manager);
-   manager->Exec("deb");
+    
+  AliRunDigitizer *pManager = new AliRunDigitizer(1,1);
+  pManager->SetInputStream(0,"galice.root");
+  pManager->Exec("deb");
   if(GetDebug()) Info("SDigits2Digits","Stop.");
 }//void AliRICH::SDigits2Digits()
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::Digits2Reco()
 {
 // Generate clusters
@@ -260,6 +213,7 @@ void AliRICH::Digits2Reco()
   if (nparticles > 0) FindClusters(0);
 
 }//void AliRICH::Digits2Reco()  
+//__________________________________________________________________________________________________
 
 
 void AliRICH::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
@@ -299,13 +253,13 @@ void AliRICH::BuildGeometry()
   
   new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
 
-  Float_t wid=fpParam->PadPlaneWidth();
-  Float_t len=fpParam->PadPlaneLength();
+  Float_t wid=fpParam->SectorSizeX();
+  Float_t len=fpParam->SectorSizeY();
   new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2);
   
-  for(int i=0;i<kNCH;i++){
+  for(int i=1;i<=kNCH;i++){
     top->cd();
-    node = new TNode(Form("RICH%i",i+1),Form("RICH%i",i+1),"S_RICH",C(i)->X(),C(i)->Y(),C(i)->Z(),C(i)->RotMatrixName());
+    node = new TNode(Form("RICH%i",i),Form("RICH%i",i),"S_RICH",C(i)->X(),C(i)->Y(),C(i)->Z(),C(i)->RotMatrixName());
     node->SetLineColor(kRed);
     node->cd();
     subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+fpParam->DeadZone(),5,len/2+fpParam->DeadZone()/2,"");
@@ -644,7 +598,7 @@ void AliRICH::CreateMaterials()
     gMC->SetCerenkov(idtmed[1009], 26, ppckov, abscoGrid, efficGrid, rIndexGrid);
     gMC->SetCerenkov(idtmed[1010], 26, ppckov, abscoOpaqueQuarz, efficAll, rIndexOpaqueQuarz);
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
 {
 
@@ -653,13 +607,10 @@ Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
     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,
                      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,
                      7.8,7.9,8.0,8.1,8.2,8.3,8.4,8.5};
-     
-
     Float_t csin[36] = {2.14,2.21,2.33,2.48,2.76,2.97,2.99,2.59,2.81,3.05,
                        2.86,2.53,2.55,2.66,2.79,2.96,3.18,3.05,2.84,2.81,2.38,2.11,
                        2.01,2.13,2.39,2.73,3.08,3.15,2.95,2.73,2.56,2.41,2.12,1.95,
                        1.72,1.53};
-      
     Float_t csik[36] = {0.,0.,0.,0.,0.,0.196,0.408,0.208,0.118,0.49,0.784,0.543,
                        0.424,0.404,0.371,0.514,0.922,1.102,1.139,1.376,1.461,1.253,0.878,
                        0.69,0.612,0.649,0.824,1.347,1.571,1.678,1.763,1.857,1.824,1.824,
@@ -703,9 +654,8 @@ Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
       
     fresn = fresn*rO;
     return(fresn);
-}
-
-//__________________________________________
+}//Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)
+//__________________________________________________________________________________________________
 Float_t AliRICH::AbsoCH4(Float_t x)
 {
 
@@ -848,11 +798,7 @@ Float_t AliRICH::AbsoCH4(Float_t x)
     Float_t alength=1./(1./abslm+1./absli+1./abslo+1./abslw);
     return (alength);
 }
-
-
-
-//___________________________________________
-//____________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::ResetDigits()
 {//Reset number of digits and the digits array for this detector
   for ( int i=0;i<kNCH;i++ ) {
@@ -860,7 +806,7 @@ void AliRICH::ResetDigits()
     if (fNdch)  fNdch[i]=0;
   }
 }
-//____________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::ResetRawClusters()
 {//Reset number of raw clusters and the raw clust array for this detector
   for ( int i=0;i<kNCH;i++ ) {
@@ -868,7 +814,7 @@ void AliRICH::ResetRawClusters()
     if (fNrawch)  fNrawch[i]=0;
   }
 }
-//____________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::ResetRecHits1D()
 {//Reset number of raw clusters and the raw clust array for this detector
   for ( int i=0;i<kNCH;i++ ) {
@@ -877,7 +823,7 @@ void AliRICH::ResetRecHits1D()
   }
 }
 
-//____________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::ResetRecHits3D()
 {// Reset number of raw clusters and the raw clust array for this detector
   for ( int i=0;i<kNCH;i++ ) {
@@ -885,14 +831,13 @@ void AliRICH::ResetRecHits3D()
     if (fNrechits3D)  fNrechits3D[i]=0;
   }
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::FindClusters(Int_t nev /*kir,Int_t lastEntry*/)
 {// Loop on chambers and on cathode planes
     for (Int_t icat=1;icat<2;icat++) {
        gAlice->ResetDigits();
        gAlice->TreeD()->GetEvent(0);
        for (Int_t ich=0;ich<kNCH;ich++) {
-      //PH       AliRICHChamber* iChamber=(AliRICHChamber*) (*fChambers)[ich];
          AliRICHChamber* iChamber=(AliRICHChamber*)fChambers->At(ich);
          TClonesArray *pRICHdigits  = this->DigitsAddress(ich);
          if (pRICHdigits == 0)       
@@ -932,1732 +877,7 @@ void AliRICH::FindClusters(Int_t nev /*kir,Int_t lastEntry*/)
     gAlice->TreeR()->Write(hname,kOverwrite,0);
     gAlice->TreeR()->Reset();    
 }//void AliRICH::FindClusters(Int_t nev)
-//______________________________________________________________________________
-AliRICHSDigit* AliRICH::FirstPad(AliRICHhit*  hit,TClonesArray *clusters ) 
-{// Initialise the pad iterator Return the address of the first sdigit for hit
-    TClonesArray *theClusters = clusters;
-    Int_t nclust = theClusters->GetEntriesFast();
-    if (nclust && hit->PHlast() > 0) {
-       sMaxIterPad=Int_t(hit->PHlast());
-       sCurIterPad=Int_t(hit->PHfirst());
-       return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
-    } else {
-       return 0;
-    }
-    
-}
-//______________________________________________________________________________
-AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters) 
-{// Iterates over pads
-  
-    sCurIterPad++;
-    if (sCurIterPad <= sMaxIterPad) {
-       return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
-    } else {
-       return 0;
-    }
-}
-
-
-void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
-{
-  
-  Int_t NpadX = 162;                 // number of pads on X
-  Int_t NpadY = 162;                 // number of pads on Y
-  
-  Int_t Pad[162][162];
-  for (Int_t i=0;i<NpadX;i++) {
-    for (Int_t j=0;j<NpadY;j++) {
-      Pad[i][j]=0;
-    }
-  }
-  
-  //  Create some histograms
-
-  TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
-  TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
-  TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
-  TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
-  TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
-  TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
-  TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
-  TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
-  TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
-  TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
-  TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
-  TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
-  TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
-  TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
-  TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
-  TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
-  TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
-  TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
-  TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
-  TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
-  TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
-  TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
-  TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
-  TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
-  TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
-  //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
-  TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
-  TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
-  TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
-  TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
-  TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
-  TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
-   
-   
-   
-
-//   Start loop over events 
-
-  Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
-  Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
-  Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
-  TRandom* random=0;
-
-   for (int nev=0; nev<= evNumber2; nev++) {
-       Int_t nparticles = gAlice->GetEvent(nev);
-       
-
-       if (nev < evNumber1) continue;
-       if (nparticles <= 0) return;
-       
-// Get pointers to RICH detector and Hits containers
-       
-       AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
-     
-       TTree *treeH = TreeH();
-       Int_t ntracks =(Int_t) treeH->GetEntries();
-            
-// Start loop on tracks in the hits containers
-       
-       for (Int_t track=0; track<ntracks;track++) {
-          printf ("Processing Track: %d\n",track);
-          gAlice->ResetHits();
-          treeH->GetEvent(track);
-                          
-          for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1); 
-              mHit;
-              mHit=(AliRICHhit*)pRICH->NextHit()) 
-            {
-              //Int_t nch  = mHit->fChamber;              // chamber number
-              //Float_t x  = mHit->X();                    // x-pos of hit
-              //Float_t y  = mHit->Z();                    // y-pos
-              //Float_t z  = mHit->Y();
-              //Float_t phi = mHit->Phi();                 //Phi angle of incidence
-              Float_t theta = mHit->Theta();             //Theta angle of incidence
-              Float_t px = mHit->MomX();
-              Float_t py = mHit->MomY();
-              Int_t index = mHit->Track();
-              Int_t particle = (Int_t)(mHit->Particle());    
-              Float_t R;
-              Float_t PTfinal;
-              Float_t PTvertex;
-
-             TParticle *current = gAlice->Particle(index);
-             
-             //Float_t energy=current->Energy(); 
-
-             R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
-             PTfinal=TMath::Sqrt(px*px + py*py);
-             PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
-             
-             
-
-             if (TMath::Abs(particle) < 10000000)
-               {
-                 hitsTheta->Fill(theta,(float) 1);
-                 if (R<5)
-                   {
-                     if (PTvertex>.5 && PTvertex<=1)
-                       {
-                         hitsTheta500MeV->Fill(theta,(float) 1);
-                       }
-                     if (PTvertex>1 && PTvertex<=2)
-                       {
-                         hitsTheta1GeV->Fill(theta,(float) 1);
-                       }
-                     if (PTvertex>2 && PTvertex<=3)
-                       {
-                         hitsTheta2GeV->Fill(theta,(float) 1);
-                       }
-                     if (PTvertex>3)
-                       {
-                         hitsTheta3GeV->Fill(theta,(float) 1);
-                       }
-                   }
-                 
-               }
-
-             //if (nch == 3)
-               //{
-             
-             if (TMath::Abs(particle) < 50000051)
-               {
-                 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
-                 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
-                   {
-                     //gMC->Rndm(&random, 1);
-                     if (random->Rndm() < .1)
-                       production->Fill(current->Vz(),R,(float) 1);
-                     if (TMath::Abs(particle) == 50000050)
-                       //if (TMath::Abs(particle) > 50000000)
-                       {
-                         photons +=1;
-                         if (R<5)
-                           {
-                             primaryphotons +=1;
-                             if (current->Energy()>0.001)
-                               highprimaryphotons +=1;
-                           }
-                       }       
-                     if (TMath::Abs(particle) == 2112)
-                       {
-                         neutron +=1;
-                         if (current->Energy()>0.0001)
-                           highneutrons +=1;
-                       }
-                   }
-                 if (TMath::Abs(particle) < 50000000)
-                   {
-                     production->Fill(current->Vz(),R,(float) 1);
-                   }
-                 //mip->Fill(x,y,(float) 1);
-               }
-             
-             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
-               {
-                 if (R<5)
-                   {
-                     pionptspectravertex->Fill(PTvertex,(float) 1);
-                     pionptspectrafinal->Fill(PTfinal,(float) 1);
-                   }
-               }
-             
-             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
-                 || TMath::Abs(particle)==311)
-               {
-                 if (R<5)
-                   {
-                     kaonptspectravertex->Fill(PTvertex,(float) 1);
-                     kaonptspectrafinal->Fill(PTfinal,(float) 1);
-                   }
-               }
-             
-             
-             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
-               {
-                 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   {
-                     pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                   }
-                 pion +=1;
-                 if (TMath::Abs(particle)==211)
-                   {
-                     chargedpions +=1;
-                     if (R<5)
-                       {
-                         primarypions +=1;
-                         if (current->Energy()>1)
-                           highprimarypions +=1;
-                       }
-                   }   
-               }
-             if (TMath::Abs(particle)==2212)
-               {
-                 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 //ptspectra->Fill(Pt,(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 proton +=1;
-               }
-             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
-                 || TMath::Abs(particle)==311)
-               {
-                 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 //ptspectra->Fill(Pt,(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 kaon +=1;
-                 if (TMath::Abs(particle)==321)
-                   {
-                     chargedkaons +=1;
-                     if (R<5)
-                       {
-                         primarykaons +=1;
-                         if (current->Energy()>1)
-                           highprimarykaons +=1;
-                       }
-                   }
-               }
-             if (TMath::Abs(particle)==11)
-               {
-                 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 //ptspectra->Fill(Pt,(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (particle == 11)
-                   electron +=1;
-                 if (particle == -11)
-                   positron +=1;
-               }
-             if (TMath::Abs(particle)==13)
-               {
-                 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 //ptspectra->Fill(Pt,(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 muon +=1;
-               }
-             if (TMath::Abs(particle)==2112)
-               {
-                 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 //ptspectra->Fill(Pt,(float) 1);
-                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                   neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                 if (R>250 && R<450)
-                   {
-                     neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                   }
-                 neutron +=1;
-               }
-             if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
-               {
-                 if (current->Energy()-current->GetCalcMass()>1)
-                   {
-                     chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                     if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
-                       chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                     if (R>250 && R<450)
-                       chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
-                   }
-               }
-             // Fill the histograms
-             //Nh1+=nhits;
-             //h->Fill(x,y,(float) 1);
-             //}
-             //}
-          }          
-          
-       }
-       
-   }
-   //   }
-
-   TStyle *mystyle=new TStyle("Plain","mystyle");
-   mystyle->SetPalette(1,0);
-   mystyle->cd();
-   
-   //Create canvases, set the view range, show histograms
-
-    TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
-    c2->Divide(2,2);
-    //c2->SetFillColor(42);
-    
-    c2->cd(1);
-    hitsTheta500MeV->SetFillColor(5);
-    hitsTheta500MeV->Draw();
-    c2->cd(2);
-    hitsTheta1GeV->SetFillColor(5);
-    hitsTheta1GeV->Draw();
-    c2->cd(3);
-    hitsTheta2GeV->SetFillColor(5);
-    hitsTheta2GeV->Draw();
-    c2->cd(4);
-    hitsTheta3GeV->SetFillColor(5);
-    hitsTheta3GeV->Draw();
-    
-            
-   
-    TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
-    c15->cd();
-    production->SetFillColor(42);
-    production->SetXTitle("z (m)");
-    production->SetYTitle("R (m)");
-    production->Draw();
-
-    TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
-    c10->Divide(2,2);
-    c10->cd(1);
-    pionptspectravertex->SetFillColor(5);
-    pionptspectravertex->SetXTitle("Pt (GeV)");
-    pionptspectravertex->Draw();
-    c10->cd(2);
-    pionptspectrafinal->SetFillColor(5);
-    pionptspectrafinal->SetXTitle("Pt (GeV)");
-    pionptspectrafinal->Draw();
-    c10->cd(3);
-    kaonptspectravertex->SetFillColor(5);
-    kaonptspectravertex->SetXTitle("Pt (GeV)");
-    kaonptspectravertex->Draw();
-    c10->cd(4);
-    kaonptspectrafinal->SetFillColor(5);
-    kaonptspectrafinal->SetXTitle("Pt (GeV)");
-    kaonptspectrafinal->Draw();
-   
-  
-   TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
-   c16->Divide(2,1);
-   
-   c16->cd(1);
-   //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
-   electronspectra1->SetFillColor(5);
-   electronspectra1->SetXTitle("log(GeV)");
-   electronspectra2->SetFillColor(46);
-   electronspectra2->SetXTitle("log(GeV)");
-   electronspectra3->SetFillColor(10);
-   electronspectra3->SetXTitle("log(GeV)");
-   //c13->SetLogx();
-   electronspectra1->Draw();
-   electronspectra2->Draw("same");
-   electronspectra3->Draw("same");
-   
-   c16->cd(2);
-   //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
-   muonspectra1->SetFillColor(5);
-   muonspectra1->SetXTitle("log(GeV)");
-   muonspectra2->SetFillColor(46);
-   muonspectra2->SetXTitle("log(GeV)");
-   muonspectra3->SetFillColor(10);
-   muonspectra3->SetXTitle("log(GeV)");
-   //c14->SetLogx();
-   muonspectra1->Draw();
-   muonspectra2->Draw("same");
-   muonspectra3->Draw("same");
-   
-   //c16->cd(3);
-   //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
-   //neutronspectra1->SetFillColor(42);
-   //neutronspectra1->SetXTitle("log(GeV)");
-   //neutronspectra2->SetFillColor(46);
-   //neutronspectra2->SetXTitle("log(GeV)");
-   //neutronspectra3->SetFillColor(10);
-   //neutronspectra3->SetXTitle("log(GeV)");
-   //c16->SetLogx();
-   //neutronspectra1->Draw();
-   //neutronspectra2->Draw("same");
-   //neutronspectra3->Draw("same");
-
-   TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
-   //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
-   c9->Divide(2,2);
-   
-   c9->cd(1);
-   pionspectra1->SetFillColor(5);
-   pionspectra1->SetXTitle("log(GeV)");
-   pionspectra2->SetFillColor(46);
-   pionspectra2->SetXTitle("log(GeV)");
-   pionspectra3->SetFillColor(10);
-   pionspectra3->SetXTitle("log(GeV)");
-   //c9->SetLogx();
-   pionspectra1->Draw();
-   pionspectra2->Draw("same");
-   pionspectra3->Draw("same");
-   
-   c9->cd(2);
-   //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
-   protonspectra1->SetFillColor(5);
-   protonspectra1->SetXTitle("log(GeV)");
-   protonspectra2->SetFillColor(46);
-   protonspectra2->SetXTitle("log(GeV)");
-   protonspectra3->SetFillColor(10);
-   protonspectra3->SetXTitle("log(GeV)");
-   //c10->SetLogx();
-   protonspectra1->Draw();
-   protonspectra2->Draw("same");
-   protonspectra3->Draw("same");
-   
-   c9->cd(3);
-   //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700); 
-   kaonspectra1->SetFillColor(5);
-   kaonspectra1->SetXTitle("log(GeV)");
-   kaonspectra2->SetFillColor(46);
-   kaonspectra2->SetXTitle("log(GeV)");
-   kaonspectra3->SetFillColor(10);
-   kaonspectra3->SetXTitle("log(GeV)");
-   //c11->SetLogx();
-   kaonspectra1->Draw();
-   kaonspectra2->Draw("same");
-   kaonspectra3->Draw("same");
-   
-   c9->cd(4);
-   //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
-   chargedspectra1->SetFillColor(5);
-   chargedspectra1->SetXTitle("log(GeV)");
-   chargedspectra2->SetFillColor(46);
-   chargedspectra2->SetXTitle("log(GeV)");
-   chargedspectra3->SetFillColor(10);
-   chargedspectra3->SetXTitle("log(GeV)");
-   //c12->SetLogx();
-   chargedspectra1->Draw();
-   chargedspectra2->Draw("same");
-   chargedspectra3->Draw("same");
-   
-
-
-   printf("*****************************************\n");
-   printf("* Particle                   *  Counts  *\n");
-   printf("*****************************************\n");
-
-   printf("* Pions:                     *   %4d   *\n",pion);
-   printf("* Charged Pions:             *   %4d   *\n",chargedpions);
-   printf("* Primary Pions:             *   %4d   *\n",primarypions);
-   printf("* Primary Pions (p>1GeV/c):  *   %4d   *\n",highprimarypions);
-   printf("* Kaons:                     *   %4d   *\n",kaon);
-   printf("* Charged Kaons:             *   %4d   *\n",chargedkaons);
-   printf("* Primary Kaons:             *   %4d   *\n",primarykaons);
-   printf("* Primary Kaons (p>1GeV/c):  *   %4d   *\n",highprimarykaons);
-   printf("* Muons:                     *   %4d   *\n",muon);
-   printf("* Electrons:                 *   %4d   *\n",electron);
-   printf("* Positrons:                 *   %4d   *\n",positron);
-   printf("* Protons:                   *   %4d   *\n",proton);
-   printf("* All Charged:               *   %4d   *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
-   printf("*****************************************\n");
-   //printf("* Photons:                   *   %3.1f   *\n",photons); 
-   //printf("* Primary Photons:           *   %3.1f   *\n",primaryphotons);
-   //printf("* Primary Photons (p>1MeV/c):*   %3.1f   *\n",highprimaryphotons);
-   //printf("*****************************************\n");
-   //printf("* Neutrons:                  *   %3.1f   *\n",neutron);
-   //printf("* Neutrons (p>100keV/c):     *   %3.1f   *\n",highneutrons);
-   //printf("*****************************************\n");
-
-   if (gAlice->TreeD())
-     {
-       gAlice->TreeD()->GetEvent(0);
-   
-       Float_t occ[7]; 
-       Float_t sum=0;
-       Float_t mean=0; 
-       printf("\n*****************************************\n");
-       printf("* Chamber   * Digits      * Occupancy   *\n");
-       printf("*****************************************\n");
-       
-       for (Int_t ich=0;ich<7;ich++)
-        {
-          TClonesArray *Digits = DigitsAddress(ich);    //  Raw clusters branch
-          Int_t ndigits = Digits->GetEntriesFast();
-          occ[ich] = Float_t(ndigits)/(160*144);
-          sum += Float_t(ndigits)/(160*144);
-          printf("*   %d      *    %d      *   %3.1f%%     *\n",ich,ndigits,occ[ich]*100);
-        }
-       mean = sum/7;
-       printf("*****************************************\n");
-       printf("* Mean occupancy          *   %3.1f%%     *\n",mean*100);
-       printf("*****************************************\n");
-     }
-  printf("\nEnd of analysis\n");
-   
-}//void AliRICH::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
-//______________________________________________________________________________
-void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
-{
-
-AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
-   AliRICHSegmentationV0*  segmentation;
-   AliRICHChamber*       chamber;
-   
-   chamber = &(pRICH->Chamber(0));
-   segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel();
-
-   Int_t NpadX = segmentation->Npx();                 // number of pads on X
-   Int_t NpadY = segmentation->Npy();                 // number of pads on Y
-    
-   //Int_t Pad[144][160];
-   /*for (Int_t i=0;i<NpadX;i++) {
-     for (Int_t j=0;j<NpadY;j++) {
-       Pad[i][j]=0;
-     }
-   } */
-
-
-   Int_t xmin= -NpadX/2;  
-   Int_t xmax=  NpadX/2;
-   Int_t ymin= -NpadY/2;
-   Int_t ymax=  NpadY/2;
-
-   Float_t PTfinal = 0;
-   Int_t pionCount = 0;
-   Int_t kaonCount = 0;
-   Int_t protonCount = 0;
-   
-   TH2F *feedback = 0;
-   TH2F *mip = 0;
-   TH2F *cerenkov = 0;
-   TH2F *h = 0;
-   TH1F *hitsX = 0;
-   TH1F *hitsY = 0;
-
-   TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-25,25,150,-45,5);
-
-   if (diaglevel == 1)
-     {
-       printf("Single Ring Hits\n");
-       feedback = new TH2F("feedback","Feedback hit distribution",150,-20,20,150,-35,5);
-       mip = new TH2F("mip","Mip hit distribution",150,-20,20,150,-35,5);
-       cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-20,20,150,-35,5);
-       h = new TH2F("h","Detector hit distribution",150,-20,20,150,-35,5);
-       hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-50,50);
-       hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,50);
-     }       
-   else
-     {
-       printf("Full Event Hits\n");
-       
-       feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
-       mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
-       cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
-       h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300); 
-       hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
-       hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
-     }
-   
-
-
-   TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-   TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
-      
-   TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
-   TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",100,.35,.8);
-   TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
-   TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
-   TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
-   TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
-   TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
-   TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
-   TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
-   //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
-   TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
-   TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
-   TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
-   TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
-   TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
-   TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
-   TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
-   TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
-   TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
-   TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
-   TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
-   TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",50,0,360);
-   TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of theta angle of incidence",50,0,15);
-   TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",50,.5,1);
-   TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",100,0,15);
-   TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",100,0,360);
-   TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",100,.35,.8);
-   TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",100,.35,.8);
-   TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
-   TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
-   TH2F *identification = new TH2F("identification","Particle Identification",100,1,5,100,0,.8);
-   TH1F *OriginalOmega = new TH1F("Original Omega","Cerenkov angle per track",100,.35,.8);
-   TH1F *OriginalPhi = new TH1F("Original Phi","Distribution of phi angle of incidence per track",100,0,360);
-   TH1F *OriginalTheta = new TH1F("Original Theta","Distribution of theta angle per track",100,0,15);
-   TH1F *OmegaError = new TH1F("Omega Error","Difference between original an reconstructed cerenkov angle",100,0,.2);
-   TH1F *PhiError = new TH1F("Phi Error","Difference between original an reconstructed phi angle",100,0,360);
-   TH1F *ThetaError = new TH1F("Theta Error","Difference between original an reconstructed phi angle",100,0,15);
-
-
-//   Start loop over events 
-
-   Int_t Nh=0;
-   Int_t pads=0;
-   Int_t Nh1=0;
-   Int_t mothers[80000];
-   Int_t mothers2[80000];
-   Float_t mom[3];
-   Int_t nraw=0;
-   Int_t phot=0;
-   Int_t feed=0;
-   Int_t padmip=0;
-   Float_t x=0,y=0;
-
-   Float_t chiSquareOmega = 0;
-   Float_t chiSquareTheta = 0;
-   Float_t chiSquarePhi = 0;
-
-   Float_t recEffEvent = 0;
-   Float_t recEffTotal = 0;
-
-   Float_t trackglob[3];
-   Float_t trackloc[3];
-
-   
-   for (Int_t i=0;i<100;i++) mothers[i]=0;
-
-   for (int nev=0; nev<= evNumber2; nev++) {
-       Int_t nparticles = gAlice->GetEvent(nev);
-       
-
-       //cout<<"nev  "<<nev<<endl;
-       printf ("\n**********************************\nProcessing Event: %d\n",nev);
-       //cout<<"nparticles  "<<nparticles<<endl;
-       printf ("Particles       : %d\n\n",nparticles);
-       if (nev < evNumber1) continue;
-       if (nparticles <= 0) return;
-       
-// Get pointers to RICH detector and Hits containers
-       
-
-       TTree *TH = TreeH(); 
-       Stat_t ntracks = TH->GetEntries();
-
-       // Start loop on tracks in the hits containers
-       //Int_t Nc=0;
-       for (Int_t track=0; track<ntracks;track++) {
-          
-        printf ("\nProcessing Track: %d\n",track);
-        gAlice->ResetHits();
-        TH->GetEvent(track);
-        Int_t nhits = pRICH->Hits()->GetEntriesFast();
-        if (nhits) Nh+=nhits;
-        printf("Hits            : %d\n",nhits);
-        for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1); 
-            mHit;
-            mHit=(AliRICHhit*)pRICH->NextHit()) 
-          {
-            Int_t nch  = mHit->Chamber();              // chamber number
-            trackglob[0] = mHit->X();                 // x-pos of hit
-            trackglob[1] = mHit->Y();
-            trackglob[2] = mHit->Z();                 // y-pos of hit
-            //x  = mHit->X();                           // x-pos of hit
-            //y  = mHit->Z();                           // y-pos
-            Float_t phi = mHit->Phi();                 //Phi angle of incidence
-            Float_t theta = mHit->Theta();             //Theta angle of incidence
-            Int_t index = mHit->Track();
-            Int_t particle = (Int_t)(mHit->Particle());        
-            //Int_t freon = (Int_t)(mHit->fLoss);    
-            Float_t px = mHit->MomX();
-            Float_t py = mHit->MomY();
-            
-            if (TMath::Abs(particle) < 10000000)
-              {
-                PTfinal=TMath::Sqrt(px*px + py*py);
-              }
-       
-            chamber = &(pRICH->Chamber(nch-1));
-            
-            
-            chamber->GlobaltoLocal(trackglob,trackloc);
-            
-            chamber->LocaltoGlobal(trackloc,trackglob);
-            
-       
-            x=trackloc[0];
-            y=trackloc[2];
-            
-            hitsX->Fill(x,(float) 1);
-            hitsY->Fill(y,(float) 1);
-              
-             
-             TParticle *current = (TParticle*)gAlice->Particle(index);
-             //printf("Particle type: %d\n",sizeoff(Particles));
-
-             hitsTheta->Fill(theta,(float) 1);
-             //hitsPhi->Fill(phi,(float) 1);
-             //if (pRICH->GetDebugLevel() == -1)
-            
-             if (current->GetPdgCode() < 10000000)
-               {
-                 mip->Fill(x,y,(float) 1);
-                 //printf("adding mip\n");
-                 //if (current->Energy() - current->GetCalcMass()>1 && freon==1)
-                 //{
-                 hitsPhi->Fill(TMath::Abs(phi),(float) 1);
-                 //hitsTheta->Fill(theta,(float) 1);
-                 //printf("Theta:%f, Phi:%f\n",theta,phi);
-                 //}
-               }
-             
-             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
-               {
-                 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
-               }
-             if (TMath::Abs(particle)==2212)
-               {
-                 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
-               }
-             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
-                 || TMath::Abs(particle)==311)
-               {
-                 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
-               }
-             if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
-               {
-                 if (current->Energy() - current->GetCalcMass()>1)
-                   chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
-               }
-             //printf("Hits:%d\n",hit);
-             //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
-              // Fill the histograms
-             Nh1+=nhits;
-             h->Fill(x,y,(float) 1);
-                 //}
-              //}
-          }
-          
-          Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
-          //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
-          //totalphotonsevent->Fill(ncerenkovs,(float) 1);
-
-          if (ncerenkovs) {
-            printf("Cerenkovs       : %d\n",ncerenkovs);
-            totalphotonsevent->Fill(ncerenkovs,(float) 1);
-            for (Int_t hit=0;hit<ncerenkovs;hit++) {
-              AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
-              Int_t nchamber = cHit->fChamber;     // chamber number
-              Int_t index =    cHit->Track();
-              //Int_t pindex =   (Int_t)(cHit->fIndex);
-              trackglob[0] = cHit->X();                 // x-pos of hit
-              trackglob[1] = cHit->Y();
-              trackglob[2] = cHit->Z();                 // y-pos of hit
-              //Float_t cx  =      cHit->X();                // x-position
-              //Float_t cy  =      cHit->Z();                // y-position
-              Int_t cmother =  cHit->fCMother;      // Index of mother particle
-              Int_t closs =    (Int_t)(cHit->fLoss);           // How did the particle get lost? 
-              Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
-              
-              chamber = &(pRICH->Chamber(nchamber-1));
-            
-              //printf("Nch:%d\n",nch);
-              
-              chamber->GlobaltoLocal(trackglob,trackloc);
-            
-              chamber->LocaltoGlobal(trackloc,trackglob);
-            
-       
-              Float_t cx=trackloc[0];
-              Float_t cy=trackloc[2];
-              
-              //printf ("Cerenkov hit number %d/%d, X:%f, Y:%f\n",hit,ncerenkovs,cx,cy); 
-
-
-              //printf("Particle:%9d\n",index);
-                                
-              TParticle *current = (TParticle*)gAlice->Particle(index);
-              Float_t energyckov = current->Energy();
-              
-              if (current->GetPdgCode() == 50000051)
-                {
-                  if (closs==4)
-                    {
-                      feedback->Fill(cx,cy,(float) 1);
-                      feed++;
-                    }
-                }
-              if (current->GetPdgCode() == 50000050)
-                {
-                  
-                  if (closs !=4)
-                    {
-                      phspectra2->Fill(energyckov*1e9,(float) 1);
-                    }
-                      
-                  if (closs==4)
-                    {
-                      cerenkov->Fill(cx,cy,(float) 1); 
-                      
-                      //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy); 
-                      
-                      //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
-                      AliRICHhit* mipHit = (AliRICHhit*) pRICH->Hits()->UncheckedAt(0);
-                      mom[0] = current->Px();
-                      mom[1] = current->Py();
-                      mom[2] = current->Pz();
-                      //mom[0] = cHit->fMomX;
-                      // mom[1] = cHit->fMomZ;
-                      //mom[2] = cHit->fMomY;
-                      //Float_t energymip = MIP->Energy();
-                      //Float_t Mip_px = mipHit->fMomFreoX;
-                      //Float_t Mip_py = mipHit->fMomFreoY;
-                      //Float_t Mip_pz = mipHit->fMomFreoZ;
-                      //Float_t Mip_px = MIP->Px();
-                      //Float_t Mip_py = MIP->Py();
-                      //Float_t Mip_pz = MIP->Pz();
-                      
-                      
-                      
-                      //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
-                      //Float_t rt = TMath::Sqrt(r);
-                      //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz; 
-                      //Float_t Mip_rt = TMath::Sqrt(Mip_r);
-                      //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
-                      //Float_t cherenkov = TMath::ACos(coscerenkov);
-                      ckovangle->Fill(cherenkov,(float) 1);                           //Cerenkov angle calculus
-                      //printf("Cherenkov: %f\n",cherenkov);
-                      Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
-                      hckphi->Fill(ckphi,(float) 1);
-                      
-                      
-                      //Float_t mix = MIP->Vx();
-                      //Float_t miy = MIP->Vy();
-                      Float_t mx = mipHit->X();
-                      Float_t my = mipHit->Z();
-                      //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
-                      Float_t dx = trackglob[0] - mx;
-                      Float_t dy = trackglob[2] - my;
-                      //printf("Dx:%f, Dy:%f\n",dx,dy);
-                      Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
-                      //printf("Final radius:%f\n",final_radius);
-                      radius->Fill(final_radius,(float) 1);
-                      
-                      phspectra1->Fill(energyckov*1e9,(float) 1);
-                      phot++;
-                    }
-                  for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
-                    if (cmother == nmothers){
-                      if (closs == 4)
-                        mothers2[cmother]++;
-                      mothers[cmother]++;
-                    }
-                  } 
-                }
-            }
-          }
-          
-
-          if(gAlice->TreeR())
-            {
-              Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
-              gAlice->TreeR()->GetEvent(nent-1);
-              TClonesArray *Rawclusters = pRICH->RawClustAddress(2);    //  Raw clusters branch
-              //printf ("Rawclusters:%p",Rawclusters);
-              Int_t nrawclusters = Rawclusters->GetEntriesFast();
-                      
-              if (nrawclusters) {
-                printf("Raw Clusters    : %d\n",nrawclusters);
-                for (Int_t hit=0;hit<nrawclusters;hit++) {
-                  AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
-                  //Int_t nchamber = rcHit->fChamber;     // chamber number
-                  //Int_t nhit = cHit->fHitNumber;        // hit number
-                  Int_t qtot = rcHit->fQ;                 // charge
-                  Float_t fx  =  rcHit->fX;                 // x-position
-                  Float_t fy  =  rcHit->fY;                 // y-position
-                  //Int_t type = rcHit->fCtype;             // cluster type ?   
-                  Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
-                  pads += mult;
-                  if (qtot > 0) {
-                    //printf ("fx: %d, fy: %d\n",fx,fy);
-                    if (fx>(x-4) && fx<(x+4)  && fy>(y-4) && fy<(y+4)) {
-                      //printf("There %d \n",mult);
-                      padmip+=mult;
-                    } else {
-                      padnumber->Fill(mult,(float) 1);
-                      nraw++;
-                      if (mult<4) Clcharge->Fill(qtot,(float) 1);
-                    }
-                    
-                  }
-                }
-              }
-              
-              
-              TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
-              Int_t nrechits1D = RecHits1D->GetEntriesFast();
-              //printf (" nrechits:%d\n",nrechits);
-              
-              if(nrechits1D)
-                {
-                  for (Int_t hit=0;hit<nrechits1D;hit++) {
-                    AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
-                    Float_t r_omega = recHit1D->fOmega;                  // Cerenkov angle
-                    Float_t *cer_pho = recHit1D->fCerPerPhoton;        // Cerenkov angle per photon
-                    Int_t *padsx = recHit1D->fPadsUsedX;           // Pads Used fo reconstruction (x)
-                    Int_t *padsy = recHit1D->fPadsUsedY;           // Pads Used fo reconstruction (y)
-                    Int_t goodPhotons = recHit1D->fGoodPhotons;    // Number of pads used for reconstruction
-                    
-                    Omega1D->Fill(r_omega,(float) 1);
-                    
-                    for (Int_t i=0; i<goodPhotons; i++)
-                      {
-                        PhotonCer->Fill(cer_pho[i],(float) 1);
-                        PadsUsed->Fill(padsx[i],padsy[i],1);
-                        //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
-                      }
-                    
-                    //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
-                  }
-                }
-
-              
-              TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
-              Int_t nrechits3D = RecHits3D->GetEntriesFast();
-              //printf (" nrechits:%d\n",nrechits);
-              
-              if(nrechits3D)
-                {
-                  recEffEvent = 0;
-                  
-                  //for (Int_t hit=0;hit<nrechits3D;hit++) {
-                  AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(track);
-                  Float_t r_omega    = recHit3D->fOmega;                  // Cerenkov angle
-                  Float_t r_theta    = recHit3D->fTheta;                  // Theta angle of incidence
-                  Float_t r_phi      = recHit3D->fPhi;                    // Phi angle if incidence
-                  Float_t meanradius = recHit3D->fMeanRadius;              // Mean radius for reconstructed point
-                  Float_t originalOmega = recHit3D->fOriginalOmega;       // Real Cerenkov angle
-                  Float_t originalTheta = recHit3D->fOriginalTheta;       // Real incidence angle
-                  Float_t originalPhi = recHit3D->fOriginalPhi;           // Real azimuthal angle
-                  
-                  
-                  //correction to track cerenkov angle
-                  originalOmega = (Float_t) ckovangle->GetMean();
-                  
-                  if(diaglevel == 4)
-                    {
-                      printf("\nMean cerenkov angle: %f\n", originalOmega);
-                      printf("Reconstructed cerenkov angle: %f\n",r_omega);
-                    }
-                  
-                  Float_t omegaError = TMath::Abs(originalOmega - r_omega);
-                  Float_t thetaError = TMath::Abs(originalTheta - r_theta);
-                  Float_t phiError   = TMath::Abs(originalPhi - r_phi);
-                  
-                  //chiSquareOmega += (omegaError/originalOmega)*(omegaError/originalOmega); 
-                  //chiSquareTheta += (thetaError/originalTheta)*(thetaError/originalTheta); 
-                  //chiSquarePhi += (phiError/originalPhi)*(phiError/originalPhi); 
-                  
-                  if(TMath::Abs(omegaError) < 0.015)
-                    recEffEvent += 1;
-                  
-                  
-                  
-                  //printf("rechit %f %f %f %f %f\n",recHit3D->fOmega,recHit3D->fTheta,recHit3D->fPhi, recHit3D->fX,recHit3D->fY);  
-                  
-                  Omega3D->Fill(r_omega,(float) 1);
-                  Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
-                  Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
-                  MeanRadius->Fill(meanradius,(float) 1);
-                  identification->Fill(PTfinal, r_omega,1);
-                  OriginalOmega->Fill(originalOmega, (float) 1);
-                  OriginalTheta->Fill(originalTheta, (float) 1);
-                  OriginalPhi->Fill(TMath::Abs(originalPhi), (float) 1);
-                  OmegaError->Fill(omegaError, (float) 1);
-                  ThetaError->Fill(thetaError, (float) 1);
-                  PhiError->Fill(phiError, (float) 1);
-                  
-                  recEffEvent = recEffEvent;
-                  recEffTotal += recEffEvent;
-                  
-                  Float_t pioncer = acos(sqrt((.139*.139+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
-                  Float_t kaoncer = acos(sqrt((.439*.439+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
-                  Float_t protoncer = acos(sqrt((.938*.938+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
-
-                  Float_t piondist = TMath::Abs(r_omega - pioncer);
-                  Float_t kaondist = TMath::Abs(r_omega - kaoncer);
-                  Float_t protondist = TMath::Abs(r_omega - protoncer);
-
-                  if(diaglevel == 4)
-                    {
-                      if(pioncer<r_omega)
-                        {
-                          printf("Identified as a PION!\n");
-                          pionCount += 1;
-                        }
-                      if(kaoncer<r_omega && pioncer>r_omega)
-                        {
-                          if(kaondist>piondist)
-                            {
-                              printf("Identified as a PION!\n");
-                              pionCount += 1;
-                            }
-                          else
-                            {
-                              printf("Identified as a KAON!\n");
-                              kaonCount += 1;
-                            }
-                        }                       }
-                      if(protoncer<r_omega && kaoncer>r_omega)
-                        {
-                          if(kaondist>protondist)
-                            {
-                              printf("Identified as a PROTON!\n");
-                              protonCount += 1;
-                            }
-                          else
-                            {
-                              printf("Identified as a KAON!\n");
-                              pionCount += 1;
-                            }
-                        }
-                      if(protoncer>r_omega)
-                        {
-                          printf("Identified as a PROTON!\n");
-                          protonCount += 1;
-                        }
-
-                      printf("\nReconstruction efficiency: %5.2f%%\n", recEffEvent*100);
-                }
-            }
-       }
-   
-       
-       for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
-        totalphotonstrack->Fill(mothers[nmothers],(float) 1);
-        mother->Fill(mothers2[nmothers],(float) 1);
-        //printf ("Entries in %d : %d\n",nmothers, mothers[nmothers]);
-       }
-       
-       clusev->Fill(nraw,(float) 1);
-       photev->Fill(phot,(float) 1);
-       feedev->Fill(feed,(float) 1);
-       padsmip->Fill(padmip,(float) 1);
-       padscl->Fill(pads,(float) 1);
-       //printf("Photons:%d\n",phot);
-       phot = 0;
-       feed = 0;
-       pads = 0;
-       nraw=0;
-       padmip=0;
-       
-       
-       
-       gAlice->ResetDigits();
-       //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
-       gAlice->TreeD()->GetEvent(0);
-       
-       if (diaglevel < 4)
-        {
-          
-          
-          TClonesArray *Digits  = pRICH->DigitsAddress(2);
-          Int_t ndigits = Digits->GetEntriesFast();
-          printf("Digits          : %d\n",ndigits);
-          padsev->Fill(ndigits,(float) 1);
-          for (Int_t hit=0;hit<ndigits;hit++) {
-            AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
-            Int_t qtot = dHit->Signal();                // charge
-            Int_t ipx  = dHit->PadX();               // pad number on X
-            Int_t ipy  = dHit->PadY();               // pad number on Y
-            //printf("%d, %d\n",ipx,ipy);
-            if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
-          }
-        }
-       
-       if (diaglevel == 5)
-        {
-          for (Int_t ich=0;ich<7;ich++)
-            {
-              TClonesArray *Digits = pRICH->DigitsAddress(ich);    //  Raw clusters branch
-              Int_t ndigits = Digits->GetEntriesFast();
-              //printf("Digits:%d\n",ndigits);
-              padsev->Fill(ndigits,(float) 1); 
-              if (ndigits) {
-                for (Int_t hit=0;hit<ndigits;hit++) {
-                  AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
-                  //Int_t nchamber = dHit->GetChamber();     // chamber number
-                  //Int_t nhit = dHit->fHitNumber;          // hit number
-                  Int_t qtot = dHit->Signal();                // charge
-                  Int_t ipx  = dHit->PadX();               // pad number on X
-                  Int_t ipy  = dHit->PadY();               // pad number on Y
-                  //Int_t iqpad  = dHit->fQpad;           // charge per pad
-                  //Int_t rpad  = dHit->fRSec;            // R-position of pad
-                  //printf ("Pad hit, PadX:%d, PadY:%d\n",ipx,ipy);
-                  if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
-                  if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
-                }
-              }
-            }
-        }
-   }
-   
-   if(diaglevel == 4)
-     {
-
-       Stat_t omegaE;
-       Stat_t thetaE;
-       Stat_t phiE;
-       
-       Stat_t omegaO;
-       Stat_t thetaO;
-       Stat_t phiO;
-       
-       for(Int_t i=0;i<99;i++)
-        {
-          omegaE = OriginalOmega->GetBinContent(i);
-          if(omegaE != 0)
-            {
-              omegaO = Omega3D->GetBinContent(i);
-              chiSquareOmega += (TMath::Power(omegaE,2) - TMath::Power(omegaO,2))/omegaO;
-            }
-
-          thetaE = OriginalTheta->GetBinContent(i);
-          if(thetaE != 0)
-            {
-              thetaO = Theta->GetBinContent(i);
-              chiSquareTheta += (TMath::Power(thetaE,2) - TMath::Power(thetaO,2))/thetaO;
-            }
-
-          phiE = OriginalPhi->GetBinContent(i);
-          if(phiE != 0)
-            {
-              phiO = Phi->GetBinContent(i);
-              chiSquarePhi += (TMath::Power(phiE,2) - TMath::Power(phiO,2))/phiO;
-            }
-          
-          //printf(" o: %f  t: %f  p: %f\n", OriginalOmega->GetBinContent(i), OriginalTheta->GetBinContent(i),OriginalPhi->GetBinContent(i));
-
-        }
-
-       
-
-       printf("\nChi square test values:   Omega - %f\n", chiSquareOmega);
-       printf("                          Theta - %f\n", chiSquareTheta);
-       printf("                          Phi   - %f\n", chiSquarePhi);
-       
-       printf("\nKolmogorov test values:   Omega - %5.4f\n", Omega3D->KolmogorovTest(OriginalOmega));
-       printf("                          Theta - %5.4f\n", Theta->KolmogorovTest(OriginalTheta));
-       printf("                          Phi   - %5.4f\n", Phi->KolmogorovTest(OriginalPhi));
-
-       recEffTotal = recEffTotal/evNumber2;
-       printf("\nTotal reconstruction efficiency: %5.2f%%\n", recEffTotal*100);
-       printf("\n Pions: %d\n Kaons: %d\n Protons:%d\n",pionCount, kaonCount, protonCount);
-
-     }
-   
-   
-   //Create canvases, set the view range, show histograms
-
-   TCanvas *c1 = 0;
-   TCanvas *c2 = 0;
-   TCanvas *c3 = 0;
-   TCanvas *c4 = 0;
-   TCanvas *c5 = 0;
-   TCanvas *c6 = 0;
-   TCanvas *c7 = 0;
-   TCanvas *c8 = 0;
-   TCanvas *c9 = 0;
-   TCanvas *c10 = 0;
-   TCanvas *c11 = 0;
-   TCanvas *c12 = 0;
-   TCanvas *c13 = 0;
-
-   //TF1* expo = 0;
-   //TF1* gaus = 0;
-   
-   TStyle *mystyle=new TStyle("Plain","mystyle");
-   mystyle->SetPalette(1,0);
-   //mystyle->SetTitleYSize(0.2);
-   //mystyle->SetStatW(0.19);
-   //mystyle->SetStatH(0.1);
-   //mystyle->SetStatFontSize(0.01);
-   //mystyle->SetTitleYSize(0.3);
-   mystyle->SetFuncColor(2);
-   //mystyle->SetOptStat(0111);
-   mystyle->SetDrawBorder(0);
-   mystyle->SetTitleBorderSize(0);
-   mystyle->SetOptFit(1111);
-   mystyle->cd();
-
-   
-   TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
-   Int_t nrechits3D = RecHits3D->GetEntriesFast();
-   TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
-   Int_t nrechits1D = RecHits1D->GetEntriesFast();
-
-  switch(diaglevel)
-     {
-     case 1:
-       
-       c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
-       hc0->SetXTitle("ix (npads)");
-       hc0->Draw("colz");
-       
-//
-       c2 = new TCanvas("c2","Hits per type",100,100,600,700);
-       c2->Divide(2,2);
-       //c4->SetFillColor(42);
-
-       c2->cd(1);
-       feedback->SetXTitle("x (cm)");
-       feedback->SetYTitle("y (cm)");
-       feedback->Draw("colz");
-       
-       c2->cd(2);
-       //mip->SetFillColor(5);
-       mip->SetXTitle("x (cm)");
-       mip->SetYTitle("y (cm)");
-       mip->Draw("colz");
-       
-       c2->cd(3);
-       //cerenkov->SetFillColor(5);
-       cerenkov->SetXTitle("x (cm)");
-       cerenkov->SetYTitle("y (cm)"); 
-       cerenkov->Draw("colz");
-       
-       c2->cd(4);
-       //h->SetFillColor(5);
-       h->SetXTitle("x (cm)");
-       h->SetYTitle("y (cm)");
-       h->Draw("colz");
-
-       c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
-       c3->Divide(2,1);
-       //c10->SetFillColor(42);
-       
-       c3->cd(1);
-       hitsX->SetFillColor(5);
-       hitsX->SetXTitle("(cm)");
-       hitsX->Draw();
-       
-       c3->cd(2);
-       hitsY->SetFillColor(5);
-       hitsY->SetXTitle("(cm)");
-       hitsY->Draw();
-       
-      
-       break;
-//
-     case 2:
-       
-       c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
-       c4->Divide(2,1);
-       //c6->SetFillColor(42);
-       
-       c4->cd(1);
-       phspectra2->SetFillColor(5);
-       phspectra2->SetXTitle("energy (eV)");
-       phspectra2->Draw();
-       c4->cd(2);
-       phspectra1->SetFillColor(5);
-       phspectra1->SetXTitle("energy (eV)");
-       phspectra1->Draw();
-       
-       c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
-       c5->Divide(2,2);
-       //c9->SetFillColor(42);
-       
-       c5->cd(1);
-       pionspectra->SetFillColor(5);
-       pionspectra->SetXTitle("(GeV)");
-       pionspectra->Draw();
-       
-       c5->cd(2);
-       protonspectra->SetFillColor(5);
-       protonspectra->SetXTitle("(GeV)");
-       protonspectra->Draw();
-       
-       c5->cd(3);
-       kaonspectra->SetFillColor(5);
-       kaonspectra->SetXTitle("(GeV)");
-       kaonspectra->Draw();
-       
-       c5->cd(4);
-       chargedspectra->SetFillColor(5);
-       chargedspectra->SetXTitle("(GeV)");
-       chargedspectra->Draw();
-
-       break;
-       
-     case 3:
-
-       
-       if(gAlice->TreeR())
-        {
-          c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
-          c6->Divide(2,2);
-          //c3->SetFillColor(42);
-          
-          c6->cd(1);
-          //TPad* c6_1;
-          //c6_1->SetLogy();
-          Clcharge->SetFillColor(5);
-          Clcharge->SetXTitle("ADC counts");
-          if (evNumber2>10)
-            {
-              Clcharge->Fit("expo");
-              //expo->SetLineColor(2);
-              //expo->SetLineWidth(3);
-            }
-          Clcharge->Draw();
-          
-          c6->cd(2);
-          padnumber->SetFillColor(5);
-          padnumber->SetXTitle("(counts)");
-          padnumber->Draw();
-          
-          c6->cd(3);
-          clusev->SetFillColor(5);
-          clusev->SetXTitle("(counts)");
-          if (evNumber2>10)
-            {
-              clusev->Fit("gaus");
-              //gaus->SetLineColor(2);
-              //gaus->SetLineWidth(3);
-            }
-          clusev->Draw();
-          
-          c6->cd(4);
-          padsmip->SetFillColor(5);
-          padsmip->SetXTitle("(counts)");
-          padsmip->Draw(); 
-        }
-       
-       if(evNumber2<1)
-        {
-          c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
-          mother->SetFillColor(5);
-          mother->SetXTitle("counts");
-          mother->Draw();
-        }
-
-       c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
-       c7->Divide(2,2);
-       //c7->SetFillColor(42);
-       
-       c7->cd(1);
-       totalphotonsevent->SetFillColor(5);
-       totalphotonsevent->SetXTitle("Photons (counts)");
-       if (evNumber2>10)
-          {
-            totalphotonsevent->Fit("gaus");
-            //gaus->SetLineColor(2);
-            //gaus->SetLineWidth(3);
-          }
-       totalphotonsevent->Draw();
-       
-       c7->cd(2);
-       photev->SetFillColor(5);
-       photev->SetXTitle("(counts)");
-       if (evNumber2>10)
-        {
-          photev->Fit("gaus");
-          //gaus->SetLineColor(2);
-          //gaus->SetLineWidth(3);
-        }
-       photev->Draw();
-       
-       c7->cd(3);
-       feedev->SetFillColor(5);
-       feedev->SetXTitle("(counts)");
-       if (evNumber2>10)
-        {
-          feedev->Fit("gaus");
-          //gaus->SetLineColor(2);
-          //gaus->SetLineWidth(3);
-        }
-       feedev->Draw();
-
-       c7->cd(4);
-       padsev->SetFillColor(5);
-       padsev->SetXTitle("(counts)");
-       if (evNumber2>10)
-        {
-          padsev->Fit("gaus");
-          //gaus->SetLineColor(2);
-          //gaus->SetLineWidth(3);
-        }
-       padsev->Draw();
-
-       break;
-
-     case 4:
-       
-
-       if(nrechits3D)
-        {
-          c8 = new TCanvas("c8","3D reconstruction of Phi angle",50,50,300,1050);
-          c8->Divide(1,3);
-          //c2->SetFillColor(42);
-          
-          
-          // data per hit
-          c8->cd(1);
-          hitsPhi->SetFillColor(5);
-          if (evNumber2>10)
-            hitsPhi->Fit("gaus");
-          hitsPhi->Draw();
-          
-           //data per track
-          c8->cd(2);
-          OriginalPhi->SetFillColor(5);
-          if (evNumber2>10)
-            OriginalPhi->Fit("gaus");
-          OriginalPhi->Draw();
-
-          //recontructed data
-          c8->cd(3);
-          Phi->SetFillColor(5);
-          if (evNumber2>10)
-            Phi->Fit("gaus");
-          Phi->Draw();
-
-          c9 = new TCanvas("c9","3D reconstruction of theta angle",75,75,300,1050);
-          c9->Divide(1,3);
-
-          // data per hit
-          c9->cd(1);
-          hitsTheta->SetFillColor(5);
-          if (evNumber2>10)
-            hitsTheta->Fit("gaus");
-          hitsTheta->Draw();
-          
-          //data per track
-          c9->cd(2);
-          OriginalTheta->SetFillColor(5);
-          if (evNumber2>10)
-            OriginalTheta->Fit("gaus");
-          OriginalTheta->Draw();
-
-          //recontructed data
-          c9->cd(3);
-          Theta->SetFillColor(5);
-          if (evNumber2>10)
-            Theta->Fit("gaus");
-          Theta->Draw();
-
-          c10 = new TCanvas("c10","3D reconstruction of cherenkov angle",100,100,300,1050);
-          c10->Divide(1,3);
-
-          // data per hit
-          c10->cd(1);
-          ckovangle->SetFillColor(5);
-          ckovangle->SetXTitle("angle (radians)");
-          if (evNumber2>10)
-            ckovangle->Fit("gaus");
-          ckovangle->Draw();
-          
-          //data per track
-          c10->cd(2);
-          OriginalOmega->SetFillColor(5);
-          OriginalOmega->SetXTitle("angle (radians)");
-          if (evNumber2>10)
-            OriginalOmega->Fit("gaus");
-          OriginalOmega->Draw();
-
-          //recontructed data
-          c10->cd(3);
-          Omega3D->SetFillColor(5);
-          Omega3D->SetXTitle("angle (radians)");
-          if (evNumber2>10)
-            Omega3D->Fit("gaus");
-          Omega3D->Draw(); 
-
-
-          c11 = new TCanvas("c11","3D reconstruction of mean radius",125,125,300,700);
-          c11->Divide(1,2);
-
-          // data per hit
-          c11->cd(1);
-          radius->SetFillColor(5);
-          radius->SetXTitle("radius (cm)");
-          radius->Draw();
-
-          //recontructed data
-          c11->cd(2);
-          MeanRadius->SetFillColor(5);
-          MeanRadius->SetXTitle("radius (cm)");
-          MeanRadius->Draw();
-
-          
-          c12 = new TCanvas("c12","Cerenkov angle vs. Momentum",150,150,550,350);
-
-          c12->cd(1);
-          identification->SetFillColor(5);
-          identification->SetXTitle("Momentum (GeV/c)");
-          identification->SetYTitle("Cherenkov angle (radians)");
-          
-          //Float_t pionmass=.139;
-          //Float_t kaonmass=.493;
-          //Float_t protonmass=.938;
-          //Float_t n=1.295;
-          
-          TF1 *pionplot = new TF1("pion","acos(sqrt((.139*.139+x*x)/(x*x*1.285*1.285)))",1,5);
-          TF1 *kaonplot = new TF1("kaon","acos(sqrt((.439*.439+x*x)/(x*x*1.285*1.285)))",1,5);
-          TF1 *protonplot = new TF1("proton","acos(sqrt((.938*.938+x*x)/(x*x*1.285*1.285)))",1,5);
-          
-          identification->Draw();
-
-          pionplot->SetLineColor(5);
-          pionplot->Draw("same");
-
-          kaonplot->SetLineColor(4);
-          kaonplot->Draw("same");
-
-          protonplot->SetLineColor(3);
-          protonplot->Draw("same");
-          //identification->Draw("same");
-
-
-
-          c13 = new TCanvas("c13","Reconstruction Errors",200,200,900,350);
-          c13->Divide(3,1);
-
-          c13->cd(1);
-          PhiError->SetFillColor(5);
-          if (evNumber2>10)
-            PhiError->Fit("gaus");
-          PhiError->Draw();
-          c13->cd(2);
-          ThetaError->SetFillColor(5);
-          if (evNumber2>10)
-            ThetaError->Fit("gaus");
-          ThetaError->Draw();
-          c13->cd(3);
-          OmegaError->SetFillColor(5);
-          OmegaError->SetXTitle("angle (radians)");
-          if (evNumber2>10)
-            OmegaError->Fit("gaus");
-          OmegaError->Draw();
-          
-        }
-       
-       if(nrechits1D)
-        {
-          c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
-          c9->Divide(3,2);
-          //c5->SetFillColor(42);
-          
-          c9->cd(1);
-          ckovangle->SetFillColor(5);
-          ckovangle->SetXTitle("angle (radians)");
-          ckovangle->Draw();
-          
-          c9->cd(2);
-          radius->SetFillColor(5);
-          radius->SetXTitle("radius (cm)");
-          radius->Draw();
-          
-          c9->cd(3);
-          hc0->SetXTitle("pads");
-          hc0->Draw("box"); 
-          
-          c9->cd(5);
-          Omega1D->SetFillColor(5);
-          Omega1D->SetXTitle("angle (radians)");
-          Omega1D->Draw();
-          
-          c9->cd(4);
-          PhotonCer->SetFillColor(5);
-          PhotonCer->SetXTitle("angle (radians)");
-          PhotonCer->Draw();
-          
-          c9->cd(6);
-          PadsUsed->SetXTitle("pads");
-          PadsUsed->Draw("box"); 
-        }
-       
-       break;
-       
-     case 5:
-       
-       printf("Drawing histograms.../n");
-
-       //if (ndigits)
-        //{
-       c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
-       c1->Divide(4,2);
-       //c1->SetFillColor(42);
-       
-       c10->cd(1);
-       hc1->SetXTitle("ix (npads)");
-       hc1->Draw("box");
-       c10->cd(2);
-       hc2->SetXTitle("ix (npads)");
-       hc2->Draw("box");
-       c10->cd(3);
-       hc3->SetXTitle("ix (npads)");
-       hc3->Draw("box");
-       c10->cd(4);
-       hc4->SetXTitle("ix (npads)");
-       hc4->Draw("box");
-       c10->cd(5);
-       hc5->SetXTitle("ix (npads)");
-       hc5->Draw("box");
-       c10->cd(6);
-       hc6->SetXTitle("ix (npads)");
-       hc6->Draw("box");
-       c10->cd(7);
-       hc7->SetXTitle("ix (npads)");
-       hc7->Draw("box");
-       c10->cd(8);
-       hc0->SetXTitle("ix (npads)");
-       hc0->Draw("box");
-        //}
-//
-       c11 = new TCanvas("c11","Hits per type",100,100,600,700);
-       c11->Divide(2,2);
-       //c4->SetFillColor(42);
-       
-       c11->cd(1);
-       feedback->SetXTitle("x (cm)");
-       feedback->SetYTitle("y (cm)");
-       feedback->Draw();
-       
-       c11->cd(2);
-       //mip->SetFillColor(5);
-       mip->SetXTitle("x (cm)");
-       mip->SetYTitle("y (cm)");
-       mip->Draw();
-       
-       c11->cd(3);
-       //cerenkov->SetFillColor(5);
-       cerenkov->SetXTitle("x (cm)");
-       cerenkov->SetYTitle("y (cm)"); 
-       cerenkov->Draw();
-       
-       c11->cd(4);
-       //h->SetFillColor(5);
-       h->SetXTitle("x (cm)");
-       h->SetYTitle("y (cm)");
-       h->Draw();
-
-       c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
-       c12->Divide(2,1);
-       //c10->SetFillColor(42);
-       
-       c12->cd(1);
-       hitsX->SetFillColor(5);
-       hitsX->SetXTitle("(cm)");
-       hitsX->Draw();
-       
-       c12->cd(2);
-       hitsY->SetFillColor(5);
-       hitsY->SetXTitle("(cm)");
-       hitsY->Draw();
-       
-       break;
-       
-     }
-       
-
-   // calculate the number of pads which give a signal
-
-
-   //Int_t Np=0;
-   /*for (Int_t i=0;i< NpadX;i++) {
-       for (Int_t j=0;j< NpadY;j++) {
-          if (Pad[i][j]>=6){
-              Np+=1;
-          }
-       }
-   }*/
-   //printf("The total number of pads which give a signal: %d %d\n",Nh,Nh1);
-   printf("\nEnd of analysis\n");
-   printf("**********************************\n");
-}//void AliRICH::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::MakeBranchInTreeD(TTree *treeD, const char *file)
 {// Create TreeD branches for the RICH.
   if(GetDebug())Info("MakeBranchInTreeD","Start.");
@@ -2676,7 +896,7 @@ void AliRICH::MakeBranchInTreeD(TTree *treeD, const char *file)
     }
   }
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::MakeBranch(Option_t* option)
 {//Create Tree branches for the RICH.
   if(GetDebug())Info("MakeBranch","Start with option= %s.",option);
@@ -2690,7 +910,6 @@ void AliRICH::MakeBranch(Option_t* option)
   const char *cR = strstr(option,"R");
   const char *cS = strstr(option,"S");
 
-
   if(cH&&TreeH()){
     if(!fHits) fHits=new TClonesArray("AliRICHhit",1000  );
     if(!fCerenkovs) fCerenkovs  = new TClonesArray("AliRICHCerenkov",1000);
@@ -2760,7 +979,7 @@ void AliRICH::MakeBranch(Option_t* option)
      }
    }//if (cR && gAlice->TreeR())
   if(GetDebug())Info("MakeBranch","Stop.");   
-}
+}//void AliRICH::MakeBranch(Option_t* option)
 //______________________________________________________________________________
 void AliRICH::SetTreeAddress()
 {//Set branch address for the Hits and Digits Tree.
@@ -2774,7 +993,6 @@ void AliRICH::SetTreeAddress()
   TTree *treeH = fLoader->TreeH();
   TTree *treeD = fLoader->TreeD();
   TTree *treeR = fLoader->TreeR();
-  TTree *treeS = fLoader->TreeS();
     
   if(treeH){
     if(GetDebug())Info("SetTreeAddress","tree H is requested.");
@@ -2794,16 +1012,13 @@ void AliRICH::SetTreeAddress()
 //kir       }
   }//if(treeH)
  
-   //this is after TreeH because we need to guarantee that fHits array is created
-  AliDetector::SetTreeAddress();
+  AliDetector::SetTreeAddress();//this is after TreeH because we need to guarantee that fHits array is created
+
     
-  if(treeS){
+  if(fLoader->TreeS()){
     if(GetDebug())Info("SetTreeAddress","tree S is requested.");
-    branch = treeS->GetBranch("RICH");
-    if(branch){
-      if(!fSDigits) fSDigits=new TClonesArray("AliRICHdigit",100000);
-      branch->SetAddress(&fSDigits);
-    }
+    if(!fSDigits) fSDigits=new TClonesArray("AliRICHdigit",100000);
+    fLoader->TreeS()->GetBranch("RICH")->SetAddress(&fSDigits);
   }
     
     
@@ -2885,14 +1100,14 @@ void AliRICH::SetTreeAddress()
   }//if(treeR)
   if(GetDebug())Info("SetTreeAddress","Stop.");
 }//void AliRICH::SetTreeAddress()
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::Print(Option_t *option)const
 {
   TObject::Print(option);
   fpParam->Dump();
-  Chambers()->Print(option);  
+  fChambers->Print(option);  
 }//void AliRICH::Print(Option_t *option)const
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::CreateGeometry()
 {//Creates detailed geometry simulation (currently GEANT volumes tree)         
   if(GetDebug())Info("CreateGeometry","Start.");
@@ -2902,8 +1117,8 @@ void AliRICH::CreateGeometry()
 //Opaque quartz thickness
   Float_t oqua_thickness = .5;
 //CsI dimensions
-  Float_t csi_width =fpParam->Nx()*fpParam->PadX()+fpParam->DeadZone();
-  Float_t csi_length=fpParam->Ny()*fpParam->PadY()+2*fpParam->DeadZone();
+  Float_t csi_width =fpParam->Nx()*fpParam->PadSizeX()+fpParam->DeadZone();
+  Float_t csi_length=fpParam->Ny()*fpParam->PadSizeY()+2*fpParam->DeadZone();
   
   Int_t *idtmed = fIdtmed->GetArray()-999;
     
@@ -2913,7 +1128,7 @@ void AliRICH::CreateGeometry()
   Float_t par[3];
     
 //External aluminium box 
-  par[0]=68.8;par[1]=13;par[2]=70.86;//Original Settings
+  par[0]=68.8*cm;par[1]=13*cm;par[2]=70.86*cm;//Original Settings
   gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
 //Air 
   par[0]=66.3;   par[1] = 13; par[2] = 68.35; //Original Settings
@@ -3110,13 +1325,11 @@ void AliRICH::CreateGeometry()
   gMC->Gspos("PCB ", 2, "SRIC ", 0.,  1.276 + fpParam->GapThickness()/2 + .5 + 1.05, -csi_width/4 - .5025 - 2.5, 0, "ONLY");
 
 //place chambers into mother volume ALIC
-  CreateChambers();
-
-  for(int i=0;i<kNCH;i++){
+  for(int i=1;i<=kNCH;i++){
     AliMatrix(idrotm[1000+i],C(i)->ThetaXd(),C(i)->PhiXd(),
                              C(i)->ThetaYd(),C(i)->PhiYd(),
                              C(i)->ThetaZd(),C(i)->PhiZd());
-    gMC->Gspos("RICH",i+1,"ALIC",C(i)->X(),C(i)->Y(),C(i)->Z(),idrotm[1000+i], "ONLY");
+    gMC->Gspos("RICH",i,"ALIC",C(i)->X(),C(i)->Y(),C(i)->Z(),idrotm[1000+i], "ONLY");
   }
 
   if(GetDebug())Info("CreateGeometry","Stop.");  
@@ -3135,3 +1348,114 @@ void AliRICH::CreateChambers()
 
   if(GetDebug())Info("CreateChambers","Stop.");
 }//void AliRICH::CreateChambers()
+//__________________________________________________________________________________________________
+void AliRICH::GenerateFeedbacks(Float_t eloss)
+{// Generate FeedBack photons
+  Int_t j;
+  Float_t cthf, phif, enfp = 0, sthf;
+  Float_t e1[3], e2[3], e3[3];
+  Float_t vmod, uswop;
+  Float_t dir[3], phi;
+  Float_t pol[3], mom[4];
+//Determine number of feedback photons
+  TLorentzVector x4;
+  gMC->TrackPosition(x4);//This sould return the current track position
+  Float_t charge=Param()->TotalCharge(gMC->TrackPid(),eloss,5*cm);//??? Totsl Charge
+  Int_t iNphotons=gMC->GetRandom()->Poisson(Param()->AlphaFeedback()*charge);    
+  Info("GenerateFeedbacks","N photons=%i",iNphotons);
+//Generate photons
+  for(Int_t i=0;i<iNphotons;i++){
+    Double_t ranf[2];
+    gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
+    cthf=ranf[0]*2-1.0;
+    if(cthf<0) continue;
+    sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
+    phif = ranf[1] * 2 * TMath::Pi();
+    
+    if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
+      enfp = 7.5e-9;
+    else if(randomNumber<=0.7)
+      enfp = 6.4e-9;
+    else
+      enfp = 7.9e-9;
+    
+
+    dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
+    gMC->Gdtom(dir, mom, 2);
+    mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
+    mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
+    
+    // Polarisation
+    e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
+    e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
+    e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
+    
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+    if (!vmod) for(j=0;j<3;j++) {
+      uswop=e1[j];
+      e1[j]=e3[j];
+      e3[j]=uswop;
+    }
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+    if (!vmod) for(j=0;j<3;j++) {
+      uswop=e2[j];
+      e2[j]=e3[j];
+      e3[j]=uswop;
+    }
+    
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+    vmod=TMath::Sqrt(1/vmod);
+    for(j=0;j<3;j++) e1[j]*=vmod;
+    
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+    vmod=TMath::Sqrt(1/vmod);
+    for(j=0;j<3;j++) e2[j]*=vmod;
+    
+    phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
+    for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
+    gMC->Gdtom(pol, pol, 2);
+    Int_t outputNtracksStored;    
+    gAlice->PushTrack(1,                             //do not transport
+                     gAlice->GetCurrentTrackNumber(),//parent track 
+                     kFeedback,                      //PID
+                    mom[0],mom[1],mom[2],mom[3],    //track momentum  
+                     x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin 
+                     pol[0],pol[1],pol[2],           //polarization
+                    kPFeedBackPhoton,outputNtracksStored,1.0);
+    
+  }
+}//Int_t AliRICH::FeedBackPhotons()
+//__________________________________________________________________________________________________
+static Int_t sMaxIterPad=0;    // Static variables for the pad-hit iterator routines
+static Int_t sCurIterPad=0;
+
+//__________________________________________________________________________________________________
+AliRICHSDigit* AliRICH::FirstPad(AliRICHhit*  hit,TClonesArray *clusters ) 
+{// Initialise the pad iterator Return the address of the first sdigit for hit
+    TClonesArray *theClusters = clusters;
+    Int_t nclust = theClusters->GetEntriesFast();
+    if (nclust && hit->PHlast() > 0) {
+       sMaxIterPad=Int_t(hit->PHlast());
+       sCurIterPad=Int_t(hit->PHfirst());
+       return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
+    } else {
+       return 0;
+    }
+    
+}
+//__________________________________________________________________________________________________
+AliRICHSDigit* AliRICH::NextPad(TClonesArray *clusters) 
+{// Iterates over pads
+  
+    sCurIterPad++;
+    if (sCurIterPad <= sMaxIterPad) {
+       return (AliRICHSDigit*) clusters->UncheckedAt(sCurIterPad-1);
+    } else {
+       return 0;
+    }
+}
+//__________________________________________________________________________________________________
index 991e238..966f89c 100644 (file)
@@ -4,51 +4,39 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id$ */
-
-
 #include <TObjArray.h>
 #include <TClonesArray.h>
+#include <TLorentzVector.h>
 #include <AliDetector.h>
 #include <AliHit.h>
 #include <AliDigit.h>
 #include "AliRICHConst.h"
 #include "AliRICHChamber.h"
 
-static const int kNCH=7;
-
-class TFile;
-
 class AliRICHRawCluster;
 class AliRICHRecHit1D;
 class AliRICHRecHit3D;
-class AliRICHClusterFinder;
-class AliRICHDetect;
-class AliRICHChamber;
-class AliRICHCerenkov;
-class AliSegmentation;
-class AliRICHResponse;
-class AliRICHGeometry;
-class AliRICHMerger;
-
-
 
+//__________________AliRICHhit______________________________________________________________________
+//__________________________________________________________________________________________________
+//__________________________________________________________________________________________________
 class AliRICHhit : public AliHit
 {
 public:
   inline   AliRICHhit();
   inline   AliRICHhit(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *hits);
+  inline   AliRICHhit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss);
   virtual ~AliRICHhit()         {;}
     
   Int_t   Chamber()             {return fChamber;}
-  Float_t Particle()            {return fParticle;}    
+  Int_t   Particle()            {return fParticle;}    
   Float_t Theta()               {return fTheta;}
   Float_t Phi()                 {return fPhi;}
   Float_t Tlength()             {return fTlength;}
   Float_t Eloss()               {return fEloss;}
   Float_t Loss()                {return fLoss;}
-  Float_t   PHfirst()           {return fPHfirst;}
-  Float_t   PHlast()            {return fPHlast;}
+  Float_t PHfirst()             {return fPHfirst;}
+  Float_t PHlast()              {return fPHlast;}
   Float_t MomX()                {return fMomX;}
   Float_t MomY()                {return fMomY;}
   Float_t MomZ()                {return fMomZ;}
@@ -57,62 +45,60 @@ public:
   Float_t MomFreoY()            {return fMomY;}
   Float_t MomFreoZ()            {return fMomZ;}
 protected:
-  Int_t     fChamber;       // Chamber number
-  Float_t   fParticle;      // Geant3 particle type
-  Float_t   fTheta ;        // Incident theta angle in degrees      
-  Float_t   fPhi   ;        // Incident phi angle in degrees
-  Float_t   fTlength;       // Track length inside the chamber
-  Float_t   fEloss;         // ionisation energy loss in gas   
-  Float_t   fPHfirst;       // first padhit
-  Float_t   fPHlast;        // last padhit
-  Float_t   fLoss;          // did it hit the freon?
-  Float_t   fMomX;          // x Momentum at photochatode entry point
-  Float_t   fMomY;          // y Momentum at photochatode entry point
-  Float_t   fMomZ;          // z Momentum at photochatode entry point
-  Float_t   fNPads;         // Pads hit
-  Float_t   fCerenkovAngle; // Dummy cerenkov angle
-  Float_t   fMomFreoX;      // x Momentum at freon entry point
-  Float_t   fMomFreoY;      // y Momentum at freon entry point
-  Float_t   fMomFreoZ;      // z Momentum at freon entry point            
-        
-  ClassDef(AliRICHhit,1)  //RICH hit class
+  Int_t     fChamber;                      //chamber number
+  Int_t     fParticle;                     //particle code
+  Float_t   fTheta,fPhi ;                  //incident theta phi angles in degrees      
+  Float_t   fTlength;                      //track length inside the chamber
+  Float_t   fEloss;                        //ionisation energy loss in gas   
+  Float_t   fPHfirst;                      //first padhit
+  Float_t   fPHlast;                       //last padhit
+  Float_t   fLoss;                         // did it hit the freon?
+  Float_t   fMomX,fMomY,fMomZ;             //momentum at photochatode entry point
+  Float_t   fNPads;                        // Pads hit
+  Float_t   fCerenkovAngle;                // Dummy cerenkov angle
+  Float_t   fMomFreoX,fMomFreoY,fMomFreoZ; //momentum at freon entry point
+  ClassDef(AliRICHhit,1)                   //RICH hit class
 };//class AliRICHhit
-//______________________________________________________________________________  
+
+  //__________________________________________________________________________________________________
 AliRICHhit::AliRICHhit()
            :AliHit() 
 {//default ctor  
-  fChamber=-1;
-  fParticle=fTheta=fPhi=fTlength=fEloss=fPHfirst=fPHlast=fLoss=-1;
-  fMomX=fMomY=fMomZ=fNPads=fCerenkovAngle=fMomFreoX=fMomFreoY=fMomFreoZ=-1;
+  fChamber=fParticle=kBad;
+  fTheta=fPhi=fTlength=fEloss=fPHfirst=fPHlast=fLoss=kBad;
+  fMomX=fMomY=fMomZ=fNPads=fCerenkovAngle=fMomFreoX=fMomFreoY=fMomFreoZ=kBad;
 }//AliRICHhit::default ctor
-//______________________________________________________________________________
-AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
+//__________________________________________________________________________________________________
+AliRICHhit::AliRICHhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hit):
             AliHit(shunt, track)
 {//ctor
   fChamber=vol[0];
-  fParticle=hits[0];
-  fX=hits[1];
-  fY=hits[2];
-  fZ=hits[3];
-  fTheta=hits[4];
-  fPhi=hits[5];
-  fTlength=hits[6];
-  fEloss=hits[7];
-  fPHfirst=(Int_t) hits[8];
-  fPHlast=(Int_t) hits[9];
-  fLoss=hits[13];
-  fMomX=hits[14];
-  fMomY=hits[15];
-  fMomZ=hits[16];
-  fNPads=hits[17];
-  fCerenkovAngle=hits[18];
-  fMomFreoX=hits[19];
-  fMomFreoY=hits[20];
-  fMomFreoZ=hits[21];
+  fParticle=(Int_t)hit[0];
+  fX=hit[1];fY=hit[2];fZ=hit[3];
+  fTheta=hit[4];fPhi=hit[5];
+  fTlength=hit[6];
+  fEloss=hit[7];
+  fPHfirst=(Int_t)hit[8];
+  fPHlast=(Int_t)hit[9];
+  fLoss=hit[13];
+  fMomX=hit[14];fMomY=hit[15];fMomZ=hit[16];
+  fNPads=hit[17];
+  fCerenkovAngle=hit[18];
+  fMomFreoX=hit[19];fMomFreoY=hit[20];fMomFreoZ=hit[21];
+}//AliRICHhit::ctor
+//__________________________________________________________________________________________________
+AliRICHhit::AliRICHhit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss):
+            AliHit(0, track)
+{//ctor
+  fChamber=iChamber;
+  fParticle=iPID;
+  fX=x4.X();fY=x4.Y();fZ=x4.Z();
+  fEloss=eloss;
 }//AliRICHhit::ctor
-//______________________________________________________________________________
-//______________________________________________________________________________
-//______________________________________________________________________________
+
+//__________________AliRICHCerenkov_________________________________________________________________
+//__________________________________________________________________________________________________
+//__________________________________________________________________________________________________
 class AliRICHCerenkov: public AliHit 
 {
 public:
@@ -120,86 +106,82 @@ public:
   inline   AliRICHCerenkov(Int_t fIshunt, Int_t track, Int_t *vol, Float_t *Cerenkovs);
   virtual ~AliRICHCerenkov() {;}
 public:
-  Int_t     fChamber;         // Chamber number
-  Float_t   fTheta ;          // Incident theta angle in degrees      
-  Float_t   fPhi   ;          // Incident phi angle in degrees
-  Float_t   fTlength;         // Track length inside the chamber
-  Float_t   fEloss;           // ionisation energy loss in gas
-  Int_t     fPHfirst;         // first padhit
-  Int_t     fPHlast;          // last padhit
-  Int_t     fCMother;         // index of mother particle
-  Float_t   fLoss;            // nature of particle loss
-  Float_t   fIndex;           // Index of photon
-  Float_t   fProduction;      // Point of production
-  Float_t   fMomX;            // Local Momentum
-  Float_t   fMomY;            // Local Momentum
-  Float_t   fMomZ;            // Local Momentum
+  Int_t     fChamber;          //chamber number
+  Float_t   fTheta,fPhi;       //incident theta phi angles in degrees      
+  Float_t   fTlength;          //track length inside the chamber
+  Float_t   fEloss;            //ionisation energy loss in gas
+  Int_t     fPHfirst;          //first padhit
+  Int_t     fPHlast;           //last padhit
+  Int_t     fCMother;          //index of mother particle
+  Float_t   fLoss;             //nature of particle loss
+  Float_t   fIndex;            //index of photon
+  Float_t   fProduction;       //point of production
+  Float_t   fMomX,fMomY,fMomZ; //local Momentum
   Float_t   fNPads;           // Pads hit
   Float_t   fCerenkovAngle;   // Cerenkov Angle
     
   ClassDef(AliRICHCerenkov,1)  //RICH cerenkov class
 };//class AliRICHCerenkov
-//______________________________________________________________________________
+
+//__________________________________________________________________________________________________
 AliRICHCerenkov::AliRICHCerenkov()
 {//ctor
-    fChamber=-1;
-    fX=fY=fZ=fTheta=fPhi=fTlength=fEloss=-1;
-    fPHfirst=fPHlast=fCMother=-1;
-    fLoss=fIndex=fProduction=fMomX=fMomY=fMomZ=fNPads=fCerenkovAngle=-1;
+    fChamber=kBad;
+    fX=fY=fZ=fTheta=fPhi=fTlength=fEloss=kBad;
+    fPHfirst=fPHlast=fCMother=kBad;
+    fLoss=fIndex=fProduction=fMomX=fMomY=fMomZ=fNPads=fCerenkovAngle=kBad;
 }//AliRICHCerenkov::ctor
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRICHCerenkov::AliRICHCerenkov(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
                 :AliHit(shunt, track)
 {//ctor
     fChamber=vol[0];
-    fX=hits[1];
-    fY=hits[2];
-    fZ=hits[3];
-    fTheta=hits[4];
-    fPhi=hits[5];
+    fX=hits[1];fY=hits[2];fZ=hits[3];
+    fTheta=hits[4];fPhi=hits[5];
     fTlength=hits[6];
     fEloss=hits[7];
-    fPHfirst=(Int_t) hits[8];
-    fPHlast=(Int_t) hits[9];
+    fPHfirst=(Int_t)hits[8];fPHlast=(Int_t)hits[9];
     fCMother=Int_t(hits[10]);
     fIndex = hits[11];
     fProduction = hits[12];  
     fLoss=hits[13];
-    fMomX=hits[14];
-    fMomY=hits[15];
-    fMomZ=hits[16];
+    fMomX=hits[14];fMomY=hits[15];fMomZ=hits[16];
     fNPads=hits[17];
     fCerenkovAngle=hits[18];
 }//AliRICHCerenkov::ctor
-//______________________________________________________________________________
-//______________________________________________________________________________
-//______________________________________________________________________________
+
+//__________________AliRICHdigit____________________________________________________________________
+//__________________________________________________________________________________________________
+//__________________________________________________________________________________________________
 class AliRICHdigit :public AliDigit
 {
 public:
-           AliRICHdigit() {fPadX=fPadY=fChamber=fAdc=fTracks[0]=fTracks[1]=fTracks[2]=-1;}
+           AliRICHdigit() {fPadX=fPadY=fChamber=fAdc=fTracks[0]=fTracks[1]=fTracks[2]=kBad;}
   inline   AliRICHdigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT1,Int_t iT2,Int_t iT3);
   virtual ~AliRICHdigit() {;}  
   Int_t C()   const{return fChamber;}
   Int_t X()   const{return fPadX;}
   Int_t Y()   const{return fPadY;}
+  Int_t Pad() const{return fPad;}
   Int_t Adc() const{return fAdc;}
 protected:
   Int_t fChamber;  //module number 
-  Int_t fPadX;    //pad number along X
-  Int_t fPadY;    //pad number along Y
-  Int_t fAdc;     //ADC value
+  Int_t fPadX;     //pad number along X
+  Int_t fPadY;     //pad number along Y
+  Int_t fPad;      //pad number 1000*X+Y
+  Int_t fAdc;      //ADC value
   ClassDef(AliRICHdigit,1) //RICH digit class       
 };//class AliRICHdigit
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRICHdigit::AliRICHdigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1,Int_t iT2)
 {
-  fChamber=iC;fPadX=iX;fPadY=iY;fAdc=iAdc;
+  fChamber=iC;fPadX=iX;fPadY=iY;fPad=1000*fPadX+fPadY;fAdc=iAdc;
   fTracks[0]=iT0;fTracks[1]=iT1;fTracks[2]=iT2;
 }//AliRICHdigit::ctor  
-//______________________________________________________________________________
-//______________________________________________________________________________
-//______________________________________________________________________________
+
+//__________________AliRICH_________________________________________________________________________
+//__________________________________________________________________________________________________
+//__________________________________________________________________________________________________
 class AliRICHParam;
 class AliRICHSDigit;
 
@@ -213,61 +195,52 @@ public:
           
   AliRICH&  operator=(const AliRICH&)                 {return *this;}
   virtual Int_t  IsVersion()const =0;            
-  AliRICHParam *Param()                               {return fpParam;}
-  inline  void  AddHit(Int_t track, Int_t *vol, Float_t *hits);//virtual
-  inline  void  AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs);
-  inline  void  AddSDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1=-1,Int_t iT2=-1);
-  inline  void  ResetHits();    //virtual
-  inline  void  ResetSDigits(); //virtual
-          void  Hits2SDigits(); //virtual 
+  
+  inline  void    AddHit(Int_t track, Int_t *vol, Float_t *hits);//virtual
+  inline  void    AddHit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss);
+  inline  void    AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs);
+  inline  void    AddSDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1=kBad,Int_t iT2=kBad);
+  inline  void    AddDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1=kBad,Int_t iT2=kBad);
+  
+  inline  void    ResetHits();    //virtual
+  inline  void    ResetSDigits(); 
+          void    ResetDigits();  //virtual
+  
+  virtual void    Hits2SDigits();   
+  virtual void    SDigits2Digits();
+  virtual void    Digits2Reco();   
           
-  TClonesArray  *SDigits()                       const{return fSDigits;}
-  TClonesArray  *Cerenkovs()                     const{return fCerenkovs;}
+  AliRICHChamber* C(Int_t i)                       const{return (AliRICHChamber*)fChambers->At(i-1);}
+  AliRICHParam*   Param()                          const{return fpParam;}
+  TClonesArray*   SDigits()                        const{return fSDigits;}
+  TClonesArray*   Cerenkovs()                      const{return fCerenkovs;}
+  
           void    CreateChambers();         
   virtual void    CreateMaterials(); //GEANT materials definition
           Float_t AbsoCH4(Float_t x);
           Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola);
   virtual void    BuildGeometry();   //TNode ROOT variant for event display
   virtual void    CreateGeometry();  //GEANT volumes tree for simulation  
-  virtual void    StepManager()=0;
   
-  AliRICHChamber* C(Int_t i)        const{return (AliRICHChamber*)fChambers->At(i);}//return pointer to requested chamber
+  virtual void    StepManager()=0;
+          void    GenerateFeedbacks(Float_t eloss);
+            
   AliRICHChamber& Chamber(Int_t id)      {return *((AliRICHChamber *) (*fChambers)[id]);}
-  TObjArray*      Chambers()        const{return fChambers;}
   
-        void  AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits);
+          void  AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits);
           void  AddRawCluster(Int_t id, const AliRICHRawCluster& cluster);
           void  AddRecHit1D(Int_t id, Float_t* rechit, Float_t* photons, Int_t* padsx, Int_t* padsy);
           void  AddRecHit3D(Int_t id, Float_t* rechit, Float_t omega, Float_t theta, Float_t phi);
-          void  ResetDigits();  //virtual
           void  ResetRawClusters();
           void  ResetRecHits1D();
           void  ResetRecHits3D();
   virtual void  FindClusters(Int_t nev);
-          Int_t Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss,Int_t id, ResponseType res);//kir ????? to be  removed
-  virtual void  SDigits2Digits();
-  virtual void  Digits2Reco();   
-          Int_t    DistancetoPrimitive(Int_t /*px*/, Int_t /*py*/)      {return 9999;}
+          Int_t DistancetoPrimitive(Int_t /*px*/, Int_t /*py*/)      {return 9999;}
    
   virtual void   MakeBranch(Option_t *opt=" ");
   virtual void   MakeBranchInTreeD(TTree *treeD, const char *file=0);
   virtual void   SetTreeAddress();
-   
-     
-   
-  AliRICHSDigit* FirstPad(AliRICHhit *hit, TClonesArray *clusters);
-  AliRICHSDigit* NextPad(TClonesArray *clusters);
-   
-
-  void     SetGeometryModel(Int_t iC,AliRICHGeometry *pRICHGeo)                    {C(iC)->SetGeometryModel(pRICHGeo);}
-  void     SetSegmentationModel(Int_t iC, AliSegmentation *pAliSeg)                {C(iC)->SetSegmentationModel(pAliSeg);}
-  void     SetResponseModel(Int_t iC, AliRICHResponse *pRICHRes)                   {C(iC)->SetResponseModel(pRICHRes);}
-  void     SetReconstructionModel(Int_t iC, AliRICHClusterFinder *pRICHReco)       {C(iC)->SetReconstructionModel(pRICHReco);}
-  AliRICHGeometry* GetGeometryModel(Int_t iC=0)                               const{return C(iC)->GetGeometryModel();}    
-  AliSegmentation* GetSegmentationModel(Int_t iC=0)                           const{return C(iC)->GetSegmentationModel();}
-  AliRICHResponse* GetResponseModel(Int_t iC)                                 const{return C(iC)->GetResponseModel();}
-
-//kir  virtual void   SetMerger(AliRICHMerger* thisMerger) {fMerger=thisMerger;}  
+              
   
   TObjArray     *Dchambers()                     {return fDchambers;}
   TObjArray     *RecHits3D()                const{return fRecHits3D;}
@@ -279,24 +252,24 @@ public:
   TClonesArray  *RecHitsAddress1D(Int_t id) const{return ((TClonesArray *) (*fRecHits1D)[id]);}
   TClonesArray  *RecHitsAddress3D(Int_t id) const{return ((TClonesArray *) (*fRecHits3D)[id]);}
   TClonesArray  *RawClustAddress(Int_t id)  const{return ((TClonesArray *) (*fRawClusters)[id]);}    
+  AliRICHSDigit* FirstPad(AliRICHhit *hit, TClonesArray *clusters);
+  AliRICHSDigit* NextPad(TClonesArray *clusters);
 
-  void DiagnosticsFE(Int_t evNumber1=0,Int_t evNumber2=0);    // Full events
-  void DiagnosticsSE(Int_t diaglevel,Int_t evNumber1=0,Int_t evNumber2=0);    // Single events
  
-  virtual void Print(Option_t *option)const; // Prints debug information
+  virtual void Print(Option_t *option)const;
     
 protected:
-  AliRICHParam         *fpParam;              //main RICH parametrization     
-  TObjArray            *fChambers;           //! List of RICH chambers
+  AliRICHParam         *fpParam;             //main RICH parametrization     
+  TObjArray            *fChambers;           //list of RICH chambers
   Int_t                 fNsdigits;           //Current number of sdigits
   Int_t                 fNcerenkovs;         //Current number of cerenkovs
   TClonesArray         *fSDigits;            //! List of sdigits
   TObjArray            *fDchambers;          //! Array of lists of digits
   TClonesArray         *fCerenkovs;          //! List of cerenkovs
   Int_t                 fNdch[kNCH];         //Array of current numbers of digits
-  TObjArray            *fRawClusters;        // !List of raw clusters
-  TObjArray            *fRecHits1D;          // !List of rec. hits
-  TObjArray            *fRecHits3D;          // !List of rec. hits
+  TObjArray            *fRawClusters;        //!List of raw clusters
+  TObjArray            *fRecHits1D;          //!List of rec. hits
+  TObjArray            *fRecHits3D;          //!List of rec. hits
   Int_t                 fNrawch[kNCH];       //Array of current numbers of raw clusters
   Int_t                 fNrechits1D[kNCH];   //Array of current numbers of rec hits 1D
   Int_t                 fNrechits3D[kNCH];   //Array of current numbers of rec hits 3D 
@@ -304,56 +277,52 @@ protected:
   Int_t fCkovNumber;                         // Number of Cerenkov photons
   Int_t fFreonProd;                          // Cerenkovs produced in freon
   Int_t fFeedbacks;                          // Number of feedback photons
-//kir  Int_t fCkovQuarz;                          // Cerenkovs crossing quartz
-//kir  Int_t fCkovGap;                            // Cerenkovs crossing gap
-//kir  Int_t fCkovCsi;                            // Cerenkovs crossing csi
-//kir  Int_t fLostRfreo;                          // Cerenkovs reflected in freon
-//kir  Int_t fLostRquar;                          // Cerenkovs reflected in quartz
-//kir  Int_t fLostAfreo;                          // Cerenkovs absorbed in freon 
-//kir  Int_t fLostAquarz;                         // Cerenkovs absorbed in quartz
-//kir  Int_t fLostAmeta;                          // Cerenkovs absorbed in methane
-//kir  Int_t fLostCsi;                            // Cerenkovs below csi quantum efficiency 
-//kir  Int_t fLostWires;                          // Cerenkovs lost in wires
-//kir  Float_t fMipx;                             // x coord. of MIP
-//kir  Float_t fMipy;                             // y coord. of MIP
-//kir  Int_t fLostFresnel;                        // Cerenkovs lost by Fresnel reflection
-
-//kir  Text_t *fFileName;                         //! File with background hits
-//kir  AliRICHMerger *fMerger;                    //! pointer to merger
     
   ClassDef(AliRICH,2)                        //Main RICH class 
 };//class AliRICH
-
-//______________________________________________________________________________
+  
+//__________________________________________________________________________________________________
 void AliRICH::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 {//Adds the current hit to the RICH hits list
   TClonesArray &tmp=*fHits;
   new(tmp[fNhits++])AliRICHhit(fIshunt,track,vol,hits);
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
+void AliRICH::AddHit(Int_t track,Int_t iPID,Int_t iChamber,TLorentzVector x4,Float_t eloss)
+{//Adds the current hit to the RICH hits list
+  TClonesArray &tmp=*fHits;
+  new(tmp[fNhits++])AliRICHhit(track,iPID,iChamber,x4,eloss);
+}
+//__________________________________________________________________________________________________
 void AliRICH::AddCerenkov(Int_t track, Int_t *vol, Float_t *cerenkovs)
 {//Adds the current RICH cerenkov hit to the Cerenkovs list   
   TClonesArray &tmp=*fCerenkovs;
   new(tmp[fNcerenkovs++]) AliRICHCerenkov(fIshunt,track,vol,cerenkovs);
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::ResetHits()
 {//Resets hits and cerenkovs
   AliDetector::ResetHits();
-  fNcerenkovs = 0;
+  fNcerenkovs=0;
   if(fCerenkovs)fCerenkovs->Clear();
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICH::AddSDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1,Int_t iT2)
 {//Adds the current Sdigit to the RICH list of Sdigits   
-  TClonesArray &tmp = *fSDigits;
+  TClonesArray &tmp=*fSDigits;
   new(tmp[fNsdigits++])AliRICHdigit(iC,iX,iY,iAdc,iT0,iT1,iT2);
 } 
-//______________________________________________________________________________    
+//__________________________________________________________________________________________________
 void AliRICH::ResetSDigits()
 {//Resets sdigits
   fNsdigits=0;
   if(fSDigits)fSDigits->Clear();
 }
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
+void AliRICH::AddDigit(Int_t iC,Int_t iX,Int_t iY,Int_t iAdc,Int_t iT0,Int_t iT1,Int_t iT2)
+{//Adds the current Digit to the RICH list of digits
+
+   TClonesArray &tmp=*fDigits;
+   new(tmp[fNdigits++]) AliRICHdigit(iC,iX,iY,iAdc,iT0,iT1,iT2);
+}
 #endif
index 67a6cab..be2d844 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id$ */
-
-
-/////////////////////////////////////////////////////////////
-//  Manager and hits classes for set: RICH default version //
-/////////////////////////////////////////////////////////////
-
-#include "Riostream.h"
-
-#include <TNode.h> 
-#include <TRandom.h> 
-#include <TTUBE.h>
-#include <TVirtualMC.h>
-
-#include "AliConst.h" 
-#include "AliPDG.h" 
-#include "AliRICHGeometry.h"
-#include "AliRICHResponse.h"
-#include "AliRICHResponseV0.h"
-#include "AliRICHSegmentationV0.h"
 #include "AliRICHv0.h"
-#include "AliRun.h"
+#include <TVirtualMC.h>
+#include <TPDGCode.h>
+#include "AliRICHConst.h" 
+#include <AliRun.h>
+#include <TLorentzVector.h>
 
 ClassImp(AliRICHv0)
-    
-//___________________________________________
+//__________________________________________________________________________________________________
 AliRICHv0::AliRICHv0(const char *name, const char *title)
           :AliRICH(name,title)
 {
-//
-// Version 0
-// Default Segmentation, no hits
-    AliRICHSegmentationV0* segmentation = new AliRICHSegmentationV0;
-//
-//  Segmentation parameters
-    segmentation->SetPadSize(0.84,0.80);
-    segmentation->SetDAnod(0.84/2);
-//
-//  Geometry parameters
-    AliRICHGeometry* geometry = new AliRICHGeometry;
-    geometry->SetGapThickness(8);
-    geometry->SetProximityGapThickness(.4);
-    geometry->SetQuartzLength(131);
-    geometry->SetQuartzWidth(126.2);
-    geometry->SetQuartzThickness(.5);
-    geometry->SetOuterFreonLength(131);
-    geometry->SetOuterFreonWidth(40.3);
-    geometry->SetInnerFreonLength(131);
-    geometry->SetInnerFreonWidth(40.3);
-    geometry->SetFreonThickness(1.5);
-//
-//  Response parameters
-    AliRICHResponseV0*  response   = new AliRICHResponseV0;
-    response->SetSigmaIntegration(5.);
-    response->SetChargeSlope(27.);
-    response->SetChargeSpread(0.18, 0.18);
-    response->SetMaxAdc(4096);
-    response->SetAlphaFeedback(0.036);
-    response->SetEIonisation(26.e-9);
-    response->SetSqrtKx3(0.77459667);
-    response->SetKx2(0.962);
-    response->SetKx4(0.379);
-    response->SetSqrtKy3(0.77459667);
-    response->SetKy2(0.962);
-    response->SetKy4(0.379);
-    response->SetPitch(0.25);
-    response->SetWireSag(0);                     // 1->On, 0->Off
-    response->SetVoltage(2150);                  // Should only be 2000, 2050, 2100 or 2150
-//
-//    AliRICH *RICH = (AliRICH *) gAlice->GetDetector("RICH"); 
-    
-    fCkovNumber=0;
-    fFreonProd=0;
-    Int_t i=0;
-    
-    fChambers = new TObjArray(kNCH);
-    for (i=0; i<kNCH; i++) {
-      
-      //PH      (*fChambers)[i] = new AliRICHChamber();  
-      fChambers->AddAt(new AliRICHChamber(), i);  
-      
-    }
-  
-    for (i=0; i<kNCH; i++) {
-      SetGeometryModel(i,geometry);
-      SetSegmentationModel(i, segmentation);
-      SetResponseModel(i, response);
-    }
+  if(GetDebug())Info("named ctor","Start.");
+  if(GetDebug())Info("named ctor","Stop.");
 }//name ctor
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRICHv0::StepManager()
 {//
+//  if(!gMC->IsNewTrack()) return;
+  char *sParticle;
+  switch(gMC->TrackPid()){
+    case kProton:
+      sParticle="p";break;
+    case kNeutron:
+      sParticle="n";break;
+    case kGamma:
+      sParticle="gamma";break;
+    case 50000050:
+      sParticle="photon";break;
+    default:
+      sParticle="not known";break;
+  }
+
+  Info("StepManager","Event=%i hunt=%i TID=%i PID=%s Mass=%f Charge=%i",
+                      gMC->CurrentEvent(),
+                                fIshunt,
+                                        gAlice->GetCurrentTrackNumber(),
+                                               sParticle,
+                                                      gMC->TrackMass(),
+                                                              gMC->TrackCharge());
+  Info("StepManager","Flags:Alive(%i) Disap(%i) Enter(%i) Exit(%i) Inside(%i) Out(%i) Stop(%i) New(%i)",
+                            gMC->IsTrackAlive(),
+                                      gMC->IsTrackDisappeared(),
+                                                gMC->IsTrackEntering(),
+                                                          gMC->IsTrackExiting(),
+                                                                    gMC->IsTrackInside(),
+                                                                          gMC->IsTrackOut(),
+                                                                                  gMC->IsTrackStop(),
+                                                                                       gMC->IsNewTrack());
+  Info("StepManager","Volume=%s of volume=%s",
+                      gMC->CurrentVolName(),gMC->CurrentVolOffName(1));
+
+//  Info("StepManager","TrackPID %i Particle %i",
+//                     gMC->TrackPid(),gAlice->Particles()[gAlice->CurrentTrack()]
+  TLorentzVector x4;
+  gMC->TrackPosition(x4);
+  Info("StepManager","x=%f y=%f z=%f r=%f theta=%f phi=%f\n",
+                      x4.X(),x4.Y(),x4.Z(),x4.Rho(),x4.Theta()*r2d,x4.Phi()*r2d);  
 }//AliRICHv0::StepManager()
-//______________________________________________________________________________
+//__________________________________________________________________________________________________
index 8570e4b..43e131c 100644 (file)
 #include <TVector3.h>
 #include <TVirtualMC.h>
 #include <TPDGCode.h> //for kNuetron
+#include <TCanvas.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TStyle.h>
 
 #include "AliConst.h"
 #include "AliMagF.h"
 #include "AliRICHSegmentationV1.h"
 #include "AliRICHv3.h"
 #include "AliRun.h"
+#include "AliRICHRecHit3D.h"
+#include "AliRICHRawCluster.h"
+#include "AliRICHDigit.h"
+#include "AliRICHRecHit1D.h"
+
 
 ClassImp(AliRICHv3)
 
@@ -55,13 +65,11 @@ AliRICHv3::AliRICHv3(const char *sName, const char *sTitle)
    AliRICHSegmentationV1 *pRICHSegmentation=new AliRICHSegmentationV1;     // ??? to be moved to AlRICHChamber::named ctor
    AliRICHResponseV0     *pRICHResponse    =new AliRICHResponseV0;         // ??? to be moved to AlRICHChamber::named ctor
      
-   fChambers = new TObjArray(kNCH);
-   for (Int_t i=0; i<kNCH; i++){    
-      fChambers->AddAt(new AliRICHChamber,i); // ??? to be changed to named ctor of AliRICHChamber
+   for (Int_t i=1; i<=kNCH; i++){    
       SetGeometryModel(i,pRICHGeometry);
       SetSegmentationModel(i,pRICHSegmentation);
       SetResponseModel(i,pRICHResponse);
-      ((AliRICHChamber*)fChambers->At(i))->Init(i); // ??? to be removed     
+      C(i)->Init(i); // ??? to be removed     
    }
   if(GetDebug())Info("named ctor","Stop.");
 }//AliRICHv3::ctor(const char *pcName, const char *pcTitle)
@@ -72,7 +80,7 @@ AliRICHv3::~AliRICHv3()
    if(GetDebug()) cout<<ClassName()<<"::dtor()>\n";
       
    if(fChambers) {
-     AliRICHChamber *ch =C(0); 
+     AliRICHChamber *ch =C(1); 
      if(ch) {
        delete ch->GetGeometryModel();
        delete ch->GetResponseModel();
@@ -1082,7 +1090,7 @@ void AliRICHv3::StepManager()
                    cherenkovLoss  += destep;
                    ckovData[7]=cherenkovLoss;
                    
-                   ckovData[17] = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kCerenkov);//for photons in CsI 
+                   ckovData[17] = Hits2SDigits(localPos[0],localPos[2],cherenkovLoss,idvol,kPhoton);//for photons in CsI 
                                    
                    if (fNsdigits > (Int_t)ckovData[8]) {
                        ckovData[8]= ckovData[8]+1;
@@ -1291,3 +1299,1853 @@ void AliRICHv3::StepManager()
       }//is MIP?
     /*************************************************End of MIP treatment**************************************/
 }//void AliRICHv3::StepManager()
+//__________________________________________________________________________________________________
+Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
+{//calls the charge disintegration method of the current chamber and adds all generated sdigits to the list of digits
+   
+   Int_t iChamber=kBad,iPadX=kBad,iPadY=kBad,iAdc=kBad,iTrack=kBad;
+   Float_t list[4][500];
+   Int_t iNdigits;
+  ((AliRICHChamber*)fChambers->At(idvol))->DisIntegration(eloss, xhit, yhit, iNdigits, list, res);
+    Int_t ic=0;
+    
+  for(Int_t i=0; i<iNdigits; i++) {
+    if(Int_t(list[0][i]) > 0) {
+           ic++;
+           iAdc = Int_t(list[0][i]);
+           iPadX = Int_t(list[1][i]);
+           iPadY = Int_t(list[2][i]);
+           iChamber = Int_t(list[3][i]);
+           AddSDigit(iChamber,iPadX,iPadY,iAdc,iTrack);
+       }
+  }
+    
+  if(fLoader->TreeS()){
+    fLoader->TreeS()->Fill();
+    fLoader->WriteSDigits("OVERWRITE");
+  }
+   return iNdigits;
+}//Int_t AliRICHv3::Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res)
+//__________________________________________________________________________________________________
+void AliRICHv3::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
+{
+  
+  Int_t NpadX = 162;                 // number of pads on X
+  Int_t NpadY = 162;                 // number of pads on Y
+  
+  Int_t Pad[162][162];
+  for (Int_t i=0;i<NpadX;i++) {
+    for (Int_t j=0;j<NpadY;j++) {
+      Pad[i][j]=0;
+    }
+  }
+  
+  //  Create some histograms
+
+  TH1F *pionspectra1 = new TH1F("pionspectra1","Pion Spectra",200,-4,2);
+  TH1F *pionspectra2 = new TH1F("pionspectra2","Pion Spectra",200,-4,2);
+  TH1F *pionspectra3 = new TH1F("pionspectra3","Pion Spectra",200,-4,2);
+  TH1F *protonspectra1 = new TH1F("protonspectra1","Proton Spectra",200,-4,2);
+  TH1F *protonspectra2 = new TH1F("protonspectra2","Proton Spectra",200,-4,2);
+  TH1F *protonspectra3 = new TH1F("protonspectra3","Proton Spectra",200,-4,2);
+  TH1F *kaonspectra1 = new TH1F("kaonspectra1","Kaon Spectra",100,-4,2);
+  TH1F *kaonspectra2 = new TH1F("kaonspectra2","Kaon Spectra",100,-4,2);
+  TH1F *kaonspectra3 = new TH1F("kaonspectra3","Kaon Spectra",100,-4,2);
+  TH1F *electronspectra1 = new TH1F("electronspectra1","Electron Spectra",100,-4,2);
+  TH1F *electronspectra2 = new TH1F("electronspectra2","Electron Spectra",100,-4,2);
+  TH1F *electronspectra3 = new TH1F("electronspectra3","Electron Spectra",100,-4,2);
+  TH1F *muonspectra1 = new TH1F("muonspectra1","Muon Spectra",100,-4,2);
+  TH1F *muonspectra2 = new TH1F("muonspectra2","Muon Spectra",100,-4,2);
+  TH1F *muonspectra3 = new TH1F("muonspectra3","Muon Spectra",100,-4,2);
+  TH1F *neutronspectra1 = new TH1F("neutronspectra1","Neutron Spectra",100,-4,2);
+  TH1F *neutronspectra2 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
+  TH1F *neutronspectra3 = new TH1F("neutronspectra2","Neutron Spectra",100,-4,2);
+  TH1F *chargedspectra1 = new TH1F("chargedspectra1","Charged particles above 1 GeV Spectra",100,-1,3);
+  TH1F *chargedspectra2 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
+  TH1F *chargedspectra3 = new TH1F("chargedspectra2","Charged particles above 1 GeV Spectra",100,-1,3);
+  TH1F *pionptspectrafinal = new TH1F("pionptspectrafinal","Primary Pions Transverse Momenta at HMPID",20,0,5);
+  TH1F *pionptspectravertex = new TH1F("pionptspectravertex","Primary Pions Transverse Momenta at vertex",20,0,5);
+  TH1F *kaonptspectrafinal = new TH1F("kaonptspectrafinal","Primary Kaons Transverse Momenta at HMPID",20,0,5);
+  TH1F *kaonptspectravertex = new TH1F("kaonptspectravertex","Primary Kaons Transverse Momenta at vertex",20,0,5);
+  //TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",100,-180,180);
+  TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of Theta angle of incidence, all tracks",100,0,50);
+  TH1F *hitsTheta500MeV = new TH1F("hitsTheta500MeV","Distribution of Theta angle of incidence, 0.5-1 GeV primary tracks",100,0,50);
+  TH1F *hitsTheta1GeV = new TH1F("hitsTheta1GeV","Distribution of Theta angle of incidence, 1-2 GeV primary tracks",100,0,50);
+  TH1F *hitsTheta2GeV = new TH1F("hitsTheta2GeV","Distribution of Theta angle of incidence, 2-3 GeV primary tracks",100,0,50);
+  TH1F *hitsTheta3GeV = new TH1F("hitsTheta3GeV","Distribution of Theta angle of incidence, >3 GeV primary tracks",100,0,50);
+  TH2F *production = new TH2F("production","Mother production vertices",100,-300,300,100,0,600);
+   
+   
+   
+
+//   Start loop over events 
+
+  Int_t pion=0, kaon=0, proton=0, electron=0, positron=0, neutron=0, highneutrons=0, muon=0;
+  Int_t chargedpions=0,primarypions=0,highprimarypions=0,chargedkaons=0,primarykaons=0,highprimarykaons=0;
+  Int_t photons=0, primaryphotons=0, highprimaryphotons=0;
+  TRandom* random=0;
+
+   for (int nev=0; nev<= evNumber2; nev++) {
+       Int_t nparticles = gAlice->GetEvent(nev);
+       
+
+       if (nev < evNumber1) continue;
+       if (nparticles <= 0) return;
+       
+// Get pointers to RICH detector and Hits containers
+       
+       AliRICH *pRICH = (AliRICH *) gAlice->GetDetector("RICH");
+     
+       TTree *treeH = TreeH();
+       Int_t ntracks =(Int_t) treeH->GetEntries();
+            
+// Start loop on tracks in the hits containers
+       
+       for (Int_t track=0; track<ntracks;track++) {
+          printf ("Processing Track: %d\n",track);
+          gAlice->ResetHits();
+          treeH->GetEvent(track);
+                          
+          for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1); 
+              mHit;
+              mHit=(AliRICHhit*)pRICH->NextHit()) 
+            {
+              //Int_t nch  = mHit->fChamber;              // chamber number
+              //Float_t x  = mHit->X();                    // x-pos of hit
+              //Float_t y  = mHit->Z();                    // y-pos
+              //Float_t z  = mHit->Y();
+              //Float_t phi = mHit->Phi();                 //Phi angle of incidence
+              Float_t theta = mHit->Theta();             //Theta angle of incidence
+              Float_t px = mHit->MomX();
+              Float_t py = mHit->MomY();
+              Int_t index = mHit->Track();
+              Int_t particle = (Int_t)(mHit->Particle());    
+              Float_t R;
+              Float_t PTfinal;
+              Float_t PTvertex;
+
+             TParticle *current = gAlice->Particle(index);
+             
+             //Float_t energy=current->Energy(); 
+
+             R=TMath::Sqrt(current->Vx()*current->Vx() + current->Vy()*current->Vy());
+             PTfinal=TMath::Sqrt(px*px + py*py);
+             PTvertex=TMath::Sqrt(current->Px()*current->Px() + current->Py()*current->Py());
+             
+             
+
+             if (TMath::Abs(particle) < 10000000)
+               {
+                 hitsTheta->Fill(theta,(float) 1);
+                 if (R<5)
+                   {
+                     if (PTvertex>.5 && PTvertex<=1)
+                       {
+                         hitsTheta500MeV->Fill(theta,(float) 1);
+                       }
+                     if (PTvertex>1 && PTvertex<=2)
+                       {
+                         hitsTheta1GeV->Fill(theta,(float) 1);
+                       }
+                     if (PTvertex>2 && PTvertex<=3)
+                       {
+                         hitsTheta2GeV->Fill(theta,(float) 1);
+                       }
+                     if (PTvertex>3)
+                       {
+                         hitsTheta3GeV->Fill(theta,(float) 1);
+                       }
+                   }
+                 
+               }
+
+             //if (nch == 3)
+               //{
+             
+             if (TMath::Abs(particle) < 50000051)
+               {
+                 //if (TMath::Abs(particle) == 50000050 || TMath::Abs(particle) == 2112)
+                 if (TMath::Abs(particle) == 2112 || TMath::Abs(particle) == 50000050)
+                   {
+                     //gMC->Rndm(&random, 1);
+                     if (random->Rndm() < .1)
+                       production->Fill(current->Vz(),R,(float) 1);
+                     if (TMath::Abs(particle) == 50000050)
+                       //if (TMath::Abs(particle) > 50000000)
+                       {
+                         photons +=1;
+                         if (R<5)
+                           {
+                             primaryphotons +=1;
+                             if (current->Energy()>0.001)
+                               highprimaryphotons +=1;
+                           }
+                       }       
+                     if (TMath::Abs(particle) == 2112)
+                       {
+                         neutron +=1;
+                         if (current->Energy()>0.0001)
+                           highneutrons +=1;
+                       }
+                   }
+                 if (TMath::Abs(particle) < 50000000)
+                   {
+                     production->Fill(current->Vz(),R,(float) 1);
+                   }
+                 //mip->Fill(x,y,(float) 1);
+               }
+             
+             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
+               {
+                 if (R<5)
+                   {
+                     pionptspectravertex->Fill(PTvertex,(float) 1);
+                     pionptspectrafinal->Fill(PTfinal,(float) 1);
+                   }
+               }
+             
+             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
+                 || TMath::Abs(particle)==311)
+               {
+                 if (R<5)
+                   {
+                     kaonptspectravertex->Fill(PTvertex,(float) 1);
+                     kaonptspectrafinal->Fill(PTfinal,(float) 1);
+                   }
+               }
+             
+             
+             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
+               {
+                 pionspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   pionspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   {
+                     pionspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                   }
+                 pion +=1;
+                 if (TMath::Abs(particle)==211)
+                   {
+                     chargedpions +=1;
+                     if (R<5)
+                       {
+                         primarypions +=1;
+                         if (current->Energy()>1)
+                           highprimarypions +=1;
+                       }
+                   }   
+               }
+             if (TMath::Abs(particle)==2212)
+               {
+                 protonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 //ptspectra->Fill(Pt,(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   protonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   protonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 proton +=1;
+               }
+             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
+                 || TMath::Abs(particle)==311)
+               {
+                 kaonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 //ptspectra->Fill(Pt,(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   kaonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   kaonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 kaon +=1;
+                 if (TMath::Abs(particle)==321)
+                   {
+                     chargedkaons +=1;
+                     if (R<5)
+                       {
+                         primarykaons +=1;
+                         if (current->Energy()>1)
+                           highprimarykaons +=1;
+                       }
+                   }
+               }
+             if (TMath::Abs(particle)==11)
+               {
+                 electronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 //ptspectra->Fill(Pt,(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   electronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   electronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (particle == 11)
+                   electron +=1;
+                 if (particle == -11)
+                   positron +=1;
+               }
+             if (TMath::Abs(particle)==13)
+               {
+                 muonspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 //ptspectra->Fill(Pt,(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   muonspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   muonspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 muon +=1;
+               }
+             if (TMath::Abs(particle)==2112)
+               {
+                 neutronspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 //ptspectra->Fill(Pt,(float) 1);
+                 if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                   neutronspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                 if (R>250 && R<450)
+                   {
+                     neutronspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                   }
+                 neutron +=1;
+               }
+             if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
+               {
+                 if (current->Energy()-current->GetCalcMass()>1)
+                   {
+                     chargedspectra1->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                     if (current->Vx()>5 && current->Vy()>5 && current->Vz()>5)
+                       chargedspectra2->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                     if (R>250 && R<450)
+                       chargedspectra3->Fill(TMath::Log10(current->Energy() - current->GetCalcMass()),(float) 1);
+                   }
+               }
+             // Fill the histograms
+             //Nh1+=nhits;
+             //h->Fill(x,y,(float) 1);
+             //}
+             //}
+          }          
+          
+       }
+       
+   }
+   //   }
+
+   TStyle *mystyle=new TStyle("Plain","mystyle");
+   mystyle->SetPalette(1,0);
+   mystyle->cd();
+   
+   //Create canvases, set the view range, show histograms
+
+    TCanvas *c2 = new TCanvas("c2","Angles of incidence",150,150,100,150);
+    c2->Divide(2,2);
+    //c2->SetFillColor(42);
+    
+    c2->cd(1);
+    hitsTheta500MeV->SetFillColor(5);
+    hitsTheta500MeV->Draw();
+    c2->cd(2);
+    hitsTheta1GeV->SetFillColor(5);
+    hitsTheta1GeV->Draw();
+    c2->cd(3);
+    hitsTheta2GeV->SetFillColor(5);
+    hitsTheta2GeV->Draw();
+    c2->cd(4);
+    hitsTheta3GeV->SetFillColor(5);
+    hitsTheta3GeV->Draw();
+    
+            
+   
+    TCanvas *c15 = new TCanvas("c15","Mothers Production Vertices",50,50,600,600);
+    c15->cd();
+    production->SetFillColor(42);
+    production->SetXTitle("z (m)");
+    production->SetYTitle("R (m)");
+    production->Draw();
+
+    TCanvas *c10 = new TCanvas("c10","Pt Spectra",50,50,600,700);
+    c10->Divide(2,2);
+    c10->cd(1);
+    pionptspectravertex->SetFillColor(5);
+    pionptspectravertex->SetXTitle("Pt (GeV)");
+    pionptspectravertex->Draw();
+    c10->cd(2);
+    pionptspectrafinal->SetFillColor(5);
+    pionptspectrafinal->SetXTitle("Pt (GeV)");
+    pionptspectrafinal->Draw();
+    c10->cd(3);
+    kaonptspectravertex->SetFillColor(5);
+    kaonptspectravertex->SetXTitle("Pt (GeV)");
+    kaonptspectravertex->Draw();
+    c10->cd(4);
+    kaonptspectrafinal->SetFillColor(5);
+    kaonptspectrafinal->SetXTitle("Pt (GeV)");
+    kaonptspectrafinal->Draw();
+   
+  
+   TCanvas *c16 = new TCanvas("c16","Particles Spectra II",150,150,600,350);
+   c16->Divide(2,1);
+   
+   c16->cd(1);
+   //TCanvas *c13 = new TCanvas("c13","Electron Spectra",400,10,600,700);
+   electronspectra1->SetFillColor(5);
+   electronspectra1->SetXTitle("log(GeV)");
+   electronspectra2->SetFillColor(46);
+   electronspectra2->SetXTitle("log(GeV)");
+   electronspectra3->SetFillColor(10);
+   electronspectra3->SetXTitle("log(GeV)");
+   //c13->SetLogx();
+   electronspectra1->Draw();
+   electronspectra2->Draw("same");
+   electronspectra3->Draw("same");
+   
+   c16->cd(2);
+   //TCanvas *c14 = new TCanvas("c14","Muon Spectra",400,10,600,700);
+   muonspectra1->SetFillColor(5);
+   muonspectra1->SetXTitle("log(GeV)");
+   muonspectra2->SetFillColor(46);
+   muonspectra2->SetXTitle("log(GeV)");
+   muonspectra3->SetFillColor(10);
+   muonspectra3->SetXTitle("log(GeV)");
+   //c14->SetLogx();
+   muonspectra1->Draw();
+   muonspectra2->Draw("same");
+   muonspectra3->Draw("same");
+   
+   //c16->cd(3);
+   //TCanvas *c16 = new TCanvas("c16","Neutron Spectra",400,10,600,700);
+   //neutronspectra1->SetFillColor(42);
+   //neutronspectra1->SetXTitle("log(GeV)");
+   //neutronspectra2->SetFillColor(46);
+   //neutronspectra2->SetXTitle("log(GeV)");
+   //neutronspectra3->SetFillColor(10);
+   //neutronspectra3->SetXTitle("log(GeV)");
+   //c16->SetLogx();
+   //neutronspectra1->Draw();
+   //neutronspectra2->Draw("same");
+   //neutronspectra3->Draw("same");
+
+   TCanvas *c9 = new TCanvas("c9","Particles Spectra",150,150,600,700);
+   //TCanvas *c9 = new TCanvas("c9","Pion Spectra",400,10,600,700);
+   c9->Divide(2,2);
+   
+   c9->cd(1);
+   pionspectra1->SetFillColor(5);
+   pionspectra1->SetXTitle("log(GeV)");
+   pionspectra2->SetFillColor(46);
+   pionspectra2->SetXTitle("log(GeV)");
+   pionspectra3->SetFillColor(10);
+   pionspectra3->SetXTitle("log(GeV)");
+   //c9->SetLogx();
+   pionspectra1->Draw();
+   pionspectra2->Draw("same");
+   pionspectra3->Draw("same");
+   
+   c9->cd(2);
+   //TCanvas *c10 = new TCanvas("c10","Proton Spectra",400,10,600,700);
+   protonspectra1->SetFillColor(5);
+   protonspectra1->SetXTitle("log(GeV)");
+   protonspectra2->SetFillColor(46);
+   protonspectra2->SetXTitle("log(GeV)");
+   protonspectra3->SetFillColor(10);
+   protonspectra3->SetXTitle("log(GeV)");
+   //c10->SetLogx();
+   protonspectra1->Draw();
+   protonspectra2->Draw("same");
+   protonspectra3->Draw("same");
+   
+   c9->cd(3);
+   //TCanvas *c11 = new TCanvas("c11","Kaon Spectra",400,10,600,700); 
+   kaonspectra1->SetFillColor(5);
+   kaonspectra1->SetXTitle("log(GeV)");
+   kaonspectra2->SetFillColor(46);
+   kaonspectra2->SetXTitle("log(GeV)");
+   kaonspectra3->SetFillColor(10);
+   kaonspectra3->SetXTitle("log(GeV)");
+   //c11->SetLogx();
+   kaonspectra1->Draw();
+   kaonspectra2->Draw("same");
+   kaonspectra3->Draw("same");
+   
+   c9->cd(4);
+   //TCanvas *c12 = new TCanvas("c12","Charged Particles Spectra",400,10,600,700);
+   chargedspectra1->SetFillColor(5);
+   chargedspectra1->SetXTitle("log(GeV)");
+   chargedspectra2->SetFillColor(46);
+   chargedspectra2->SetXTitle("log(GeV)");
+   chargedspectra3->SetFillColor(10);
+   chargedspectra3->SetXTitle("log(GeV)");
+   //c12->SetLogx();
+   chargedspectra1->Draw();
+   chargedspectra2->Draw("same");
+   chargedspectra3->Draw("same");
+   
+
+
+   printf("*****************************************\n");
+   printf("* Particle                   *  Counts  *\n");
+   printf("*****************************************\n");
+
+   printf("* Pions:                     *   %4d   *\n",pion);
+   printf("* Charged Pions:             *   %4d   *\n",chargedpions);
+   printf("* Primary Pions:             *   %4d   *\n",primarypions);
+   printf("* Primary Pions (p>1GeV/c):  *   %4d   *\n",highprimarypions);
+   printf("* Kaons:                     *   %4d   *\n",kaon);
+   printf("* Charged Kaons:             *   %4d   *\n",chargedkaons);
+   printf("* Primary Kaons:             *   %4d   *\n",primarykaons);
+   printf("* Primary Kaons (p>1GeV/c):  *   %4d   *\n",highprimarykaons);
+   printf("* Muons:                     *   %4d   *\n",muon);
+   printf("* Electrons:                 *   %4d   *\n",electron);
+   printf("* Positrons:                 *   %4d   *\n",positron);
+   printf("* Protons:                   *   %4d   *\n",proton);
+   printf("* All Charged:               *   %4d   *\n",(chargedpions+chargedkaons+muon+electron+positron+proton));
+   printf("*****************************************\n");
+   //printf("* Photons:                   *   %3.1f   *\n",photons); 
+   //printf("* Primary Photons:           *   %3.1f   *\n",primaryphotons);
+   //printf("* Primary Photons (p>1MeV/c):*   %3.1f   *\n",highprimaryphotons);
+   //printf("*****************************************\n");
+   //printf("* Neutrons:                  *   %3.1f   *\n",neutron);
+   //printf("* Neutrons (p>100keV/c):     *   %3.1f   *\n",highneutrons);
+   //printf("*****************************************\n");
+
+   if (gAlice->TreeD())
+     {
+       gAlice->TreeD()->GetEvent(0);
+   
+       Float_t occ[7]; 
+       Float_t sum=0;
+       Float_t mean=0; 
+       printf("\n*****************************************\n");
+       printf("* Chamber   * Digits      * Occupancy   *\n");
+       printf("*****************************************\n");
+       
+       for (Int_t ich=0;ich<7;ich++)
+        {
+          TClonesArray *Digits = DigitsAddress(ich);    //  Raw clusters branch
+          Int_t ndigits = Digits->GetEntriesFast();
+          occ[ich] = Float_t(ndigits)/(160*144);
+          sum += Float_t(ndigits)/(160*144);
+          printf("*   %d      *    %d      *   %3.1f%%     *\n",ich,ndigits,occ[ich]*100);
+        }
+       mean = sum/7;
+       printf("*****************************************\n");
+       printf("* Mean occupancy          *   %3.1f%%     *\n",mean*100);
+       printf("*****************************************\n");
+     }
+  printf("\nEnd of analysis\n");
+   
+}//void AliRICHv3::DiagnosticsFE(Int_t evNumber1,Int_t evNumber2)
+//__________________________________________________________________________________________________
+void AliRICHv3::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
+{
+
+AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
+   AliRICHSegmentationV0*  segmentation;
+   AliRICHChamber*       chamber;
+   
+   chamber = &(pRICH->Chamber(0));
+   segmentation=(AliRICHSegmentationV0*) chamber->GetSegmentationModel();
+
+   Int_t NpadX = segmentation->Npx();                 // number of pads on X
+   Int_t NpadY = segmentation->Npy();                 // number of pads on Y
+    
+   Int_t xmin= -NpadX/2;  
+   Int_t xmax=  NpadX/2;
+   Int_t ymin= -NpadY/2;
+   Int_t ymax=  NpadY/2;
+
+   Float_t PTfinal = 0;
+   Int_t pionCount = 0;
+   Int_t kaonCount = 0;
+   Int_t protonCount = 0;
+   
+   TH2F *feedback = 0;
+   TH2F *mip = 0;
+   TH2F *cerenkov = 0;
+   TH2F *h = 0;
+   TH1F *hitsX = 0;
+   TH1F *hitsY = 0;
+
+   TH2F *hc0 = new TH2F("hc0","Zoom on center of central chamber",150,-25,25,150,-45,5);
+
+   if (diaglevel == 1)
+     {
+       printf("Single Ring Hits\n");
+       feedback = new TH2F("feedback","Feedback hit distribution",150,-20,20,150,-35,5);
+       mip = new TH2F("mip","Mip hit distribution",150,-20,20,150,-35,5);
+       cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-20,20,150,-35,5);
+       h = new TH2F("h","Detector hit distribution",150,-20,20,150,-35,5);
+       hitsX = new TH1F("hitsX","Distribution of hits along x-axis",150,-50,50);
+       hitsY = new TH1F("hitsY","Distribution of hits along z-axis",150,-50,50);
+     }       
+   else
+     {
+       printf("Full Event Hits\n");
+       
+       feedback = new TH2F("feedback","Feedback hit distribution",150,-300,300,150,-300,300);
+       mip = new TH2F("mip","Mip hit distribution",150,-300,300,150,-300,300);
+       cerenkov = new TH2F("cerenkov","Cerenkov hit distribution",150,-300,300,150,-300,300);
+       h = new TH2F("h","Detector hit distribution",150,-300,300,150,-300,300); 
+       hitsX = new TH1F("digitsX","Distribution of hits along x-axis",200,-300,300);
+       hitsY = new TH1F("digitsY","Distribution of hits along z-axis",200,-300,300);
+     }
+   
+
+
+   TH2F *hc1 = new TH2F("hc1","Chamber 1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc2 = new TH2F("hc2","Chamber 2 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc3 = new TH2F("hc3","Chamber 3 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc4 = new TH2F("hc4","Chamber 4 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc5 = new TH2F("hc5","Chamber 5 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc6 = new TH2F("hc6","Chamber 6 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+   TH2F *hc7 = new TH2F("hc7","Chamber 7 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
+      
+   TH1F *Clcharge = new TH1F("Clcharge","Cluster Charge Distribution",500,0.,500.);
+   TH1F *ckovangle = new TH1F("ckovangle","Cerenkov angle per photon",100,.35,.8);
+   TH1F *hckphi = new TH1F("hckphi","Cerenkov phi angle per photon",620,-3.1,3.1);
+   TH1F *mother = new TH1F("mother","Cerenkovs per Mip",75,0.,75.);
+   TH1F *radius = new TH1F("radius","Mean distance to Mip",100,0.,20.);
+   TH1F *phspectra1 = new TH1F("phspectra1","Detected Photon Spectra",200,5.,10.);
+   TH1F *phspectra2 = new TH1F("phspectra2","Produced Photon Spectra",200,5.,10.);
+   TH1F *totalphotonstrack = new TH1F("totalphotonstrack","Produced Photons per Mip",100,200,700.);
+   TH1F *totalphotonsevent = new TH1F("totalphotonsevent","Produced Photons per Mip",100,200,700.);
+   //TH1F *feedbacks = new TH1F("feedbacks","Produced Feedbacks per Mip",50,0.5,50.);
+   TH1F *padnumber = new TH1F("padnumber","Number of pads per cluster",50,-0.5,50.);
+   TH1F *padsev = new TH1F("padsev","Number of pads hit per MIP",50,0.5,100.);
+   TH1F *clusev = new TH1F("clusev","Number of clusters per MIP",50,0.5,50.);
+   TH1F *photev = new TH1F("photev","Number of detected photons per MIP",50,0.5,50.);
+   TH1F *feedev = new TH1F("feedev","Number of feedbacks per MIP",50,0.5,50.);
+   TH1F *padsmip = new TH1F("padsmip","Number of pads per event inside MIP region",50,0.5,50.);
+   TH1F *padscl = new TH1F("padscl","Number of pads per event from cluster count",50,0.5,100.);
+   TH1F *pionspectra = new TH1F("pionspectra","Pion Spectra",200,.5,10.);
+   TH1F *protonspectra = new TH1F("protonspectra","Proton Spectra",200,.5,10.);
+   TH1F *kaonspectra = new TH1F("kaonspectra","Kaon Spectra",100,.5,10.);
+   TH1F *chargedspectra = new TH1F("chargedspectra","Charged particles above 1 GeV Spectra",100,.5,10.);
+   TH1F *hitsPhi = new TH1F("hitsPhi","Distribution of phi angle of incidence",50,0,360);
+   TH1F *hitsTheta = new TH1F("hitsTheta","Distribution of theta angle of incidence",50,0,15);
+   TH1F *Omega1D = new TH1F("omega","Reconstructed Cerenkov angle per track",50,.5,1);
+   TH1F *Theta = new TH1F("theta","Reconstructed theta incidence angle per track",100,0,15);
+   TH1F *Phi = new TH1F("phi","Reconstructed phi incidence per track",100,0,360);
+   TH1F *Omega3D = new TH1F("omega","Reconstructed Cerenkov angle per track",100,.35,.8);
+   TH1F *PhotonCer = new TH1F("photoncer","Reconstructed Cerenkov angle per photon",100,.35,.8);
+   TH2F *PadsUsed = new TH2F("padsused","Pads Used for Reconstruction",100,-30,30,100,-30,30);
+   TH1F *MeanRadius = new TH1F("radius","Mean Radius for reconstructed track",100,0.,20.);
+   TH2F *identification = new TH2F("identification","Particle Identification",100,1,5,100,0,.8);
+   TH1F *OriginalOmega = new TH1F("Original Omega","Cerenkov angle per track",100,.35,.8);
+   TH1F *OriginalPhi = new TH1F("Original Phi","Distribution of phi angle of incidence per track",100,0,360);
+   TH1F *OriginalTheta = new TH1F("Original Theta","Distribution of theta angle per track",100,0,15);
+   TH1F *OmegaError = new TH1F("Omega Error","Difference between original an reconstructed cerenkov angle",100,0,.2);
+   TH1F *PhiError = new TH1F("Phi Error","Difference between original an reconstructed phi angle",100,0,360);
+   TH1F *ThetaError = new TH1F("Theta Error","Difference between original an reconstructed phi angle",100,0,15);
+
+
+//   Start loop over events 
+
+   Int_t Nh=0;
+   Int_t pads=0;
+   Int_t Nh1=0;
+   Int_t mothers[80000];
+   Int_t mothers2[80000];
+   Float_t mom[3];
+   Int_t nraw=0;
+   Int_t phot=0;
+   Int_t feed=0;
+   Int_t padmip=0;
+   Float_t x=0,y=0;
+
+   Float_t chiSquareOmega = 0;
+   Float_t chiSquareTheta = 0;
+   Float_t chiSquarePhi = 0;
+
+   Float_t recEffEvent = 0;
+   Float_t recEffTotal = 0;
+
+   Float_t trackglob[3];
+   Float_t trackloc[3];
+
+   
+   for (Int_t i=0;i<100;i++) mothers[i]=0;
+
+   for (int nev=0; nev<= evNumber2; nev++) {
+       Int_t nparticles = gAlice->GetEvent(nev);
+       
+
+       //cout<<"nev  "<<nev<<endl;
+       printf ("\n**********************************\nProcessing Event: %d\n",nev);
+       //cout<<"nparticles  "<<nparticles<<endl;
+       printf ("Particles       : %d\n\n",nparticles);
+       if (nev < evNumber1) continue;
+       if (nparticles <= 0) return;
+       
+// Get pointers to RICH detector and Hits containers
+       
+
+       TTree *TH = TreeH(); 
+       Stat_t ntracks = TH->GetEntries();
+
+       // Start loop on tracks in the hits containers
+       //Int_t Nc=0;
+       for (Int_t track=0; track<ntracks;track++) {
+          
+        printf ("\nProcessing Track: %d\n",track);
+        gAlice->ResetHits();
+        TH->GetEvent(track);
+        Int_t nhits = pRICH->Hits()->GetEntriesFast();
+        if (nhits) Nh+=nhits;
+        printf("Hits            : %d\n",nhits);
+        for(AliRICHhit* mHit=(AliRICHhit*)pRICH->FirstHit(-1); 
+            mHit;
+            mHit=(AliRICHhit*)pRICH->NextHit()) 
+          {
+            Int_t nch  = mHit->Chamber();              // chamber number
+            trackglob[0] = mHit->X();                 // x-pos of hit
+            trackglob[1] = mHit->Y();
+            trackglob[2] = mHit->Z();                 // y-pos of hit
+            //x  = mHit->X();                           // x-pos of hit
+            //y  = mHit->Z();                           // y-pos
+            Float_t phi = mHit->Phi();                 //Phi angle of incidence
+            Float_t theta = mHit->Theta();             //Theta angle of incidence
+            Int_t index = mHit->Track();
+            Int_t particle = (Int_t)(mHit->Particle());        
+            //Int_t freon = (Int_t)(mHit->fLoss);    
+            Float_t px = mHit->MomX();
+            Float_t py = mHit->MomY();
+            
+            if (TMath::Abs(particle) < 10000000)
+              {
+                PTfinal=TMath::Sqrt(px*px + py*py);
+              }
+       
+            chamber = &(pRICH->Chamber(nch-1));
+            
+            
+            chamber->GlobaltoLocal(trackglob,trackloc);
+            
+            chamber->LocaltoGlobal(trackloc,trackglob);
+            
+       
+            x=trackloc[0];
+            y=trackloc[2];
+            
+            hitsX->Fill(x,(float) 1);
+            hitsY->Fill(y,(float) 1);
+              
+             
+             TParticle *current = (TParticle*)gAlice->Particle(index);
+
+             hitsTheta->Fill(theta,(float) 1);
+            
+             if (current->GetPdgCode() < 10000000)
+               {
+                 mip->Fill(x,y,(float) 1);
+                 hitsPhi->Fill(TMath::Abs(phi),(float) 1);
+               }
+             
+             if (TMath::Abs(particle)==211 || TMath::Abs(particle)==111)
+               {
+                 pionspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
+               }
+             if (TMath::Abs(particle)==2212)
+               {
+                 protonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
+               }
+             if (TMath::Abs(particle)==321 || TMath::Abs(particle)==130 || TMath::Abs(particle)==310 
+                 || TMath::Abs(particle)==311)
+               {
+                 kaonspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
+               }
+             if(TMath::Abs(particle)==211 || TMath::Abs(particle)==2212 || TMath::Abs(particle)==321)
+               {
+                 if (current->Energy() - current->GetCalcMass()>1)
+                   chargedspectra->Fill(current->Energy() - current->GetCalcMass(),(float) 1);
+               }
+             //printf("Hits:%d\n",hit);
+             //printf ("Chamber number:%d x:%f y:%f\n",nch,x,y);
+              // Fill the histograms
+             Nh1+=nhits;
+             h->Fill(x,y,(float) 1);
+                 //}
+              //}
+          }
+          
+          Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
+          //if (current->GetPdgCode() < 50000051 && current->GetPdgCode() > 50000040)
+          //totalphotonsevent->Fill(ncerenkovs,(float) 1);
+
+          if (ncerenkovs) {
+            printf("Cerenkovs       : %d\n",ncerenkovs);
+            totalphotonsevent->Fill(ncerenkovs,(float) 1);
+            for (Int_t hit=0;hit<ncerenkovs;hit++) {
+              AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
+              Int_t nchamber = cHit->fChamber;     // chamber number
+              Int_t index =    cHit->Track();
+              //Int_t pindex =   (Int_t)(cHit->fIndex);
+              trackglob[0] = cHit->X();                 // x-pos of hit
+              trackglob[1] = cHit->Y();
+              trackglob[2] = cHit->Z();                 // y-pos of hit
+              //Float_t cx  =      cHit->X();                // x-position
+              //Float_t cy  =      cHit->Z();                // y-position
+              Int_t cmother =  cHit->fCMother;      // Index of mother particle
+              Int_t closs =    (Int_t)(cHit->fLoss);           // How did the particle get lost? 
+              Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
+              
+              chamber = &(pRICH->Chamber(nchamber-1));
+            
+              //printf("Nch:%d\n",nch);
+              
+              chamber->GlobaltoLocal(trackglob,trackloc);
+            
+              chamber->LocaltoGlobal(trackloc,trackglob);
+            
+       
+              Float_t cx=trackloc[0];
+              Float_t cy=trackloc[2];
+              
+              //printf ("Cerenkov hit number %d/%d, X:%f, Y:%f\n",hit,ncerenkovs,cx,cy); 
+
+
+              //printf("Particle:%9d\n",index);
+                                
+              TParticle *current = (TParticle*)gAlice->Particle(index);
+              Float_t energyckov = current->Energy();
+              
+              if (current->GetPdgCode() == 50000051)
+                {
+                  if (closs==4)
+                    {
+                      feedback->Fill(cx,cy,(float) 1);
+                      feed++;
+                    }
+                }
+              if (current->GetPdgCode() == 50000050)
+                {
+                  
+                  if (closs !=4)
+                    {
+                      phspectra2->Fill(energyckov*1e9,(float) 1);
+                    }
+                      
+                  if (closs==4)
+                    {
+                      cerenkov->Fill(cx,cy,(float) 1); 
+                      
+                      //printf ("Cerenkov hit number %d/%d, X:%d, Y:%d\n",hit,ncerenkovs,cx,cy); 
+                      
+                      //TParticle *MIP = (TParticle*)gAlice->Particle(cmother);
+                      AliRICHhit* mipHit = (AliRICHhit*) pRICH->Hits()->UncheckedAt(0);
+                      mom[0] = current->Px();
+                      mom[1] = current->Py();
+                      mom[2] = current->Pz();
+                      //mom[0] = cHit->fMomX;
+                      // mom[1] = cHit->fMomZ;
+                      //mom[2] = cHit->fMomY;
+                      //Float_t energymip = MIP->Energy();
+                      //Float_t Mip_px = mipHit->fMomFreoX;
+                      //Float_t Mip_py = mipHit->fMomFreoY;
+                      //Float_t Mip_pz = mipHit->fMomFreoZ;
+                      //Float_t Mip_px = MIP->Px();
+                      //Float_t Mip_py = MIP->Py();
+                      //Float_t Mip_pz = MIP->Pz();
+                      
+                      
+                      
+                      //Float_t r = mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2];
+                      //Float_t rt = TMath::Sqrt(r);
+                      //Float_t Mip_r = Mip_px*Mip_px + Mip_py*Mip_py + Mip_pz*Mip_pz; 
+                      //Float_t Mip_rt = TMath::Sqrt(Mip_r);
+                      //Float_t coscerenkov = (mom[0]*Mip_px + mom[1]*Mip_py + mom[2]*Mip_pz)/(rt*Mip_rt+0.0000001);
+                      //Float_t cherenkov = TMath::ACos(coscerenkov);
+                      ckovangle->Fill(cherenkov,(float) 1);                           //Cerenkov angle calculus
+                      //printf("Cherenkov: %f\n",cherenkov);
+                      Float_t ckphi=TMath::ATan2(mom[0], mom[2]);
+                      hckphi->Fill(ckphi,(float) 1);
+                      
+                      
+                      //Float_t mix = MIP->Vx();
+                      //Float_t miy = MIP->Vy();
+                      Float_t mx = mipHit->X();
+                      Float_t my = mipHit->Z();
+                      //printf("FX %e, FY %e, VX %e, VY %e\n",cx,cy,mx,my);
+                      Float_t dx = trackglob[0] - mx;
+                      Float_t dy = trackglob[2] - my;
+                      //printf("Dx:%f, Dy:%f\n",dx,dy);
+                      Float_t final_radius = TMath::Sqrt(dx*dx+dy*dy);
+                      //printf("Final radius:%f\n",final_radius);
+                      radius->Fill(final_radius,(float) 1);
+                      
+                      phspectra1->Fill(energyckov*1e9,(float) 1);
+                      phot++;
+                    }
+                  for (Int_t nmothers=0;nmothers<=ntracks;nmothers++){
+                    if (cmother == nmothers){
+                      if (closs == 4)
+                        mothers2[cmother]++;
+                      mothers[cmother]++;
+                    }
+                  } 
+                }
+            }
+          }
+          
+
+          if(gAlice->TreeR())
+            {
+              Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
+              gAlice->TreeR()->GetEvent(nent-1);
+              TClonesArray *Rawclusters = pRICH->RawClustAddress(2);    //  Raw clusters branch
+              //printf ("Rawclusters:%p",Rawclusters);
+              Int_t nrawclusters = Rawclusters->GetEntriesFast();
+                      
+              if (nrawclusters) {
+                printf("Raw Clusters    : %d\n",nrawclusters);
+                for (Int_t hit=0;hit<nrawclusters;hit++) {
+                  AliRICHRawCluster* rcHit = (AliRICHRawCluster*) pRICH->RawClustAddress(2)->UncheckedAt(hit);
+                  //Int_t nchamber = rcHit->fChamber;     // chamber number
+                  //Int_t nhit = cHit->fHitNumber;        // hit number
+                  Int_t qtot = rcHit->fQ;                 // charge
+                  Float_t fx  =  rcHit->fX;                 // x-position
+                  Float_t fy  =  rcHit->fY;                 // y-position
+                  //Int_t type = rcHit->fCtype;             // cluster type ?   
+                  Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
+                  pads += mult;
+                  if (qtot > 0) {
+                    //printf ("fx: %d, fy: %d\n",fx,fy);
+                    if (fx>(x-4) && fx<(x+4)  && fy>(y-4) && fy<(y+4)) {
+                      //printf("There %d \n",mult);
+                      padmip+=mult;
+                    } else {
+                      padnumber->Fill(mult,(float) 1);
+                      nraw++;
+                      if (mult<4) Clcharge->Fill(qtot,(float) 1);
+                    }
+                    
+                  }
+                }
+              }
+              
+              
+              TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
+              Int_t nrechits1D = RecHits1D->GetEntriesFast();
+              //printf (" nrechits:%d\n",nrechits);
+              
+              if(nrechits1D)
+                {
+                  for (Int_t hit=0;hit<nrechits1D;hit++) {
+                    AliRICHRecHit1D* recHit1D = (AliRICHRecHit1D*) pRICH->RecHitsAddress1D(2)->UncheckedAt(hit);
+                    Float_t r_omega = recHit1D->fOmega;                  // Cerenkov angle
+                    Float_t *cer_pho = recHit1D->fCerPerPhoton;        // Cerenkov angle per photon
+                    Int_t *padsx = recHit1D->fPadsUsedX;           // Pads Used fo reconstruction (x)
+                    Int_t *padsy = recHit1D->fPadsUsedY;           // Pads Used fo reconstruction (y)
+                    Int_t goodPhotons = recHit1D->fGoodPhotons;    // Number of pads used for reconstruction
+                    
+                    Omega1D->Fill(r_omega,(float) 1);
+                    
+                    for (Int_t i=0; i<goodPhotons; i++)
+                      {
+                        PhotonCer->Fill(cer_pho[i],(float) 1);
+                        PadsUsed->Fill(padsx[i],padsy[i],1);
+                        //printf("Angle:%f, pad: %d %d\n",cer_pho[i],padsx[i],padsy[i]);
+                      }
+                    
+                    //printf("Omega: %f, Theta: %f, Phi: %f\n",r_omega,r_theta,r_phi);
+                  }
+                }
+
+              
+              TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
+              Int_t nrechits3D = RecHits3D->GetEntriesFast();
+              //printf (" nrechits:%d\n",nrechits);
+              
+              if(nrechits3D)
+                {
+                  recEffEvent = 0;
+                  
+                  //for (Int_t hit=0;hit<nrechits3D;hit++) {
+                  AliRICHRecHit3D* recHit3D = (AliRICHRecHit3D*) pRICH->RecHitsAddress3D(2)->UncheckedAt(track);
+                  Float_t r_omega    = recHit3D->fOmega;                  // Cerenkov angle
+                  Float_t r_theta    = recHit3D->fTheta;                  // Theta angle of incidence
+                  Float_t r_phi      = recHit3D->fPhi;                    // Phi angle if incidence
+                  Float_t meanradius = recHit3D->fMeanRadius;              // Mean radius for reconstructed point
+                  Float_t originalOmega = recHit3D->fOriginalOmega;       // Real Cerenkov angle
+                  Float_t originalTheta = recHit3D->fOriginalTheta;       // Real incidence angle
+                  Float_t originalPhi = recHit3D->fOriginalPhi;           // Real azimuthal angle
+                  
+                  
+                  //correction to track cerenkov angle
+                  originalOmega = (Float_t) ckovangle->GetMean();
+                  
+                  if(diaglevel == 4)
+                    {
+                      printf("\nMean cerenkov angle: %f\n", originalOmega);
+                      printf("Reconstructed cerenkov angle: %f\n",r_omega);
+                    }
+                  
+                  Float_t omegaError = TMath::Abs(originalOmega - r_omega);
+                  Float_t thetaError = TMath::Abs(originalTheta - r_theta);
+                  Float_t phiError   = TMath::Abs(originalPhi - r_phi);
+                  
+                  
+                  if(TMath::Abs(omegaError) < 0.015)
+                    recEffEvent += 1;
+                  
+                  Omega3D->Fill(r_omega,(float) 1);
+                  Theta->Fill(r_theta*180/TMath::Pi(),(float) 1);
+                  Phi->Fill(r_phi*180/TMath::Pi()-180,(float) 1);
+                  MeanRadius->Fill(meanradius,(float) 1);
+                  identification->Fill(PTfinal, r_omega,1);
+                  OriginalOmega->Fill(originalOmega, (float) 1);
+                  OriginalTheta->Fill(originalTheta, (float) 1);
+                  OriginalPhi->Fill(TMath::Abs(originalPhi), (float) 1);
+                  OmegaError->Fill(omegaError, (float) 1);
+                  ThetaError->Fill(thetaError, (float) 1);
+                  PhiError->Fill(phiError, (float) 1);
+                  
+                  recEffEvent = recEffEvent;
+                  recEffTotal += recEffEvent;
+                  
+                  Float_t pioncer = acos(sqrt((.139*.139+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
+                  Float_t kaoncer = acos(sqrt((.439*.439+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
+                  Float_t protoncer = acos(sqrt((.938*.938+PTfinal*PTfinal)/(PTfinal*PTfinal*1.285*1.285)));
+
+                  Float_t piondist = TMath::Abs(r_omega - pioncer);
+                  Float_t kaondist = TMath::Abs(r_omega - kaoncer);
+                  Float_t protondist = TMath::Abs(r_omega - protoncer);
+
+                  if(diaglevel == 4)
+                    {
+                      if(pioncer<r_omega)
+                        {
+                          printf("Identified as a PION!\n");
+                          pionCount += 1;
+                        }
+                      if(kaoncer<r_omega && pioncer>r_omega)
+                        {
+                          if(kaondist>piondist)
+                            {
+                              printf("Identified as a PION!\n");
+                              pionCount += 1;
+                            }
+                          else
+                            {
+                              printf("Identified as a KAON!\n");
+                              kaonCount += 1;
+                            }
+                        }                       }
+                      if(protoncer<r_omega && kaoncer>r_omega)
+                        {
+                          if(kaondist>protondist)
+                            {
+                              printf("Identified as a PROTON!\n");
+                              protonCount += 1;
+                            }
+                          else
+                            {
+                              printf("Identified as a KAON!\n");
+                              pionCount += 1;
+                            }
+                        }
+                      if(protoncer>r_omega)
+                        {
+                          printf("Identified as a PROTON!\n");
+                          protonCount += 1;
+                        }
+
+                      printf("\nReconstruction efficiency: %5.2f%%\n", recEffEvent*100);
+                }
+            }
+       }
+   
+       
+       for (Int_t nmothers=0;nmothers<ntracks;nmothers++){
+        totalphotonstrack->Fill(mothers[nmothers],(float) 1);
+        mother->Fill(mothers2[nmothers],(float) 1);
+       }
+       
+       clusev->Fill(nraw,(float) 1);
+       photev->Fill(phot,(float) 1);
+       feedev->Fill(feed,(float) 1);
+       padsmip->Fill(padmip,(float) 1);
+       padscl->Fill(pads,(float) 1);
+       phot = 0;
+       feed = 0;
+       pads = 0;
+       nraw=0;
+       padmip=0;
+       
+       
+       
+       gAlice->ResetDigits();
+       gAlice->TreeD()->GetEvent(0);
+       
+       if (diaglevel < 4)
+        {
+          
+          
+          TClonesArray *Digits  = pRICH->DigitsAddress(2);
+          Int_t ndigits = Digits->GetEntriesFast();
+          printf("Digits          : %d\n",ndigits);
+          padsev->Fill(ndigits,(float) 1);
+          for (Int_t hit=0;hit<ndigits;hit++) {
+            AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
+            Int_t qtot = dHit->Signal();                // charge
+            Int_t ipx  = dHit->PadX();               // pad number on X
+            Int_t ipy  = dHit->PadY();               // pad number on Y
+            //printf("%d, %d\n",ipx,ipy);
+            if( ipx<=100 && ipy <=100) hc0->Fill(ipx,ipy,(float) qtot);
+          }
+        }
+       
+       if (diaglevel == 5)
+        {
+          for (Int_t ich=0;ich<7;ich++)
+            {
+              TClonesArray *Digits = pRICH->DigitsAddress(ich);    //  Raw clusters branch
+              Int_t ndigits = Digits->GetEntriesFast();
+              //printf("Digits:%d\n",ndigits);
+              padsev->Fill(ndigits,(float) 1); 
+              if (ndigits) {
+                for (Int_t hit=0;hit<ndigits;hit++) {
+                  AliRICHDigit* dHit = (AliRICHDigit*) Digits->UncheckedAt(hit);
+                  Int_t qtot = dHit->Signal();                // charge
+                  Int_t ipx  = dHit->PadX();               // pad number on X
+                  Int_t ipy  = dHit->PadY();               // pad number on Y
+                  if( ipx<=100 && ipy <=100 && ich==2) hc0->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==0) hc1->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==1) hc2->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==2) hc3->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==3) hc4->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==4) hc5->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==5) hc6->Fill(ipx,ipy,(float) qtot);
+                  if( ipx<=162 && ipy <=162 && ich==6) hc7->Fill(ipx,ipy,(float) qtot);
+                }
+              }
+            }
+        }
+   }
+   
+   if(diaglevel == 4)
+     {
+
+       Stat_t omegaE;
+       Stat_t thetaE;
+       Stat_t phiE;
+       
+       Stat_t omegaO;
+       Stat_t thetaO;
+       Stat_t phiO;
+       
+       for(Int_t i=0;i<99;i++)
+        {
+          omegaE = OriginalOmega->GetBinContent(i);
+          if(omegaE != 0)
+            {
+              omegaO = Omega3D->GetBinContent(i);
+              chiSquareOmega += (TMath::Power(omegaE,2) - TMath::Power(omegaO,2))/omegaO;
+            }
+
+          thetaE = OriginalTheta->GetBinContent(i);
+          if(thetaE != 0)
+            {
+              thetaO = Theta->GetBinContent(i);
+              chiSquareTheta += (TMath::Power(thetaE,2) - TMath::Power(thetaO,2))/thetaO;
+            }
+
+          phiE = OriginalPhi->GetBinContent(i);
+          if(phiE != 0)
+            {
+              phiO = Phi->GetBinContent(i);
+              chiSquarePhi += (TMath::Power(phiE,2) - TMath::Power(phiO,2))/phiO;
+            }
+        }
+
+       
+
+       printf("\nChi square test values:   Omega - %f\n", chiSquareOmega);
+       printf("                          Theta - %f\n", chiSquareTheta);
+       printf("                          Phi   - %f\n", chiSquarePhi);
+       
+       printf("\nKolmogorov test values:   Omega - %5.4f\n", Omega3D->KolmogorovTest(OriginalOmega));
+       printf("                          Theta - %5.4f\n", Theta->KolmogorovTest(OriginalTheta));
+       printf("                          Phi   - %5.4f\n", Phi->KolmogorovTest(OriginalPhi));
+
+       recEffTotal = recEffTotal/evNumber2;
+       printf("\nTotal reconstruction efficiency: %5.2f%%\n", recEffTotal*100);
+       printf("\n Pions: %d\n Kaons: %d\n Protons:%d\n",pionCount, kaonCount, protonCount);
+
+     }
+   
+   
+   //Create canvases, set the view range, show histograms
+
+   TCanvas *c1 = 0;
+   TCanvas *c2 = 0;
+   TCanvas *c3 = 0;
+   TCanvas *c4 = 0;
+   TCanvas *c5 = 0;
+   TCanvas *c6 = 0;
+   TCanvas *c7 = 0;
+   TCanvas *c8 = 0;
+   TCanvas *c9 = 0;
+   TCanvas *c10 = 0;
+   TCanvas *c11 = 0;
+   TCanvas *c12 = 0;
+   TCanvas *c13 = 0;
+
+   
+   TStyle *mystyle=new TStyle("Plain","mystyle");
+   mystyle->SetPalette(1,0);
+   mystyle->SetFuncColor(2);
+   mystyle->SetDrawBorder(0);
+   mystyle->SetTitleBorderSize(0);
+   mystyle->SetOptFit(1111);
+   mystyle->cd();
+
+   
+   TClonesArray *RecHits3D = pRICH->RecHitsAddress3D(2);
+   Int_t nrechits3D = RecHits3D->GetEntriesFast();
+   TClonesArray *RecHits1D = pRICH->RecHitsAddress1D(2);
+   Int_t nrechits1D = RecHits1D->GetEntriesFast();
+
+  switch(diaglevel)
+     {
+     case 1:
+       
+       c1 = new TCanvas("c1","Alice RICH digits",50,50,300,350);
+       hc0->SetXTitle("ix (npads)");
+       hc0->Draw("colz");
+       
+       c2 = new TCanvas("c2","Hits per type",100,100,600,700);
+       c2->Divide(2,2);
+       //c4->SetFillColor(42);
+
+       c2->cd(1);
+       feedback->SetXTitle("x (cm)");
+       feedback->SetYTitle("y (cm)");
+       feedback->Draw("colz");
+       
+       c2->cd(2);
+       //mip->SetFillColor(5);
+       mip->SetXTitle("x (cm)");
+       mip->SetYTitle("y (cm)");
+       mip->Draw("colz");
+       
+       c2->cd(3);
+       //cerenkov->SetFillColor(5);
+       cerenkov->SetXTitle("x (cm)");
+       cerenkov->SetYTitle("y (cm)"); 
+       cerenkov->Draw("colz");
+       
+       c2->cd(4);
+       //h->SetFillColor(5);
+       h->SetXTitle("x (cm)");
+       h->SetYTitle("y (cm)");
+       h->Draw("colz");
+
+       c3 = new TCanvas("c3","Hits distribution",150,150,600,350);
+       c3->Divide(2,1);
+       //c10->SetFillColor(42);
+       
+       c3->cd(1);
+       hitsX->SetFillColor(5);
+       hitsX->SetXTitle("(cm)");
+       hitsX->Draw();
+       
+       c3->cd(2);
+       hitsY->SetFillColor(5);
+       hitsY->SetXTitle("(cm)");
+       hitsY->Draw();
+       
+      
+       break;
+     case 2:
+       
+       c4 = new TCanvas("c4","Photon Spectra",50,50,600,350);
+       c4->Divide(2,1);
+       
+       c4->cd(1);
+       phspectra2->SetFillColor(5);
+       phspectra2->SetXTitle("energy (eV)");
+       phspectra2->Draw();
+       c4->cd(2);
+       phspectra1->SetFillColor(5);
+       phspectra1->SetXTitle("energy (eV)");
+       phspectra1->Draw();
+       
+       c5 = new TCanvas("c5","Particles Spectra",100,100,600,700);
+       c5->Divide(2,2);
+       
+       c5->cd(1);
+       pionspectra->SetFillColor(5);
+       pionspectra->SetXTitle("(GeV)");
+       pionspectra->Draw();
+       
+       c5->cd(2);
+       protonspectra->SetFillColor(5);
+       protonspectra->SetXTitle("(GeV)");
+       protonspectra->Draw();
+       
+       c5->cd(3);
+       kaonspectra->SetFillColor(5);
+       kaonspectra->SetXTitle("(GeV)");
+       kaonspectra->Draw();
+       
+       c5->cd(4);
+       chargedspectra->SetFillColor(5);
+       chargedspectra->SetXTitle("(GeV)");
+       chargedspectra->Draw();
+
+       break;
+       
+     case 3:
+
+       
+       if(gAlice->TreeR())
+        {
+          c6=new TCanvas("c6","Clusters Statistics",50,50,600,700);
+          c6->Divide(2,2);
+          
+          c6->cd(1);
+          Clcharge->SetFillColor(5);
+          Clcharge->SetXTitle("ADC counts");
+          if (evNumber2>10)
+            {
+              Clcharge->Fit("expo");
+            }
+          Clcharge->Draw();
+          
+          c6->cd(2);
+          padnumber->SetFillColor(5);
+          padnumber->SetXTitle("(counts)");
+          padnumber->Draw();
+          
+          c6->cd(3);
+          clusev->SetFillColor(5);
+          clusev->SetXTitle("(counts)");
+          if (evNumber2>10)
+            {
+              clusev->Fit("gaus");
+              //gaus->SetLineColor(2);
+              //gaus->SetLineWidth(3);
+            }
+          clusev->Draw();
+          
+          c6->cd(4);
+          padsmip->SetFillColor(5);
+          padsmip->SetXTitle("(counts)");
+          padsmip->Draw(); 
+        }
+       
+       if(evNumber2<1)
+        {
+          c11 = new TCanvas("c11","Cherenkov per Mip",400,10,600,700);
+          mother->SetFillColor(5);
+          mother->SetXTitle("counts");
+          mother->Draw();
+        }
+
+       c7 = new TCanvas("c7","Production Statistics",100,100,600,700);
+       c7->Divide(2,2);
+       //c7->SetFillColor(42);
+       
+       c7->cd(1);
+       totalphotonsevent->SetFillColor(5);
+       totalphotonsevent->SetXTitle("Photons (counts)");
+       if (evNumber2>10)
+          {
+            totalphotonsevent->Fit("gaus");
+            //gaus->SetLineColor(2);
+            //gaus->SetLineWidth(3);
+          }
+       totalphotonsevent->Draw();
+       
+       c7->cd(2);
+       photev->SetFillColor(5);
+       photev->SetXTitle("(counts)");
+       if (evNumber2>10)
+        {
+          photev->Fit("gaus");
+          //gaus->SetLineColor(2);
+          //gaus->SetLineWidth(3);
+        }
+       photev->Draw();
+       
+       c7->cd(3);
+       feedev->SetFillColor(5);
+       feedev->SetXTitle("(counts)");
+       if (evNumber2>10)
+        {
+          feedev->Fit("gaus");
+        }
+       feedev->Draw();
+
+       c7->cd(4);
+       padsev->SetFillColor(5);
+       padsev->SetXTitle("(counts)");
+       if (evNumber2>10)
+        {
+          padsev->Fit("gaus");
+        }
+       padsev->Draw();
+
+       break;
+
+     case 4:
+       
+
+       if(nrechits3D)
+        {
+          c8 = new TCanvas("c8","3D reconstruction of Phi angle",50,50,300,1050);
+          c8->Divide(1,3);
+          //c2->SetFillColor(42);
+          
+          
+          // data per hit
+          c8->cd(1);
+          hitsPhi->SetFillColor(5);
+          if (evNumber2>10)
+            hitsPhi->Fit("gaus");
+          hitsPhi->Draw();
+          
+           //data per track
+          c8->cd(2);
+          OriginalPhi->SetFillColor(5);
+          if (evNumber2>10)
+            OriginalPhi->Fit("gaus");
+          OriginalPhi->Draw();
+
+          //recontructed data
+          c8->cd(3);
+          Phi->SetFillColor(5);
+          if (evNumber2>10)
+            Phi->Fit("gaus");
+          Phi->Draw();
+
+          c9 = new TCanvas("c9","3D reconstruction of theta angle",75,75,300,1050);
+          c9->Divide(1,3);
+
+          // data per hit
+          c9->cd(1);
+          hitsTheta->SetFillColor(5);
+          if (evNumber2>10)
+            hitsTheta->Fit("gaus");
+          hitsTheta->Draw();
+          
+          //data per track
+          c9->cd(2);
+          OriginalTheta->SetFillColor(5);
+          if (evNumber2>10)
+            OriginalTheta->Fit("gaus");
+          OriginalTheta->Draw();
+
+          //recontructed data
+          c9->cd(3);
+          Theta->SetFillColor(5);
+          if (evNumber2>10)
+            Theta->Fit("gaus");
+          Theta->Draw();
+
+          c10 = new TCanvas("c10","3D reconstruction of cherenkov angle",100,100,300,1050);
+          c10->Divide(1,3);
+
+          // data per hit
+          c10->cd(1);
+          ckovangle->SetFillColor(5);
+          ckovangle->SetXTitle("angle (radians)");
+          if (evNumber2>10)
+            ckovangle->Fit("gaus");
+          ckovangle->Draw();
+          
+          //data per track
+          c10->cd(2);
+          OriginalOmega->SetFillColor(5);
+          OriginalOmega->SetXTitle("angle (radians)");
+          if (evNumber2>10)
+            OriginalOmega->Fit("gaus");
+          OriginalOmega->Draw();
+
+          //recontructed data
+          c10->cd(3);
+          Omega3D->SetFillColor(5);
+          Omega3D->SetXTitle("angle (radians)");
+          if (evNumber2>10)
+            Omega3D->Fit("gaus");
+          Omega3D->Draw(); 
+
+
+          c11 = new TCanvas("c11","3D reconstruction of mean radius",125,125,300,700);
+          c11->Divide(1,2);
+
+          // data per hit
+          c11->cd(1);
+          radius->SetFillColor(5);
+          radius->SetXTitle("radius (cm)");
+          radius->Draw();
+
+          //recontructed data
+          c11->cd(2);
+          MeanRadius->SetFillColor(5);
+          MeanRadius->SetXTitle("radius (cm)");
+          MeanRadius->Draw();
+
+          
+          c12 = new TCanvas("c12","Cerenkov angle vs. Momentum",150,150,550,350);
+
+          c12->cd(1);
+          identification->SetFillColor(5);
+          identification->SetXTitle("Momentum (GeV/c)");
+          identification->SetYTitle("Cherenkov angle (radians)");
+          
+          TF1 *pionplot = new TF1("pion","acos(sqrt((.139*.139+x*x)/(x*x*1.285*1.285)))",1,5);
+          TF1 *kaonplot = new TF1("kaon","acos(sqrt((.439*.439+x*x)/(x*x*1.285*1.285)))",1,5);
+          TF1 *protonplot = new TF1("proton","acos(sqrt((.938*.938+x*x)/(x*x*1.285*1.285)))",1,5);
+          
+          identification->Draw();
+
+          pionplot->SetLineColor(5);
+          pionplot->Draw("same");
+
+          kaonplot->SetLineColor(4);
+          kaonplot->Draw("same");
+
+          protonplot->SetLineColor(3);
+          protonplot->Draw("same");
+
+          c13 = new TCanvas("c13","Reconstruction Errors",200,200,900,350);
+          c13->Divide(3,1);
+
+          c13->cd(1);
+          PhiError->SetFillColor(5);
+          if (evNumber2>10)
+            PhiError->Fit("gaus");
+          PhiError->Draw();
+          c13->cd(2);
+          ThetaError->SetFillColor(5);
+          if (evNumber2>10)
+            ThetaError->Fit("gaus");
+          ThetaError->Draw();
+          c13->cd(3);
+          OmegaError->SetFillColor(5);
+          OmegaError->SetXTitle("angle (radians)");
+          if (evNumber2>10)
+            OmegaError->Fit("gaus");
+          OmegaError->Draw();
+          
+        }
+       
+       if(nrechits1D)
+        {
+          c9 = new TCanvas("c9","1D Reconstruction",100,100,1100,700);
+          c9->Divide(3,2);
+          //c5->SetFillColor(42);
+          
+          c9->cd(1);
+          ckovangle->SetFillColor(5);
+          ckovangle->SetXTitle("angle (radians)");
+          ckovangle->Draw();
+          
+          c9->cd(2);
+          radius->SetFillColor(5);
+          radius->SetXTitle("radius (cm)");
+          radius->Draw();
+          
+          c9->cd(3);
+          hc0->SetXTitle("pads");
+          hc0->Draw("box"); 
+          
+          c9->cd(5);
+          Omega1D->SetFillColor(5);
+          Omega1D->SetXTitle("angle (radians)");
+          Omega1D->Draw();
+          
+          c9->cd(4);
+          PhotonCer->SetFillColor(5);
+          PhotonCer->SetXTitle("angle (radians)");
+          PhotonCer->Draw();
+          
+          c9->cd(6);
+          PadsUsed->SetXTitle("pads");
+          PadsUsed->Draw("box"); 
+        }
+       
+       break;
+       
+     case 5:
+       
+       printf("Drawing histograms.../n");
+
+       c10 = new TCanvas("c10","Alice RICH digits",50,50,1200,700);
+       c1->Divide(4,2);
+       
+       c10->cd(1);
+       hc1->SetXTitle("ix (npads)");
+       hc1->Draw("box");
+       c10->cd(2);
+       hc2->SetXTitle("ix (npads)");
+       hc2->Draw("box");
+       c10->cd(3);
+       hc3->SetXTitle("ix (npads)");
+       hc3->Draw("box");
+       c10->cd(4);
+       hc4->SetXTitle("ix (npads)");
+       hc4->Draw("box");
+       c10->cd(5);
+       hc5->SetXTitle("ix (npads)");
+       hc5->Draw("box");
+       c10->cd(6);
+       hc6->SetXTitle("ix (npads)");
+       hc6->Draw("box");
+       c10->cd(7);
+       hc7->SetXTitle("ix (npads)");
+       hc7->Draw("box");
+       c10->cd(8);
+       hc0->SetXTitle("ix (npads)");
+       hc0->Draw("box");
+       c11 = new TCanvas("c11","Hits per type",100,100,600,700);
+       c11->Divide(2,2);
+       
+       c11->cd(1);
+       feedback->SetXTitle("x (cm)");
+       feedback->SetYTitle("y (cm)");
+       feedback->Draw();
+       
+       c11->cd(2);
+       mip->SetXTitle("x (cm)");
+       mip->SetYTitle("y (cm)");
+       mip->Draw();
+       
+       c11->cd(3);
+       cerenkov->SetXTitle("x (cm)");
+       cerenkov->SetYTitle("y (cm)"); 
+       cerenkov->Draw();
+       
+       c11->cd(4);
+       h->SetXTitle("x (cm)");
+       h->SetYTitle("y (cm)");
+       h->Draw();
+
+       c12 = new TCanvas("c12","Hits distribution",150,150,600,350);
+       c12->Divide(2,1);
+       
+       c12->cd(1);
+       hitsX->SetFillColor(5);
+       hitsX->SetXTitle("(cm)");
+       hitsX->Draw();
+       
+       c12->cd(2);
+       hitsY->SetFillColor(5);
+       hitsY->SetXTitle("(cm)");
+       hitsY->Draw();
+       
+       break;
+       
+     }
+       
+
+   printf("\nEnd of analysis\n");
+   printf("**********************************\n");
+}//void AliRICHv3::DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2)
+
+//__________________________________________________________________________________________________
+void AliRICHv3::MakeBranch(Option_t* option)
+{//Create Tree branches for the RICH.
+  if(GetDebug())Info("MakeBranch","Start with option= %s.",option);
+    
+  const Int_t kBufferSize = 4000;
+  char branchname[20];
+      
+   
+  const char *cH = strstr(option,"H");
+  const char *cD = strstr(option,"D");
+  const char *cR = strstr(option,"R");
+  const char *cS = strstr(option,"S");
+
+
+  if(cH&&TreeH()){
+    if(!fHits) fHits=new TClonesArray("AliRICHhit",1000  );
+    if(!fCerenkovs) fCerenkovs  = new TClonesArray("AliRICHCerenkov",1000);
+    MakeBranchInTree(TreeH(),"RICHCerenkov", &fCerenkovs, kBufferSize, 0) ;
+
+    if(!fSDigits) fSDigits    = new TClonesArray("AliRICHdigit",100000);
+    MakeBranchInTree(TreeH(),"RICHSDigits", &fSDigits, kBufferSize, 0) ;
+  }     
+  AliDetector::MakeBranch(option);//this is after cH because we need to guarantee that fHits array is created
+      
+  if(cS&&fLoader->TreeS()){  
+    if(!fSDigits) fSDigits=new TClonesArray("AliRICHdigit",100000);
+    MakeBranchInTree(fLoader->TreeS(),"RICH",&fSDigits,kBufferSize,0) ;
+  }
+   
+  int i;
+  if (cD&&fLoader->TreeD()){
+    if(!fDchambers){
+      fDchambers=new TObjArray(kNCH);    // one branch for digits per chamber
+      for(i=0;i<kNCH;i++){ 
+        fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i); 
+      }       
+    }
+    for (i=0; i<kNCH ;i++) 
+      {
+        sprintf(branchname,"%sDigits%d",GetName(),i+1);        
+        MakeBranchInTree(fLoader->TreeD(),branchname, &((*fDchambers)[i]), kBufferSize, 0);
+      }
+   }
+
+  if (cR&&gAlice->TreeR()){//one branch for raw clusters per chamber
+    Int_t i;
+    if (fRawClusters == 0x0 ) 
+     {
+       fRawClusters = new TObjArray(kNCH);
+       for (i=0; i<kNCH ;i++) 
+         {
+           fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i); 
+         }
+     }
+     
+    if (fRecHits1D == 0x0) 
+     {
+        fRecHits1D = new TObjArray(kNCH);
+        for (i=0; i<kNCH ;i++) 
+         {
+          fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
+         }
+     }
+
+    if (fRecHits3D == 0x0) 
+     {
+        fRecHits3D = new TObjArray(kNCH);
+        for (i=0; i<kNCH ;i++) 
+         {
+          fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
+         }
+     }
+       
+    for (i=0; i<kNCH ;i++){
+       sprintf(branchname,"%sRawClusters%d",GetName(),i+1);      
+       MakeBranchInTree(gAlice->TreeR(),branchname, &((*fRawClusters)[i]), kBufferSize, 0);
+       sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
+       MakeBranchInTree(fLoader->TreeR(),branchname, &((*fRecHits1D)[i]), kBufferSize, 0);
+       sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);  
+       MakeBranchInTree(fLoader->TreeR(),branchname, &((*fRecHits3D)[i]), kBufferSize, 0);
+     }
+   }//if (cR && gAlice->TreeR())
+  if(GetDebug())Info("MakeBranch","Stop.");   
+}
+//______________________________________________________________________________
+void AliRICHv3::SetTreeAddress()
+{//Set branch address for the Hits and Digits Tree.
+  if(GetDebug())Info("SetTreeAddress","Start.");
+  
+  char branchname[20];
+  Int_t i;
+
+    
+  TBranch *branch;
+  TTree *treeH = fLoader->TreeH();
+  TTree *treeD = fLoader->TreeD();
+  TTree *treeR = fLoader->TreeR();
+  TTree *treeS = fLoader->TreeS();
+    
+  if(treeH){
+    if(GetDebug())Info("SetTreeAddress","tree H is requested.");
+    if(fHits==0x0) fHits=new TClonesArray("AliRICHhit",1000); 
+    
+    branch = treeH->GetBranch("RICHCerenkov");
+    if(branch){
+      if (fCerenkovs == 0x0) fCerenkovs  = new TClonesArray("AliRICHCerenkov",1000); 
+        branch->SetAddress(&fCerenkovs);
+    }
+       
+      branch = treeH->GetBranch("RICHSDigits");
+      if (branch) 
+       {
+         if (fSDigits == 0x0) fSDigits    = new TClonesArray("AliRICHdigit",100000);
+         branch->SetAddress(&fSDigits);
+       }
+  }//if(treeH)
+   //this is after TreeH because we need to guarantee that fHits array is created
+  AliDetector::SetTreeAddress();
+    
+  if(treeS){
+    if(GetDebug())Info("SetTreeAddress","tree S is requested.");
+    branch = treeS->GetBranch("RICH");
+    if(branch){
+      if(!fSDigits) fSDigits=new TClonesArray("AliRICHdigit",100000);
+      branch->SetAddress(&fSDigits);
+    }
+  }
+    
+    
+  if(treeD){
+    if(GetDebug())Info("SetTreeAddress","tree D is requested.");
+
+      if (fDchambers == 0x0) 
+        {
+           fDchambers = new TObjArray(kNCH);
+           for (i=0; i<kNCH ;i++) 
+             {
+               fDchambers->AddAt(new TClonesArray("AliRICHDigit",10000), i); 
+             }
+        }
+      
+      for (i=0; i<kNCH; i++) {
+        sprintf(branchname,"%sDigits%d",GetName(),i+1);
+        if (fDchambers) {
+           branch = treeD->GetBranch(branchname);
+           if (branch) branch->SetAddress(&((*fDchambers)[i]));
+        }
+      }
+    }
+    
+  if(treeR){
+    if(GetDebug())Info("SetTreeAddress","tree R is requested.");
+
+    if (fRawClusters == 0x0 ) 
+     {
+       fRawClusters = new TObjArray(kNCH);
+       for (i=0; i<kNCH ;i++) 
+         {
+           fRawClusters->AddAt(new TClonesArray("AliRICHRawCluster",10000), i); 
+         }
+     }
+     
+    if (fRecHits1D == 0x0) 
+     {
+        fRecHits1D = new TObjArray(kNCH);
+        for (i=0; i<kNCH ;i++) 
+         {
+          fRecHits1D->AddAt(new TClonesArray("AliRICHRecHit1D",1000), i);
+         }
+     }
+
+    if (fRecHits3D == 0x0) 
+     {
+        fRecHits3D = new TObjArray(kNCH);
+        for (i=0; i<kNCH ;i++) 
+         {
+          fRecHits3D->AddAt(new TClonesArray("AliRICHRecHit3D",1000), i);
+         }
+     }
+    
+    for (i=0; i<kNCH; i++) {
+         sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
+         if (fRawClusters) {
+             branch = treeR->GetBranch(branchname);
+             if (branch) branch->SetAddress(&((*fRawClusters)[i]));
+         }
+    }
+      
+    for (i=0; i<kNCH; i++) {
+       sprintf(branchname,"%sRecHits1D%d",GetName(),i+1);
+       if (fRecHits1D) {
+         branch = treeR->GetBranch(branchname);
+         if (branch) branch->SetAddress(&((*fRecHits1D)[i]));
+         }
+     }
+      
+     for (i=0; i<kNCH; i++) {
+       sprintf(branchname,"%sRecHits3D%d",GetName(),i+1);
+       if (fRecHits3D) {
+         branch = treeR->GetBranch(branchname);
+         if (branch) branch->SetAddress(&((*fRecHits3D)[i]));
+         }
+      } 
+      
+  }//if(treeR)
+  if(GetDebug())Info("SetTreeAddress","Stop.");
+}//void AliRICHv3::SetTreeAddress()
index 0ffa15e..673bd8e 100644 (file)
@@ -6,6 +6,8 @@
 
 #include "AliRICH.h"
 
+class AliRICHSDigit;
+
 class AliRICHv3 : public AliRICH 
 {    
 public:
@@ -18,6 +20,18 @@ public:
   virtual void   CreateGeometry();  
   virtual void   BuildGeometry();   
   virtual void   Init();            // Makes nothing for a while          
+          void   DiagnosticsFE(Int_t evNumber1,Int_t evNumber2);
+          void   DiagnosticsSE(Int_t diaglevel,Int_t evNumber1,Int_t evNumber2);
+         Int_t   Hits2SDigits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol, ResponseType res);
+  void     SetGeometryModel(Int_t c,AliRICHGeometry *pRICHGeo)                    {C(c)->SetGeometryModel(pRICHGeo);}
+  void     SetSegmentationModel(Int_t c, AliSegmentation *pAliSeg)                {C(c)->SetSegmentationModel(pAliSeg);}
+  void     SetResponseModel(Int_t c, AliRICHResponse *pRICHRes)                   {C(c)->SetResponseModel(pRICHRes);}
+  void     SetReconstructionModel(Int_t c, AliRICHClusterFinder *pRICHReco)       {C(c)->SetReconstructionModel(pRICHReco);}
+  AliRICHGeometry* GetGeometryModel(Int_t c=1)                               const{return C(c)->GetGeometryModel();}    
+  AliSegmentation* GetSegmentationModel(Int_t c=1)                           const{return C(c)->GetSegmentationModel();}
+  AliRICHResponse* GetResponseModel(Int_t c=1)                               const{return C(c)->GetResponseModel();}
+  void MakeBranch(Option_t* option);//virtual
+  void SetTreeAddress();//virtual
 private:
   Double_t* RotateXY(const Double_t* r, Double_t a);//Rotation in the X-Y plane in G3 notation
   ClassDef(AliRICHv3,1)  //RICH full version, configurable with azimuthal rotation