New Geom, support for ESD- major change
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 25 Sep 2004 09:48:53 +0000 (09:48 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 25 Sep 2004 09:48:53 +0000 (09:48 +0000)
29 files changed:
RICH/AliGenRadioactive.cxx [new file with mode: 0644]
RICH/AliGenRadioactive.h [new file with mode: 0644]
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHChamber.cxx
RICH/AliRICHChamber.h
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHDisplFast.cxx
RICH/AliRICHHelix.cxx [new file with mode: 0644]
RICH/AliRICHHelix.h [new file with mode: 0644]
RICH/AliRICHParam.cxx
RICH/AliRICHParam.h
RICH/AliRICHRecon.cxx
RICH/AliRICHRecon.h
RICH/AliRICHReconstructor.cxx
RICH/AliRICHReconstructor.h
RICH/AliRICHTracker.cxx [new file with mode: 0644]
RICH/AliRICHTracker.h [new file with mode: 0644]
RICH/AliRICHv0.cxx
RICH/AliRICHv0.h
RICH/AliRICHv1.cxx
RICH/CreateConfig.C
RICH/Geom.C
RICH/Opticals.h
RICH/RICHLinkDef.h
RICH/RichBatch.C
RICH/api.txt
RICH/libRICH.pkg
RICH/menu.C

diff --git a/RICH/AliGenRadioactive.cxx b/RICH/AliGenRadioactive.cxx
new file mode 100644 (file)
index 0000000..1244b35
--- /dev/null
@@ -0,0 +1,45 @@
+#include "AliGenRadioactive.h"
+#include <TPDGCode.h>
+#include <AliRun.h>
+#include <TDatabasePDG.h>
+#include <AliLog.h>
+
+//__________________________________________________________________________________________________
+AliGenRadioactive::AliGenRadioactive(Int_t iSrcNucleus,Int_t iNsecondaries):AliGenerator()
+{
+//Main ctor. Used to define radioactive source.
+  Double_t e[100],a[100];//arrays to store experimental points
+  Int_t nPoints; 
+  switch(iSrcNucleus){
+    case kSr90: fPartId=kElectron; nPoints=46;  //experimental part
+  a[ 0]=0.08605; a[ 1]=0.0878;  a[ 2]=0.08705; a[ 3]=0.07855; a[ 4]=0.0709; a[ 5]=0.0647;  a[ 6]=0.05015; a[ 7]=0.0372; a[ 8]=0.0268;  a[ 9]=0.0215;
+  a[10]=0.0157;  a[11]=0.01685; a[12]=0.01745; a[13]=0.01645; a[14]=0.0175; a[15]=0.01635; a[16]=0.01825; a[17]=0.0177; a[18]=0.01735; a[19]=0.0161;
+  a[20]=0.0159;  a[21]=0.0176;  a[22]=0.01605; a[23]=0.0161;  a[24]=0.01495;a[25]=0.01595; a[26]=0.01525; a[27]=0.0138; a[28]=0.0121;  a[29]=0.0101;
+  a[30]=0.01175; a[31]=0.01095; a[32]=0.0089;  a[33]=0.0091;  a[34]=0.0625; a[35]=0.0505;  a[36]=0.0475;  a[37]=0.0039; a[38]=0.0031;  a[39]=0.0028;
+  a[40]=0.0025;  a[41]=0.0017;  a[42]=4.5e-4;  a[43]=4.5e-4;  a[44]=1.5e-4; a[45]=0;     
+  break;
+    default:    AliError("Wrong source nucleus specified");  return;   
+  }
+  for(Int_t i=0;i<nPoints;i++) e[i]=0.001*(i*0.05+0.025); //kinetic energy GeV
+  fGenH1=new TH1F("Sr90","Sr90 generator hist",nPoints-1,0,e[nPoints-1]);  
+  for(Int_t i=0;i<nPoints;i++) fGenH1->Fill(e[i],a[i]);
+  fNpart=iNsecondaries;
+}
+//__________________________________________________________________________________________________
+void AliGenRadioactive::Generate()
+{
+// Generate one trigger
+  Int_t nt=0;
+  Double_t ekin=0,p=0,theta=0,phi=0,x=0,y=0,z=0,px=0,py=0,pz=0,polx=0,poly=0,polz=0;
+  Double_t m=gAlice->PDGDB()->GetParticle(fPartId)->Mass();
+  for(Int_t i=0;i<fNpart;i++){
+    x=fOrigin.At(0)+fOsigma.At(0)*(Rndm()-0.5); y=fOrigin.At(1)+fOsigma.At(1)*(Rndm()-0.5); z=fOrigin.At(2)+fOsigma.At(2)*(Rndm()-0.5);
+    ekin=fGenH1->GetRandom();  p=TMath::Sqrt(ekin*(2*m+ekin)); 
+    theta=Rndm()*fThetaMax; phi=Rndm()*fPhiMax;
+    px=p*TMath::Cos(theta)*TMath::Cos(phi);  py=p*TMath::Cos(theta)*TMath::Sin(phi);  pz=p*TMath::Sin(theta);
+//   AliDebug(1,Form("Origin=(%5.2f,%5.2f,%5.2f) Ekin=%5.3fMeV,P=%5.3fMeV (%5.2f,%5.2f,%5.2f)",x,y,z,1000*ekin,1000*p,1000*px,1000*py,1000*pz));
+    PushTrack(fTrackIt,-1,fPartId,px,py,pz,ekin+m,
+                                   x, y, z,0,
+                                    polx=0,poly=0,polz=0,kPPrimary,nt);//cm, GeV
+  }
+}
diff --git a/RICH/AliGenRadioactive.h b/RICH/AliGenRadioactive.h
new file mode 100644 (file)
index 0000000..0a257f4
--- /dev/null
@@ -0,0 +1,28 @@
+#ifndef AliGenRadioactive_h
+#define AliGenRadioactive_h
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+// Simple radioactive source generator.
+// Currently Sr90 only
+
+#include <AliGenerator.h>
+#include <TH1.h>
+
+typedef enum {kSr90} RadioNuclei_t;     
+
+class AliGenRadioactive : public AliGenerator
+{
+ public:
+           AliGenRadioactive() : AliGenerator(),fPartId(0),fGenH1(0) {;}//default ctor
+           AliGenRadioactive(Int_t iSource,Int_t iNparts);                  //main ctor
+  virtual ~AliGenRadioactive()                                           {if(fGenH1) delete fGenH1;}//dtor
+  virtual void Generate();                                                              //interface from AliGenerator, generates current event 
+
+
+protected:
+  Int_t    fPartId;          //ID of secondary from radioactive nucleus
+  TH1F    *fGenH1;           //Histogram containg exp spectrum to be used to sample the energy
+  ClassDef(AliGenRadioactive,1) // Radioactive source generator
+};
+#endif
index 72e34ed..842bd21 100644 (file)
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TStopwatch.h>
+#include <AliLog.h>
  
 ClassImp(AliRICHhit)
 //__________________________________________________________________________________________________
 void AliRICHhit::Print(Option_t*)const
 {
-  ::Info("hit","Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)"
-      ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z());
+  AliInfo(Form("Ch=%1i, TID=%6i, eloss=%9.3f eV, in-out dist=%9.4f, OUT(%7.2f,%7.2f,%7.2f)"
+      ,fChamber,fTrack,fEloss*1e9,Length(),fOutX3.X(),fOutX3.Y(),fOutX3.Z()));
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICHdigit)
 //__________________________________________________________________________________________________
 void AliRICHdigit::Print(Option_t*)const
 {
-  ::Info("digit","cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i",
-                  fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]);
+  AliInfo(Form("cfm=%9i, cs=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i",
+                  fCFM,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]));
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICHcluster)
@@ -66,8 +67,8 @@ void AliRICHcluster::Print(Option_t*)const
     ::Info("cluster","cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits",
                              fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,fDigits->GetEntriesFast());
   else
-    ::Info("cluster","cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits",
-                             fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,0);
+    AliInfo(Form("cfm=%10i, cs=%2i, SiMa=%6i, Shape=%5i, x=%7.3f, y=%7.3f, Q=%6i, %s with %i digits",
+                             fCFM,fChamber,fSize,fShape,fX,fY,fQdc,status,0));
     
 }
 //__________________________________________________________________________________________________
@@ -93,17 +94,17 @@ AliRICH::AliRICH(const char *name, const char *title)
         :AliDetector(name,title),fpParam(new AliRICHParam),fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0)
 {
 //Named ctor
-  if(GetDebug())Info("named ctor","Start.");
+  AliDebug(1,"Start.");
 //AliDetector ctor deals with Hits and Digits (reset them to 0, does not create them)
   CreateHits();          gAlice->GetMCApp()->AddHitList(fHits);
   fCounters.ResizeTo(20); fCounters.Zero();
-  if(GetDebug())Info("named ctor","Stop.");
+  AliDebug(1,"Stop.");
 }//AliRICH::AliRICH(const char *name, const char *title)
 //__________________________________________________________________________________________________
 AliRICH::~AliRICH()
 {
 //dtor
-  if(GetDebug()) Info("dtor","Start.");
+  AliDebug(1,"Start.");
 
   if(fpParam)    delete fpParam;
   
@@ -112,13 +113,13 @@ AliRICH::~AliRICH()
   if(fDigits)    delete fDigits;
   if(fDigitsNew) {fDigitsNew->Delete();   delete fDigitsNew;}
   if(fClusters)  {fClusters->Delete();    delete fClusters;}
-  if(GetDebug()) Info("dtor","Stop.");    
+  AliDebug(1,"Stop.");    
 }//AliRICH::~AliRICH()
 //__________________________________________________________________________________________________
 void AliRICH::Hits2SDigits()
 {
 // Create a list of sdigits corresponding to list of hits. Every hit generates one or more sdigits.
-  if(GetDebug()) Info("Hit2SDigits","Start.");
+  AliDebug(1,"Start.");
   for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
     GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
   
@@ -130,16 +131,16 @@ void AliRICH::Hits2SDigits()
       GetLoader()->TreeH()->GetEntry(iPrimN);
       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop 
         AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit                
-        TVector2 x2 = C(pHit->C())->Glob2Loc(pHit->OutX3());//hit position in the chamber local system
+        TVector2 x2 = C(pHit->C())->Mrs2Pc(pHit->OutX3());//hit position in photocathode plane
         Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
         if(iTotQdc==0) continue;
         TVector area=P()->Loc2Area(x2);//determine affected pads, dead zones analysed inside
-        if(GetDebug()) Info("Hits2SDigits","hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc);
+        AliDebug(1,Form("hit(%6.2f,%6.2f)->area(%3.0f,%3.0f)-(%3.0f,%3.0f) QDC=%4i",x2.X(),x2.Y(),area[0],area[1],area[2],area[3],iTotQdc));
         TVector pad(2);
         for(pad[1]=area[1];pad[1]<=area[3];pad[1]++)//affected pads loop
           for(pad[0]=area[0];pad[0]<=area[2];pad[0]++){                    
             Double_t padQdc=iTotQdc*P()->FracQdc(x2,pad);
-            if(GetDebug()) Info("Hits2SDigits","current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc);
+            AliDebug(1,Form("current pad(%3.0f,%3.0f) with QDC  =%6.2f",pad[0],pad[1],padQdc));
             if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
           }//affected pads loop 
       }//hits loop
@@ -150,127 +151,135 @@ void AliRICH::Hits2SDigits()
   }//events loop  
   GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
   GetLoader()->UnloadSDigits();  
-  if(GetDebug()) Info("Hit2SDigits","Stop.");
+  AliDebug(1,"Stop.");
 }//Hits2SDigits()
 //__________________________________________________________________________________________________
 void AliRICH::BuildGeometry() 
 {
 //Builds a TNode geometry for event display
-  if(GetDebug())Info("BuildGeometry","Start.");
+  AliInfo("Start.");
   
   TNode *node, *subnode, *top;
   top=gAlice->GetGeometry()->GetNode("alice");
-  
-  new TBRIK("S_RICH","S_RICH","void",71.09999,11.5,73.15);
 
-  Float_t wid=P()->SectorSizeX();
-  Float_t len=P()->SectorSizeY();
-  new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2);
-  
+  Float_t widx=P()->SectorSizeX();
+  Float_t leny=P()->SectorSizeY();
+  Float_t dz = P()->Zfreon()+P()->Zwin()+P()->Pc2Win();
+  Float_t dead = P()->DeadZone();
+
+  new TBRIK("RICH","RICH","void",widx+dead/2,leny+leny/2+dead,dz+0.1); //RICH chamber
+  new TBRIK("RPC" ,"RPC" ,"void",widx/2,leny/2,0.01);                  //RICH sector 
+
   for(int i=1;i<=kNchambers;i++){
     top->cd();
-    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 = new TNode(Form("RICH%i",i),Form("RICH%i",i),"RICH",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),C(i)->RotMatrixName());
     node->SetLineColor(kRed);
     node->cd();
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,len/2+P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2,-leny-dead/2,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,len/2+P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2,-leny-dead/2,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,len/2+P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2,           0,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",wid+P()->DeadZone(),5,-len/2-P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2,           0,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",0,5,-len/2 -P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC",-widx/2-dead/2, leny+dead/2,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
-    subnode = new TNode("PHOTO1","PHOTO1","PHOTO",-wid-P()->DeadZone(),5,-len/2 - P()->DeadZone()/2,"");
+    subnode = new TNode("PHOTO1","PHOTO1","RPC", widx/2+dead/2, leny+dead/2,dz,"");
     subnode->SetLineColor(kGreen);
     fNodes->Add(subnode);
     fNodes->Add(node);
-  }  
-  if(GetDebug())Info("BuildGeometry","Stop.");    
+  }
+
+  AliDebug(1,"Stop.");    
 }//void AliRICH::BuildGeometry()
 
 //______________________________________________________________________________
 void AliRICH::CreateMaterials()
 {
 // Definition of available RICH materials  
-#include "Opticals.h"
         
-  Float_t a=0,z=0,den=0,radl=0,absl=0;
+  Int_t   material=0; //tmp material id number
+  Float_t a=0,z=0,den=0,radl=0,absl=0; //tmp material parameters
   Float_t tmaxfd=-10.0, deemax=-0.2, stemax=-0.1,epsil=0.001, stmin=-0.001; 
   Int_t   isxfld = gAlice->Field()->Integ();
   Float_t sxmgmx = gAlice->Field()->Max();
-  Int_t material;
     
-  AliMaterial(material=1, "Air     $",a=14.61,z=7.3, den=0.001205,radl=30420.0,absl=67500);//(Air)
-  AliMedium(1, "DEFAULT MEDIUM AIR$",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  Float_t aAir[4]={12.,14.,16.,36.};  Float_t zAir[4]={6.,7.,8.,18.}; Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};//total 0.9999999
+  AliMixture(++material, "RichAir",aAir,zAir,den=0.00120479,4,wAir);                                          //1 (Air) 0.01% C 75% N  23% O 1% Ar
+  AliMedium(kAir, "RichAir",material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMaterial( 6, "HON",      a=12.01,z=6.0, den=0.1,     radl=18.8,   absl=0);    //(C)-equivalent radl
-  AliMedium(2, "HONEYCOMB$", 6, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMaterial(++material, "RichRohacell", a=12.01,z=6.0, den=0.1,     radl=18.8,   absl=0);                   //2 Rohacell 51 C-equiv radl rad cover
+  AliMedium(kRoha, "RichRohacell", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMaterial(16, "CSI",      a=12.01,z=6.0, den=0.1,     radl=18.8,   absl=0);    //CsI-radl equivalent
-  AliMedium(kCSI, "CSI$", 16, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
+  Float_t  aQuartz[2]={28.09,16.0};  Float_t  zQuartz[2]={14.00, 8.0};  Float_t  wQuartz[2]={1,2};
+  AliMixture(++material, "RichSiO2",aQuartz,zQuartz,den=2.64,-2, wQuartz);                                    //3 Quarz (SiO2) -trasparent rad window
+  AliMedium(kSiO2, "RichSiO2",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMaterial(11, "GRI",      a=63.54,z=29.0,den=8.96,    radl=1.43,   absl=0);    //anode grid (Cu) 
-  AliMedium(7, "GRID$", 11, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  Float_t  aFreon[2]={12,19};  Float_t  zFreon[2]={6,9};  Float_t wmatFreon[2]={6,14};
+  AliMixture(++material, "RichC6F14",aFreon,zFreon,den=1.68,-2,wmatFreon);                                    //4 Freon (C6F14) 
+  AliMedium(kC6F14, "RichC6F14",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  
+  Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wMethane[2]={1,4};
+  AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane);                        //5,9 methane (CH4) normal and for Gap    
+  AliMedium(kCH4, "RichCH4"   , material, 1, isxfld, sxmgmx, tmaxfd, stemax,  deemax, epsil,  stmin);  
+  AliMixture (++material, "RichCH4", aMethane, zMethane, den=7.17e-4,-2, wMethane);                        //5,9 methane (CH4) normal and for Gap    
+  AliMedium(kGap, "RichCH4gap$", material, 1, isxfld, sxmgmx, tmaxfd, 0.1   , -deemax, epsil, -stmin);
+    
+  AliMaterial(++material, "RichCsI",      a=12.01,z=6.0, den=0.1,     radl=18.8,   absl=0);                   //6 CsI-radl equivalent
+  AliMedium(kCsI, "RichCsI$", material, 1, isxfld, sxmgmx,tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMaterial(50, "ALUM",     a=26.98,z=13.0,den=2.699,     radl=8.9,    absl=0);    //aluminium sheet (Al)
-  AliMedium(10, "ALUMINUM$", 50, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMaterial(++material, "GridCu",    a=63.54,z=29.0,den=8.96,    radl=1.43,   absl=0);                   //7 anode grid (Cu) 
+  AliMedium(kGridCu, "GRID$", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMaterial(material=31, "COPPER$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);    //(Cu)
-  AliMedium(12, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMixture (++material, "OpSiO2",aQuartz, zQuartz, den=2.64, -2, wQuartz);                             //8 Quarz (SiO2) - opaque
+  AliMedium(kOpSiO2, "QUARTZO$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  Float_t  aQuartz[2]={28.09,16.0};  Float_t  zQuartz[2]={14.00, 8.0};  Float_t  wmatQuartz[2]={1,2};
-  AliMixture (20, "QUA",aQuartz,zQuartz,den=2.64,-2, wmatQuartz);//Quarz (SiO2) - trasnparent 
-  AliMedium(3, "QUARTZ$", 20, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMaterial(++material, "ALUM",     a=26.98,z=13.0,den=2.699,     radl=8.9,    absl=0);                 //10 aluminium sheet (Al)
+  AliMedium(kAl, "ALUMINUM$", material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMixture (21, "QUAO",aQuartz, zQuartz, den=2.64, -2, wmatQuartz);//Quarz (SiO2) - opaque
-  AliMedium(8, "QUARTZO$", 21, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  Float_t aGlass[5]={12.01,28.09,16,10.8,23}; Float_t zGlass[5]={6,14,8,5,11};  Float_t wGlass[5]={0.5,0.105,0.355,0.03,0.01};
+  AliMixture(++material,"GLASS",aGlass, zGlass, den=1.74, 5, wGlass);                                    //11 Glass 50%-C 10.5%-Si 35.5%-O 3%-B 1%-Na
+  AliMedium(kGlass, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  Float_t  aFreon[2]={12,19};  Float_t  zFreon[2]={6,9};  Float_t wmatFreon[2]={6,14};
-  AliMixture (material=30, "C6F14",aFreon,zFreon,den=1.68,-2,wmatFreon);//Freon (C6F14) 
-  AliMedium(4, "C6F14$",material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  AliMaterial(++material, "COPPER$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);                   //12 Cu
+  AliMedium(kCu, "PCB_COPPER", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  
+  
+  AliMaterial(++material, "W$",  a=183.84,z=74.0,den=19.3,    radl=0.35,    absl=185.0/den);              //13 W - anod wires
+  AliMedium(kW, "W", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wmatMethane[2]={1,4};
-  AliMixture (material=40, "CH4", aMethane, zMethane, den=7.17e-4,-2, wmatMethane);//methane (CH4)     
-  AliMedium(kCH4, "CH4$"   , material, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);  
-  AliMedium(kGAP, "CH4GAP$", material, 1, isxfld, sxmgmx,tmaxfd, 0.1, -deemax, epsil, -stmin);
   
   if(P()->IsRadioSrc()){
-    AliMaterial(material=45, "STEEL$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);    //Steel
+    AliInfo("Special radioactive source materials");
+    AliMaterial(++material, "STEEL$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);                  //14 Steel
     AliMedium(kSteel, "STEEL", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-    AliMaterial(material=46, "PERPEX$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);    //Perpex
+    AliMaterial(++material, "PERPEX$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);                 //15 Perpex
     AliMedium(kPerpex, "PERPEX", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
     
-    AliMaterial(material=47, "STRONZIUM$",  a=87.62,z=38.0,den=2.54,    radl=4.24,    absl=0); //Sr90
+    AliMaterial(++material, "STRONZIUM$",  a=87.62,z=38.0,den=2.54,    radl=4.24,    absl=0);             //16 Sr90
     AliMedium(kSr90, "STRONZIUM", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   }
   
-  Float_t aGlass[5]={12.01, 28.09, 16.,   10.8,  23.};
-  Float_t zGlass[5]={ 6.,   14.,    8.,    5.,   11.};
-  Float_t wGlass[5]={ 0.5,  0.105, 0.355, 0.03,  0.01};
-  AliMixture(material=32, "GLASS",aGlass, zGlass, den=1.74, 5, wGlass);//Glass 50%C+10.5%Si+35.5%O+3% + 1%
-  AliMedium(11, "GLASS", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
-            
-  Int_t *idtmed = fIdtmed->GetArray()-999;
-  gMC->SetCerenkov(idtmed[1000], kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);
-  gMC->SetCerenkov(idtmed[1001], kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);
-  gMC->SetCerenkov(idtmed[1002], kNbins, aPckov, aAbsSiO2,   aQeAll, aIdxSiO2);
-  gMC->SetCerenkov(idtmed[1003], kNbins, aPckov, aAbsC6F14,  aQeAll, aIdxC6F14);
-  gMC->SetCerenkov(idtmed[1004], kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);
-  gMC->SetCerenkov(idtmed[1005], kNbins, aPckov, aAbsCsI,    aQeCsI, aIdxCH4);
-  gMC->SetCerenkov(idtmed[1006], kNbins, aPckov, aAbsGrid,   aQeAll, aIdxGrid);
-  gMC->SetCerenkov(idtmed[1007], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);
-  gMC->SetCerenkov(idtmed[1008], kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);
-  gMC->SetCerenkov(idtmed[1009], kNbins, aPckov, aAbsGrid,   aQeAll, aIdxGrid);
-  gMC->SetCerenkov(idtmed[1010], kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);
-    
+//Optical properties:
+#include "Opticals.h"
+  gMC->SetCerenkov((*fIdtmed)[kAir]      , kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);       //1 Air
+  gMC->SetCerenkov((*fIdtmed)[kRoha]     , kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);       //2 Honeycomb  
+  gMC->SetCerenkov((*fIdtmed)[kSiO2]     , kNbins, aPckov, aAbsSiO2,   aQeAll, aIdxSiO2);      //3 Quartz SiO2 
+  gMC->SetCerenkov((*fIdtmed)[kC6F14]    , kNbins, aPckov, aAbsC6F14,  aQeAll, aIdxC6F14);     //4 Freon C6F14
+  gMC->SetCerenkov((*fIdtmed)[kCH4]      , kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);       //5 Methane CH4 
+  gMC->SetCerenkov((*fIdtmed)[kCsI]      , kNbins, aPckov, aAbsCsI,    aQeCsI, aIdxCH4);       //6 CsI
+  gMC->SetCerenkov((*fIdtmed)[kGridCu]   , kNbins, aPckov, aAbsGrid,   aQeAll, aIdxGrid);      //7 grid Cu
+  gMC->SetCerenkov((*fIdtmed)[kOpSiO2]   , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);    //8 Opaque quartz SiO2
+  gMC->SetCerenkov((*fIdtmed)[kGap]      , kNbins, aPckov, aAbsCH4,    aQeAll, aIdxCH4);       //9 Special methane gap
+  gMC->SetCerenkov((*fIdtmed)[kAl]       , kNbins, aPckov, aAbsGrid,   aQeAll, aIdxGrid);      //10 Aluminium
+  gMC->SetCerenkov((*fIdtmed)[kGlass]    , kNbins, aPckov, aAbsOpSiO2, aQeAll, aIdxOpSiO2);    //11 Glass    
 }//void AliRICH::CreateMaterials()
 //__________________________________________________________________________________________________
 Float_t AliRICH::Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)const
@@ -365,7 +374,7 @@ Float_t AliRICH::AbsoCH4(Float_t x)const
 void AliRICH::MakeBranch(Option_t* option)
 {
 //Create Tree branches for the RICH.
-  if(GetDebug())Info("MakeBranch","Start with option= %s.",option);
+  AliDebug(1,Form("Start with option= %s.",option));
     
   const Int_t kBufferSize = 4000;
       
@@ -395,29 +404,29 @@ void AliRICH::MakeBranch(Option_t* option)
     for(Int_t i=0;i<kNchambers;i++)
       MakeBranchInTree(fLoader->TreeR(),Form("%sClusters%d",GetName(),i+1), &((*fClusters)[i]), kBufferSize, 0);    
   }//R
-  if(GetDebug())Info("MakeBranch","Stop.");   
+  AliDebug(1,"Stop.");   
 }//void AliRICH::MakeBranch(Option_t* option)
 //__________________________________________________________________________________________________
 void AliRICH::SetTreeAddress()
 {
 //Set branch address for the Hits and Digits Tree.
-  if(GetDebug())Info("SetTreeAddress","Start.");
+  AliDebug(1,"Start.");
       
   TBranch *branch;
     
   if(fLoader->TreeH()){//H
-    if(GetDebug())Info("SetTreeAddress","tree H is requested.");
+    AliDebug(1,"tree H is requested.");
     CreateHits();//branch map will be in AliDetector::SetTreeAddress    
   }//H
   AliDetector::SetTreeAddress();//this is after TreeH because we need to guarantee that fHits array is created
 
   if(fLoader->TreeS()){//S
-    if(GetDebug())Info("SetTreeAddress","tree S is requested.");
+    AliDebug(1,"tree S is requested.");
     branch=fLoader->TreeS()->GetBranch(GetName());        if(branch){CreateSDigits();   branch->SetAddress(&fSdigits);}
   }//S
     
   if(fLoader->TreeD()){//D    
-    if(GetDebug())Info("SetTreeAddress","tree D is requested.");
+    AliDebug(1,"tree D is requested.");
     for(int i=0;i<kNchambers;i++){      
       branch=fLoader->TreeD()->GetBranch(Form("%s%d",GetName(),i+1)); 
       if(branch){CreateDigits(); branch->SetAddress(&((*fDigitsNew)[i]));}
@@ -425,13 +434,13 @@ void AliRICH::SetTreeAddress()
   }//D
     
   if(fLoader->TreeR()){//R
-    if(GetDebug())Info("SetTreeAddress","tree R is requested.");
+    AliDebug(1,"tree R is requested.");
     for(int i=0;i<kNchambers;i++){         
       branch=fLoader->TreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1));
       if(branch){CreateClusters(); branch->SetAddress(&((*fClusters)[i]));}
     }
   }//R
-  if(GetDebug())Info("SetTreeAddress","Stop.");
+  AliDebug(1,"Stop.");
 }//void AliRICH::SetTreeAddress()
 //__________________________________________________________________________________________________
 void AliRICH::Print(Option_t *option)const
@@ -445,216 +454,100 @@ void AliRICH::Print(Option_t *option)const
 void AliRICH::CreateGeometry()
 {
 //Creates detailed geometry simulation (currently GEANT volumes tree)         
-  if(GetDebug())Info("CreateGeometry","Start main.");
-//Opaque quartz thickness
-  Float_t oquaThickness = .5;
-//CsI dimensions
-  Float_t pcX=P()->PcSizeX();
-  Float_t pcY=P()->PcSizeY();
-  
-  Int_t *idtmed = fIdtmed->GetArray()-999;
-    
-  Int_t i;
-  Float_t zs;
-  Int_t idrotm[1099];
+  AliDebug(1,"Start main.");
+  Double_t cm=1,mm=0.1*cm,mkm=0.001*cm;
   Float_t par[3];
-  par[0]=68.8;par[1]=13   ;par[2]=70.86;  gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kAl], par,3);//External aluminium box 
-  par[0]=66.3;par[1]=13   ;par[2]=68.35;  gMC->Gsvolu("SRIC","BOX ",(*fIdtmed)[kAir],par,3);//Air 
-  par[0]=66.3;par[1]=0.025;par[2]=68.35;  gMC->Gsvolu("ALUM","BOX ",(*fIdtmed)[kAl], par,3);//Aluminium sheet 
-//Air 2 (cutting the lower part of the box)
-//  par[0]=1.25;    par[1] = 3;    par[2] = 70.86;   gMC->Gsvolu("AIR2", "BOX ", idtmed[1000], par, 3);
-//Air 3 (cutting the lower part of the box)
-//  par[0]=66.3;    par[1] = 3;  par[2] = 1.2505;    gMC->Gsvolu("AIR3", "BOX ", idtmed[1000], par, 3);
-//Honeycomb 
-  par[0]=66.3;par[1]=0.188;  par[2] = 68.35;       gMC->Gsvolu("HONE", "BOX ", idtmed[1001], par, 3);
-
-  //par[0] = 66.5; par[1] = .025; par[2] = 63.1;
-//Quartz 
-  par[0]=P()->QuartzWidth()/2;par[1]=P()->QuartzThickness()/2;par[2]=P()->QuartzLength()/2;
-  gMC->Gsvolu("QUAR", "BOX ", idtmed[1002], par, 3);
-//Spacers (cylinders) 
-  par[0]=0.;par[1]=.5;par[2]=P()->FreonThickness()/2;  gMC->Gsvolu("SPAC", "TUBE", idtmed[1002], par, 3);    
-//Feet (freon slabs supports)
-  par[0] = .7;  par[1] = .3;  par[2] = 1.9;        gMC->Gsvolu("FOOT", "BOX", idtmed[1009], par, 3);
-//Opaque quartz 
-  par[0]=P()->QuartzWidth()/2;par[1]= .2;par[2]=P()->QuartzLength()/2;
-  gMC->Gsvolu("OQUA", "BOX ", idtmed[1007], par, 3);
-//Frame of opaque quartz
-  par[0]=P()->OuterFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->OuterFreonLength()/2; 
-  gMC->Gsvolu("OQF1", "BOX ", idtmed[1007], par, 3);
-  par[0]=P()->InnerFreonWidth()/2;par[1]=P()->FreonThickness()/2;par[2]=P()->InnerFreonLength()/2; 
-  gMC->Gsvolu("OQF2", "BOX ", idtmed[1007], par, 3);
-//Freon 
-  par[0]=P()->OuterFreonWidth()/2 - oquaThickness;
-  par[1]=P()->FreonThickness()/2;
-  par[2]=P()->OuterFreonLength()/2 - 2*oquaThickness; 
-  gMC->Gsvolu("FRE1", "BOX ", idtmed[1003], par, 3);
-
-  par[0]=P()->InnerFreonWidth()/2 - oquaThickness;
-  par[1]=P()->FreonThickness()/2;
-  par[2]=P()->InnerFreonLength()/2 - 2*oquaThickness; 
-  gMC->Gsvolu("FRE2", "BOX ", idtmed[1003], par, 3);    
-//Methane 
-  par[0]=pcX/2;par[1]=P()->GapThickness()/2;par[2]=pcY/2;         gMC->Gsvolu("META","BOX ",idtmed[1004], par, 3);
-//Methane gap 
-  par[0]=pcX/2;par[1]=P()->GapAmp()/2;par[2]=pcY/2;gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kGAP],par,3);
-//CsI PC
-  par[0]=pcX/2;par[1]=.25;par[2]=pcY/2;  gMC->Gsvolu("CSI ", "BOX ", (*fIdtmed)[kCSI], par, 3);
-//Anode grid 
-  par[0] = 0.;par[1] = .001;par[2] = 20.;  gMC->Gsvolu("GRID", "TUBE", idtmed[1006], par, 3);
-
-//Wire supports
-//Bar of metal
-  par[0]=pcX/2;par[1]=1.05;par[2]=1.05;  gMC->Gsvolu("WSMe", "BOX ", idtmed[1009], par, 3);
-//Ceramic pick up (base)
-  par[0]=pcX/2;par[1]= .25;par[2]=1.05;  gMC->Gsvolu("WSG1", "BOX ", idtmed[1010], par, 3);
-//Ceramic pick up (head)
-  par[0] = pcX/2;par[1] = .1;par[2] = .1;  gMC->Gsvolu("WSG2", "BOX ", idtmed[1010], par, 3);
-
-//Aluminium supports for methane and CsI
-//Short bar
-  par[0]=pcX/2;par[1]=P()->GapThickness()/2 + .25; par[2] = (68.35 - pcY/2)/2;
-  gMC->Gsvolu("SMSH", "BOX", idtmed[1009], par, 3);
-//Long bar
-  par[0]=(66.3 - pcX/2)/2;par[1]=P()->GapThickness()/2+.25;par[2]=pcY/2+68.35-pcY/2;
-  gMC->Gsvolu("SMLG", "BOX", idtmed[1009], par, 3);
-    
-//Aluminium supports for freon
-//Short bar
-  par[0] = P()->QuartzWidth()/2; par[1] = .3; par[2] = (68.35 - P()->QuartzLength()/2)/2;
-  gMC->Gsvolu("SFSH", "BOX", idtmed[1009], par, 3);    
-//Long bar
-  par[0] = (66.3 - P()->QuartzWidth()/2)/2; par[1] = .3;
-  par[2] = P()->QuartzLength()/2 + 68.35 - P()->QuartzLength()/2;
-  gMC->Gsvolu("SFLG", "BOX", idtmed[1009], par, 3);    
-//PCB backplane
-  par[0] = pcX/2;par[1] = .25; par[2] = pcY/4 -.5025;  gMC->Gsvolu("PCB ", "BOX", idtmed[1011], par, 3);
-
-//Backplane supports
-//Aluminium slab
-  par[0] = 33.15;par[1] = 2;par[2] = 21.65;  gMC->Gsvolu("BACK", "BOX", idtmed[1009], par, 3);    
-//Big hole
-  par[0] = 9.05; par[1] = 2; par[2] = 4.4625;  gMC->Gsvolu("BKHL", "BOX", idtmed[1000], par, 3);
-//Small hole
-  par[0] = 5.7;par[1] = 2;par[2] = 4.4625;  gMC->Gsvolu("BKHS", "BOX", idtmed[1000], par, 3);
-//Place holes inside backplane support
-  gMC->Gspos("BKHS", 1, "BACK", .8 + 5.7,0., .6 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 2, "BACK", -.8 - 5.7,0., .6 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 3, "BACK", .8 + 5.7,0., -.6 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 4, "BACK", -.8 - 5.7,0., -.6 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 5, "BACK", .8 + 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 6, "BACK", -.8 - 5.7,0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 7, "BACK", .8 + 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHS", 8, "BACK", -.8 - 5.7,0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 1, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., .6 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 2, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 3, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 4, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 5, "BACK", .8 + 11.4+ 1.6 + 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 6, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., .6 + 8.925 + 1.2 + 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 7, "BACK", .8 + 11.4 + 1.6 + 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
-  gMC->Gspos("BKHL", 8, "BACK", -.8 - 11.4 - 1.6 - 9.05, 0., -.6 - 8.925 - 1.2 - 4.4625, 0, "ONLY");
-//Place material inside RICH 
-  gMC->Gspos("SRIC", 1, "RICH", 0.,0., 0., 0, "ONLY");
-//  gMC->Gspos("AIR2", 1, "RICH", 66.3 + 1.2505, 1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
-//  gMC->Gspos("AIR2", 2, "RICH", -66.3 - 1.2505,1.276-P()->GapThickness()/2-P()->QuartzThickness()-P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, 0., 0, "ONLY");
-//  gMC->Gspos("AIR3", 1, "RICH", 0., 1.276-P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35, -68.35 - 1.25, 0, "ONLY");
-//  gMC->Gspos("AIR3", 2, "RICH", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.5 - 3.35,  68.35 + 1.25, 0, "ONLY");
-  gMC->Gspos("ALUM", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .376 -.025, 0., 0, "ONLY");
-  gMC->Gspos("HONE", 1, "SRIC", 0., 1.276- P()->GapThickness()/2  - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .05 - .188, 0., 0, "ONLY");
-  gMC->Gspos("ALUM", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .6 - .025, 0., 0, "ONLY");
-  gMC->Gspos("FOOT", 1, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 2, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3 , 36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 3, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 4, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, 36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 5, "SRIC", 64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 6, "SRIC", 21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 7, "SRIC", -21.65, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
-  gMC->Gspos("FOOT", 8, "SRIC", -64.95, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .4 - .3, -36.9, 0, "ONLY");
-  gMC->Gspos("OQUA", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()- .2, 0., 0, "ONLY");
-// Methane supports
-  gMC->Gspos("SMLG", 1, "SRIC", pcX/2 + (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY");
-  gMC->Gspos("SMLG", 2, "SRIC", - pcX/2 - (66.3 - pcX/2)/2, 1.276 + .25, 0., 0, "ONLY");
-  gMC->Gspos("SMSH", 1, "SRIC", 0., 1.276 + .25, pcY/2 + (68.35 - pcY/2)/2, 0, "ONLY");
-  gMC->Gspos("SMSH", 2, "SRIC", 0., 1.276 + .25, - pcY/2 - (68.35 - pcY/2)/2, 0, "ONLY");
-//Freon supports
-  Float_t suppY = 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness() - .2 + .3; //y position of freon supports
-  gMC->Gspos("SFLG", 1, "SRIC", P()->QuartzWidth()/2 + (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY");
-  gMC->Gspos("SFLG", 2, "SRIC", - P()->QuartzWidth()/2 - (66.3 - P()->QuartzWidth()/2)/2, suppY, 0., 0, "ONLY");
-  gMC->Gspos("SFSH", 1, "SRIC", 0., suppY, P()->QuartzLength()/2 + (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY");
-  gMC->Gspos("SFSH", 2, "SRIC", 0., suppY, - P()->QuartzLength()/2 - (68.35 - P()->QuartzLength()/2)/2, 0, "ONLY");
-  AliMatrix(idrotm[1019], 0., 0., 90., 0., 90., 90.);
-//Place spacers
-  Int_t nspacers = 30;
-  for (i = 0; i < nspacers/3; i++) {
-    zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE1", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings 
-  }
-  for (i = nspacers/3; i < (nspacers*2)/3; i++) {
-    zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE1", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings 
-  }
-  for (i = (nspacers*2)/3; i < nspacers; ++i) {
-    zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE1", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings  
-  }
-  for (i = 0; i < nspacers/3; i++) {
-    zs = -11.6/2 + (TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE2", 10.5, 0., zs, idrotm[1019], "ONLY");  //Original settings 
-  }
-  for (i = nspacers/3; i < (nspacers*2)/3; i++) {
-    zs = -11.6/2 + (nspacers/3 + TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE2", 0, 0., zs, idrotm[1019], "ONLY");  //Original settings 
-  }
-  for (i = (nspacers*2)/3; i < nspacers; ++i) {
-    zs = -11.6/2 + ((nspacers*2)/3 + TMath::Abs(nspacers/6) - i) * 12.2;
-    gMC->Gspos("SPAC", i, "FRE2", -10.5, 0., zs, idrotm[1019], "ONLY"); //Original settings  
-  }
-  gMC->Gspos("FRE1", 1, "OQF1", 0., 0., 0., 0, "ONLY");
-  gMC->Gspos("FRE2", 1, "OQF2", 0., 0., 0., 0, "ONLY");
-  gMC->Gspos("OQF1", 1, "SRIC", P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2 + 2, 1.276 - P()->GapThickness()/2- P()->QuartzThickness() -P()->FreonThickness()/2, 0., 0, "ONLY"); //Original settings (31.3)
-  gMC->Gspos("OQF2", 2, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY");          //Original settings 
-  gMC->Gspos("OQF1", 3, "SRIC", - (P()->OuterFreonWidth()/2 + P()->InnerFreonWidth()/2) - 2, 1.276 - P()->GapThickness()/2 - P()->QuartzThickness() - P()->FreonThickness()/2, 0., 0, "ONLY");       //Original settings (-31.3)
-  gMC->Gspos("QUAR", 1, "SRIC", 0., 1.276 - P()->GapThickness()/2 - P()->QuartzThickness()/2, 0., 0, "ONLY");
-  gMC->Gspos("GAP ", 1, "META", 0., P()->GapThickness()/2 - P()->GapAmp()/2 - 0.0001, 0., 0, "ONLY");
-  gMC->Gspos("META", 1, "SRIC", 0., 1.276, 0., 0, "ONLY");
-  gMC->Gspos("CSI ", 1, "SRIC", 0., 1.276 + P()->GapThickness()/2 + .25, 0., 0, "ONLY");
-//Wire support placing
-  gMC->Gspos("WSG2", 1, "GAP ", 0., P()->GapAmp()/2 - .1, 0., 0, "ONLY");
-  gMC->Gspos("WSG1", 1, "CSI ", 0., 0., 0., 0, "ONLY");
-  gMC->Gspos("WSMe", 1, "SRIC ", 0., 1.276 + P()->GapThickness()/2 + .5 + 1.05, 0., 0, "ONLY");
-//Backplane placing
-  gMC->Gspos("BACK", 1, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 43.3, 0, "ONLY");
-  gMC->Gspos("BACK", 2, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2 , 43.3, 0, "ONLY");
-  gMC->Gspos("BACK", 3, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
-  gMC->Gspos("BACK", 4, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, 0., 0, "ONLY");
-  gMC->Gspos("BACK", 5, "SRIC ", 33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
-  gMC->Gspos("BACK", 6, "SRIC ", -33.15, 1.276 + P()->GapThickness()/2 + .5 + 2.1 + 2, -43.3, 0, "ONLY");
-//PCB placing
-  gMC->Gspos("PCB ", 1, "SRIC ", 0.,  1.276 + P()->GapThickness()/2 + .5 + 1.05, pcX/4 + .5025 + 2.5, 0, "ONLY");
-  gMC->Gspos("PCB ", 2, "SRIC ", 0.,  1.276 + P()->GapThickness()/2 + .5 + 1.05, -pcX/4 - .5025 - 2.5, 0, "ONLY");
-  
+  Int_t id=0;
+       
 //place chambers into mother volume ALIC
-  for(int i=1;i<=kNchambers;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,"ALIC",C(i)->X(),C(i)->Y(),C(i)->Z(),idrotm[1000+i], "ONLY");
-  }
+  par[0]=(6*mm+1681*mm+6*mm)/2;par[1]=(6*mm+1466*mm+6*mm)/2;par[2]=(80*mm+40*mm)*2/2;gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3);//2033P1  z p84 TDR
+  if(P()->IsRadioSrc()){
+    AliInfo("Special test beam geometry");
+    gMC->Gspos("RICH",1,"ALIC",0,0,0,0, "ONLY");  //single RICH chamber without rotation test beam geometry
+  }else
+    for(int i=1;i<=kNchambers;i++){ //full RICH with 7 chambers
+      AliMatrix(id,C(i)->ThetaXd(),C(i)->PhiXd(),  C(i)->ThetaYd(),C(i)->PhiYd(),  C(i)->ThetaZd(),C(i)->PhiZd());
+      gMC->Gspos("RICH",i,"ALIC",C(i)->Center().X(),C(i)->Center().Y(),C(i)->Center().Z(),id, "ONLY");
+    }
+//Pad Panel frame  6 sectors
+  par[0]=648*mm/2;par[1]=  411*mm/2;par[2]=40  *mm/2;gMC->Gsvolu("RPPF","BOX ",(*fIdtmed)[kAl]  ,par,3);//PPF 2001P2 inner size of the slab by 1mm more
+  par[0]=181*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFL","BOX ",(*fIdtmed)[kAir] ,par,3);//large whole
+  par[0]=114*mm/2;par[1]=89.25*mm/2;par[2]=38.3*mm/2;gMC->Gsvolu("PPFS","BOX ",(*fIdtmed)[kAir] ,par,3);//small whole
+  par[0]=644*mm/2;par[1]=  407*mm/2;par[2]= 1.7*mm/2;gMC->Gsvolu("RPC ","BOX ",(*fIdtmed)[kCsI] ,par,3);//by 0.2 mm more then actual size (PCB 2006P1)
+  
+  gMC->Gspos("RPPF",1,"RICH",    -335*mm,      -433*mm,  8*cm+20*mm,  0,"ONLY");//F1 2040P1 z p.84 TDR
+  gMC->Gspos("RPPF",2,"RICH",    +335*mm,      -433*mm,  8*cm+20*mm,  0,"ONLY");
+  gMC->Gspos("RPPF",3,"RICH",    -335*mm,         0*mm,  8*cm+20*mm,  0,"ONLY");
+  gMC->Gspos("RPPF",4,"RICH",    +335*mm,         0*mm,  8*cm+20*mm,  0,"ONLY");
+  gMC->Gspos("RPPF",5,"RICH",    -335*mm,      +433*mm,  8*cm+20*mm,  0,"ONLY");
+  gMC->Gspos("RPPF",6,"RICH",    +335*mm,      +433*mm,  8*cm+20*mm,  0,"ONLY");  
+  gMC->Gspos("RPC ",1,"RPPF",       0*mm,         0*mm,   -19.15*mm,  0,"ONLY");//PPF 2001P2 
+  gMC->Gspos("PPFL",1,"RPPF",  -224.5*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",2,"RPPF",  -224.5*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",3,"RPPF",  -224.5*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",4,"RPPF",  -224.5*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",1,"RPPF",  - 65.0*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",2,"RPPF",  - 65.0*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",3,"RPPF",  - 65.0*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",4,"RPPF",  - 65.0*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",5,"RPPF",  + 65.0*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",6,"RPPF",  + 65.0*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",7,"RPPF",  + 65.0*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFS",8,"RPPF",  + 65.0*mm,  +151.875*mm,     0.85*mm,  0,"ONLY"); 
+  gMC->Gspos("PPFL",5,"RPPF",  +224.5*mm,  -151.875*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",6,"RPPF",  +224.5*mm,  - 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",7,"RPPF",  +224.5*mm,  + 50.625*mm,     0.85*mm,  0,"ONLY");
+  gMC->Gspos("PPFL",8,"RPPF",  +224.5*mm,  +151.875*mm,     0.85*mm,  0,"ONLY");
+//Gap - anod wires 6 copies to RICH
+  par[0]=648*mm/2;par[1]=  411*mm/2 ;par[2]=4.45*mm/2;gMC->Gsvolu("RGAP ","BOX ",(*fIdtmed)[kCH4] ,par,3);//xy as PPF 2001P2 z WP 2099P1
+  par[0]=  0*mm  ;par[1]=  20*mkm/2 ;par[2]= 648*mm/2;gMC->Gsvolu("RANO","TUBE",(*fIdtmed)[kW]   ,par,3);//WP 2099P1 z = gap x PPF 2001P2
+  AliMatrix(id,180,0, 90,90, 90,0); //wires along x
+  
+  gMC->Gspos("RGAP",1,"RICH",    -335*mm,      -433*mm,8*cm-2.225*mm, 0,"ONLY"); //F1 2040P1 z WP 2099P1
+  gMC->Gspos("RGAP",2,"RICH",    +335*mm,      -433*mm,8*cm-2.225*mm, 0,"ONLY"); 
+  gMC->Gspos("RGAP",3,"RICH",    -335*mm,         0*mm,8*cm-2.225*mm, 0,"ONLY"); 
+  gMC->Gspos("RGAP",4,"RICH",    +335*mm,         0*mm,8*cm-2.225*mm, 0,"ONLY"); 
+  gMC->Gspos("RGAP",5,"RICH",    -335*mm,      +433*mm,8*cm-2.225*mm, 0,"ONLY"); 
+  gMC->Gspos("RGAP",6,"RICH",    +335*mm,      +433*mm,8*cm-2.225*mm, 0,"ONLY"); 
+  for(int i=1;i<=96;i++)
+    gMC->Gspos("RANO",i,"RGAP",     0*mm, -411/2*mm+i*4*mm, 0.185*mm, id,"ONLY"); //WP 2099P1  
+//Radiator 3 modules
+  par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=  24*mm/2;  gMC->Gsvolu("RRAD","BOX ",(*fIdtmed)[kC6F14]     ,par,3); // Rad 2011P1
+  par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=   4*mm/2;  gMC->Gsvolu("RRFR","BOX ",(*fIdtmed)[kRoha]      ,par,3); //front 
+  par[0]=1330*mm/2 ;par[1]= 413*mm/2  ;par[2]=   5*mm/2;  gMC->Gsvolu("RRWI","BOX ",(*fIdtmed)[kSiO2]      ,par,3); //window
+  par[0]=1330*mm/2 ;par[1]=   5*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RRLO","BOX ",(*fIdtmed)[kRoha]      ,par,3); //long side  
+  par[0]=  10*mm/2 ;par[1]= 403*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RRSH","BOX ",(*fIdtmed)[kRoha]      ,par,3); //short side 
+  par[0]=   0      ;par[1]=  10*mm/2  ;par[2]=  15*mm/2;  gMC->Gsvolu("RRSP","TUBE",(*fIdtmed)[kSiO2]      ,par,3); //spacer        
+    
+  gMC->Gspos("RRAD",1,"RICH",   0*mm,-434*mm,   -12*mm,  0,"ONLY"); //radiator to RICH
+  gMC->Gspos("RRAD",2,"RICH",   0*mm,   0*mm,   -12*mm,  0,"ONLY"); 
+  gMC->Gspos("RRAD",3,"RICH",   0*mm,+434*mm,   -12*mm,  0,"ONLY"); 
+    gMC->Gspos("RRFR",1,"RRAD",   0*mm,   0*mm, -10.0*mm,  0,"ONLY"); //front cover 
+    gMC->Gspos("RRWI",1,"RRAD",   0*mm,   0*mm,   9.5*mm,  0,"ONLY"); //quartz window (back cover)
+    gMC->Gspos("RRLO",1,"RRAD",   0*mm,-204*mm,  -0.5*mm,  0,"ONLY"); //long side
+    gMC->Gspos("RRLO",2,"RRAD",   0*mm,+204*mm,  -0.5*mm,  0,"ONLY"); //long side
+    gMC->Gspos("RRSH",1,"RRAD",-660*mm,   0*mm,  -0.5*mm,  0,"ONLY"); //short side
+    gMC->Gspos("RRSH",2,"RRAD",+660*mm,   0*mm,  -0.5*mm,  0,"ONLY"); //short side 
+    for(int i=0;i<3;i++)
+      for(int j=0;j<10;j++)
+        gMC->Gspos("RRSP",10*i+j,"RRAD",-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm,0,"ONLY");
 
-         
-  if(GetDebug())Info("CreateGeometry","Stop main.");  
+//Sandbox  
+  par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]=50.5*mm/2; gMC->Gsvolu("RSNB","BOX ",(*fIdtmed)[kAir]  ,par,3);  //2072P1   
+  par[0]=1419*mm/2 ;par[1]=1378*mm/2;par[2]= 0.5*mm/2; gMC->Gsvolu("RSCO","BOX ",(*fIdtmed)[kAl]   ,par,3);  //cover
+  par[0]=1359*mm/2 ;par[1]=1318*mm/2;par[2]=49.5*mm/2; gMC->Gsvolu("RSHO","BOX ",(*fIdtmed)[kRoha] ,par,3); //honeycomb structure 
+  
+  gMC->Gspos("RSNB",1,"RICH",   0*mm, 0*mm, -73.75*mm, 0,"ONLY"); //p.84 TDR sandbox to rich
+    gMC->Gspos("RSHO",1,"RSNB", 0*mm, 0*mm,      0*mm, 0,"ONLY"); //2072P1 honeycomv to sandbox
+    gMC->Gspos("RSCO",1,"RSNB", 0*mm, 0*mm,    +25*mm, 0,"ONLY"); //cover to sandbox
+    gMC->Gspos("RSCO",2,"RSNB", 0*mm, 0*mm,    -25*mm, 0,"ONLY"); //cover to sandbox
+             
+  AliDebug(1,"Stop main.");  
 }//void AliRICH::CreateGeometry()
 //__________________________________________________________________________________________________
-void AliRICH::Reconstruct()const
-{
-  AliRICHClusterFinder finder(const_cast<AliRICH*>(this));
-  finder.Exec();
-}
-//__________________________________________________________________________________________________
 void AliRICH::ControlPlots()
 { 
-//Creates a set of hists to control the results of simulation
+// Creates a set of hists to control the results of simulation. Hists are in file $HOME/RCP.root
      
   TH1F *pHxD=0,*pHyD=0,*pNumClusH1=0,
                    *pQdcH1=0,       *pSizeH1=0,
@@ -667,14 +560,14 @@ void AliRICH::ControlPlots()
   Bool_t isDig =!GetLoader()->LoadDigits();
   Bool_t isClus=!GetLoader()->LoadRecPoints();
 
-  if(!isDig && !isClus){Error("ControlPlots","No digits and clusters! Nothing to do.");return;}
+  if(!isDig && !isClus){AliError("No digits and clusters! Nothing to do.");return;}
   
   TStopwatch sw;TDatime time;
     
   TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE");   
   
   if(isDig){
-    cout<<"Digits available\n";
+    AliInfo("Digits available");
     pHxD=new TH1F("HitDigitDiffX","Hit-Digits diff X all chambers;diff [cm]",100,-10,10); 
     pHyD=new TH1F("HitDigitDiffY","Hit-Digits diff Y all chambers;diff [cm]",100,-10,10); 
   }//isDig
@@ -742,8 +635,8 @@ void AliRICH::ControlPlots()
       if(isDig){
         for(Int_t iDigN=0;iDigN<R()->Digits(iChamN)->GetEntries();iDigN++){//digits loop
           AliRICHdigit *pDig=(AliRICHdigit*)R()->Digits(iChamN)->At(iDigN);
-          AliRICHhit   *pHit=Hit(pDig->GetTrack(0));
-          TVector2 hitV2=R()->C(iChamN)->Glob2Loc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad());
+          AliRICHhit   *pHit=Hit(pDig->GetTrack(0));//get first hit of this digit
+          TVector2 hitV2=R()->C(iChamN)->Mrs2Pc(pHit->OutX3()); TVector2 digV2=R()->P()->Pad2Loc(pDig->Pad());//center of pad for digit
           pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y());
         }//digits loop
       }//isDig
@@ -776,7 +669,7 @@ AliRICHhit* AliRICH::Hit(Int_t tid)
 void AliRICH::PrintHits(Int_t iEvtN)
 {
 //Prints a list of RICH hits for a given event. Default is event number 0.
-  Info("PrintHits","List of RICH hits for event %i",iEvtN);
+  AliInfo(Form("List of RICH hits for event %i",iEvtN));
   R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
   if(R()->GetLoader()->LoadHits()) return;
   
@@ -787,7 +680,7 @@ void AliRICH::PrintHits(Int_t iEvtN)
     iTotalHits+=R()->Hits()->GetEntries();
   }
   R()->GetLoader()->UnloadHits();
-  Info("PrintHits","totally %i hits",iTotalHits);
+  AliInfo(Form("totally %i hits",iTotalHits));
 }
 //__________________________________________________________________________________________________
 void AliRICH::PrintSDigits(Int_t iEvtN)
@@ -837,40 +730,13 @@ void AliRICH::PrintClusters(Int_t iEvtN)
   Info("PrintClusters","totally %i clusters",iTotalClusters);
 }
 //__________________________________________________________________________________________________
-void AliRICH::FillESD(AliESD *pESD)const
+void AliRICH::PrintTracks(Int_t iEvtN)
 {
-//This methode fills AliESDtrack with information from RICH  
-  Info("FillESD","Start with %i tracks",pESD->GetNumberOfTracks());
-  const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton
-  const Double_t refIndex = 1.29052;
-
-  Double_t thetaExp = 0.7;
-  Double_t thetaTh[5];
-  Double_t sinThetaThNorm;
-  Double_t sigmaThetaTh[5];
-  Double_t height[5];
-  Double_t totalHeight=0;
-  
-  for(Int_t iTrackN=0;iTrackN<pESD->GetNumberOfTracks();iTrackN++){//ESD tracks loop
-    AliESDtrack *pTrack = pESD->GetTrack(iTrackN);
-    
-    pTrack->Print("");
-//    TVector2 x2=P()->HelixCross(pTrack);//returns cross point of track with RICH PC in LRS
-    Double_t pmod = pTrack->GetP();
-    
-    for(Int_t iPart=4;iPart>=0;iPart--){
-      Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod);
-      if(cosThetaTh>=1) {break;}
-      thetaTh[iPart] = TMath::ACos(cosThetaTh);
-      sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex));
-      sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
-      height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]);
-      totalHeight +=height[iPart];
-    }
-    
-    pTrack->SetRICHsignal(thetaExp);
-    Double_t richPID[5];
-    for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;    
-    pTrack->SetRICHpid(richPID); 
-  }//ESD tracks loop
+//prints a list of tracks (including secondaru) for a given event
+  Info("PrintTracks","List of all tracks for event %i",iEvtN);
+  R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
+  if(R()->GetLoader()->GetRunLoader()->LoadHeader()) return;
+  if(R()->GetLoader()->GetRunLoader()->LoadKinematics()) return;
+  Int_t iTracksCounter=0;
+  Info("PrintTracks","totally %i tracks",iTracksCounter);
 }
index bba4b60..6ab61c6 100644 (file)
@@ -77,9 +77,9 @@ class AliRICHcluster :public TObject
 public:
   enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved,kEmpty=kBad};
                     AliRICHcluster():TObject(),fCFM(0),fSize(0),fShape(0),fQdc(0),fChamber(0),fX(0),fY(0),fStatus(kEmpty),fDigits(0) {}    
-  virtual          ~AliRICHcluster()                 {Reset();}  
+  virtual          ~AliRICHcluster()                 {AliDebug(1,"Start");/*Reset();*/}  
          void       Reset()                          {DeleteDigits();fCFM=fSize=fShape=fQdc=fChamber=0;fX=fY=0;fStatus=kEmpty;} //cleans the cluster
-         void       DeleteDigits()                   {if(fDigits) delete fDigits; fDigits=0;}           //deletes the list of digits  
+         void       DeleteDigits()                   {if(fDigits) {delete fDigits;} fDigits=0;}           //deletes the list of digits  
   AliRICHcluster&   operator=(const AliRICHcluster&) {return *this;}                                 
          Int_t      Nlocals()                   const{return fSize-10000*(fSize/10000);}                //number of local maximums
          Int_t      Size()                      const{return fSize/10000;}                              //number of digits in cluster
@@ -108,11 +108,12 @@ public:
          Int_t      CombiPid()                  const{return fCFM;}                                        //
          void       CFM(Int_t c,Int_t f,Int_t m)     {fCFM=1000000*c+1000*f+m;}                            //cluster contributors
          TObjArray* Digits()                    const{return fDigits;}                                     //  
-  virtual void      Print(Option_t *option="")const;                                                                //
-  inline  void      AddDigit(AliRICHdigit *pDig);                                                                   //
-  inline  void      CoG(Int_t nLocals);                                                                             //calculates center of gravity
-          void      Fill(AliRICHcluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm)                           //form new resolved cluster from raw one
-                    {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} //
+  virtual void      Print(Option_t *option="")const;                                                   //
+  inline  void      AddDigit(AliRICHdigit *pDig);                                                      //
+  inline  void      CoG(Int_t nLocals);                                                                //calculates center of gravity
+          void      Fill(AliRICHcluster *pRaw,Double_t x,Double_t y,Double_t q,Int_t cfm)              //form new resolved cluster from raw one
+                    {fCFM=cfm;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} 
+         Double_t   DistTo(TVector2 x)          const{return TMath::Sqrt((x.X()-fX)*(x.X()-fX)+(x.Y()-fY)*(x.Y()-fY));} //distance to given point 
 protected:
   Int_t         fCFM;         //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
   Int_t         fSize;        //10000*(how many digits belong to this cluster) + nLocalMaxima     
@@ -171,9 +172,8 @@ public:
   virtual void          StepManager()                               =0;                                  //interface from AliMC
   virtual void          Hits2SDigits();                                                                  //interface from AliSimulation
   virtual AliDigitizer* CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);}  //interface from AliSimulation
-  virtual void          Reconstruct()                         const;                                     //interface from AliReconstruction
-  virtual void          FillESD(AliESD *pESD)                 const;                                     //interface from AliReconstruction          
-  virtual void          Print(Option_t *option)               const;                                     //prints current RICH status
+//  virtual void          Reconstruct()                         const;                                     //interface from AliReconstruction
+//  virtual void          FillESD(AliESD *pESD)                 const;                                     //interface from AliReconstruction          
   virtual void          SetTreeAddress();                                                                //interface from AliLoader
   virtual void          MakeBranch(Option_t *opt=" ");                                                   //interface from AliLoader
   virtual void          CreateMaterials();                                                               //interface from AliMC
@@ -198,10 +198,12 @@ public:
   AliRICH*        R()                        {return this;}                                 //provides pointer to RICH main object
   TVector         Counters()            const{return fCounters;}                            //provides a set of counters
   void            ControlPlots();                                                           //utility
+  virtual void    Print(Option_t *option="")               const;                           //prints current RICH status
   void            PrintHits    (Int_t iEvent=0);                                            //utility
   void            PrintSDigits (Int_t iEvent=0);                                            //utility
   void            PrintDigits  (Int_t iEvent=0);                                            //utility
   void            PrintClusters(Int_t iEvent=0);                                            //utility
+  void            PrintTracks  (Int_t iEvent=0);                                            //utility
             
   void AddHit(Int_t c,Int_t tid,TVector3 i3,TVector3 o3,Double_t eloss=0){TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(c,tid,i3,o3,eloss);}
   inline void AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid); 
@@ -213,7 +215,7 @@ public:
        {Int_t c=cl.C()-1;TClonesArray &tmp=*((TClonesArray*)fClusters->At(c));new(tmp[fNclusters[c]++])AliRICHcluster(cl);}
   AliRICHhit* Hit(Int_t tid);           //returns pointer ot RICH hit for a given tid
 protected:  
-  enum                  {kAir=1,kCSI=6,kGAP=9,kAl=10,kCH4=5,kSteel=15,kPerpex=16,kSr90=17};
+  enum                  {kAir=1,kRoha,kSiO2,kC6F14,kCH4,kCsI,kGridCu,kOpSiO2,kGap,kAl,kGlass,kCu,kW,kSteel,kPerpex,kSr90};
   AliRICHParam         *fpParam;             //main RICH parametrization     
                                              //fHits and fDigits belong to AliDetector
   TClonesArray         *fSdigits;            //! list of sdigits  
@@ -233,21 +235,21 @@ protected:
 void AliRICH::CreateHits()
 {
   if(fHits) return;
-  if(GetDebug())Info("CreateHits","creating hits container.");
+  AliDebug(1,"creating hits container.");
   fHits=new TClonesArray("AliRICHhit",10000);   fNhits=0;
 }
 //__________________________________________________________________________________________________
 void AliRICH::CreateSDigits()
 {
   if(fSdigits) return;
-  if(GetDebug())Info("CreateSDigits","creating sdigits container.");
+  AliDebug(1,"creating sdigits container.");
   fSdigits=new TClonesArray("AliRICHdigit",10000); fNsdigits=0;
 }
 //__________________________________________________________________________________________________
 void AliRICH::CreateDigits()
 {
   if(fDigitsNew) return;
-  if(GetDebug())Info("CreateDigits","creating digits containers.");
+  AliDebug(1,"creating digits containers.");
   fDigitsNew = new TObjArray(kNchambers);  
   for(Int_t i=0;i<kNchambers;i++) {fDigitsNew->AddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;}
 }
@@ -255,7 +257,7 @@ void AliRICH::CreateDigits()
 void AliRICH::CreateClusters()
 {
   if(fClusters) return;
-  if(GetDebug())Info("CreateClusters","creating clusters containers.");
+  AliDebug(1,"creating clusters containers.");
   fClusters = new TObjArray(kNchambers);  
   for(Int_t i=0;i<kNchambers;i++) {fClusters->AddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;}
 }
index 7f7f86f..c945d3a 100644 (file)
 
 #include "AliRICHChamber.h"
 #include <TRotMatrix.h>
+#include <AliLog.h>
 
 ClassImp(AliRICHChamber)       
 //______________________________________________________________________________
 AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed()
 {
 //main ctor. Defines all geometry parameters for the given module.
-  SetToZenith();//put to up position   
+// 7 6     ^      
+// 5 4 3   |      
+//   2 1   |
+// <-----z y   . x     
+//  horizontal angle between chambers  19.5 grad
+//  vertical angle between chambers    20   grad     
+  RotY(90);//rotate around y
+  fCenterX3.SetXYZ(490,0,0);fPcX3.SetXYZ(490+8,0,0);fRadX3.SetXYZ(490-2,0,0); //shift center along x by 490 cm
+  
   switch(iChamber){
-    case 1:
-      RotateX(-AliRICHParam::AngleYZ());   
-      RotateZ(-AliRICHParam::AngleXY());      
-      fName="RICHc1";fTitle="RICH chamber 1";
+    case 0:                    //special test beam configuration without rotation.
+      break;
+    case 1:        
+      RotY(19.5); RotZ(-20);   //right and down 
       break;      
     case 2:
-      RotateZ(-AliRICHParam::AngleXY());      
-      fName="RICHc2";fTitle="RICH chamber 2";
+      RotZ(-20);              //down
       break;      
     case 3:
-      RotateX(-AliRICHParam::AngleYZ());
-      fName="RICHc3";fTitle="RICH chamber 3";
+      RotY(19.5);             //right 
       break;      
     case 4:          
-      fName="RICHc4";fTitle="RICH chamber 4";  //no turns
-      break;      
+      break;                  //no rotation
     case 5:
-      RotateX(AliRICHParam::AngleYZ());
-      fName="RICHc5";fTitle="RICH chamber 5";
+      RotY(-19.5);            //left   
       break;      
     case 6:
-      RotateZ(AliRICHParam::AngleXY());      
-      fName="RICHc6";fTitle="RICH chamber 6";
+      RotZ(20);               //up
       break;      
     case 7:
-      RotateX(AliRICHParam::AngleYZ());            
-      RotateZ(AliRICHParam::AngleXY());      
-      fName="RICHc7";fTitle="RICH chamber 7";
+      RotY(-19.5); RotZ(20);  //left and up 
       break;      
     default:
       Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iChamber);
   }//switch(iModuleN)
-  RotateZ(AliRICHParam::AngleRot());//apply common rotation  
+  fName=Form("RICHc%i",iChamber);fTitle=Form("RICH chamber %i",iChamber);
+  RotZ(30);     //apply common rotation  
   fpRotMatrix=new TRotMatrix("rot"+fName,"rot"+fName, Rot().ThetaX()*TMath::RadToDeg(), Rot().PhiX()*TMath::RadToDeg(),
                                                       Rot().ThetaY()*TMath::RadToDeg(), Rot().PhiY()*TMath::RadToDeg(),
                                                       Rot().ThetaZ()*TMath::RadToDeg(), Rot().PhiZ()*TMath::RadToDeg());
@@ -64,9 +67,6 @@ AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed()
 void AliRICHChamber::Print(Option_t *opt) const
 {
 // Debug printout
-//  printf("%s r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f  %6.2f,%6.2f %6.2f,%6.2f %6.2f,%6.2f\n",fName.Data(),
-//                     Rho(), ThetaD(),PhiD(),   X(),    Y(),    Z(),
-//                     ThetaXd(),PhiXd(),ThetaYd(),PhiYd(),ThetaZd(),PhiZd());
-  fCenterV3.Print(opt);fPcX3.Print(opt);
+  ToAliInfo(fCenterX3.Print(opt));
 }//Print()
 //__________________________________________________________________________________________________
index 08ebdb9..2942481 100644 (file)
@@ -20,49 +20,43 @@ public:
            AliRICHChamber(const AliRICHChamber &chamber):TNamed(chamber) {;}
   virtual ~AliRICHChamber()                                              {;}
            AliRICHChamber& operator=(const AliRICHChamber&)              {return *this;}
-  
+
+  static Double_t AlphaFeedback(Int_t )      {return 0.030;}                              //determines number of feedback photons
   TRotMatrix* RotMatrix()          const{return fpRotMatrix;}
-  TString RotMatrixName()          const{return "rot"+fName;}
+  TString     RotMatrixName()      const{return "rot"+fName;}
   TRotation   Rot()                const{return fRot;}
-  Double_t    Rho()                const{return fCenterV3.Mag();}                                //gives  distance to chamber center in MRS
-  Double_t    ThetaD()             const{return fCenterV3.Theta()*TMath::RadToDeg();}            //gives polar angle of chamber center in MRS
-  Double_t    PhiD()               const{return fCenterV3.Phi()  *TMath::RadToDeg();}            //gives azimuthal angle of chamber center in MRS
+  Double_t    Rho()                const{return fCenterX3.Mag();}                                //gives  distance to chamber center in MRS
+  Double_t    ThetaD()             const{return fCenterX3.Theta()*TMath::RadToDeg();}            //gives polar angle of chamber center in MRS
+  Double_t    PhiD()               const{return fCenterX3.Phi()  *TMath::RadToDeg();}            //gives azimuthal angle of chamber center in MRS
   Double_t    ThetaXd()            const{return fRot.ThetaX()    *TMath::RadToDeg();}    
   Double_t    PhiXd()              const{return fRot.PhiX()      *TMath::RadToDeg();}    
   Double_t    ThetaYd()            const{return fRot.ThetaY()    *TMath::RadToDeg();}    
   Double_t    PhiYd()              const{return fRot.PhiY()      *TMath::RadToDeg();}    
   Double_t    ThetaZd()            const{return fRot.ThetaZ()    *TMath::RadToDeg();}    
   Double_t    PhiZd()              const{return fRot.PhiZ()      *TMath::RadToDeg();}    
-  void        RotateX(Double_t a)       {fRot.RotateX(a);fCenterV3.RotateX(a);fPcX3.RotateX(a);} //rotate chamber around X by "a" degrees
-  void        RotateY(Double_t a)       {fRot.RotateY(a);fCenterV3.RotateY(a);fPcX3.RotateY(a);} //rotate chamber around Y by "a" degrees
-  void        RotateZ(Double_t a)       {fRot.RotateZ(a);fCenterV3.RotateZ(a);fPcX3.RotateZ(a);} //rotate chamber around Z by "a" degrees
-  Double_t    X()                  const{return fCenterV3.X();}  
-  Double_t    Y()                  const{return fCenterV3.Y();}   
-  Double_t    Z()                  const{return fCenterV3.Z();}
-  TVector2    Glob2Loc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(x3.Z()+0.5*AliRICHParam::PcSizeX(),-x3.X()+0.5*AliRICHParam::PcSizeY());}//Y and Z are misplaced?????
-  TVector3    Loc2Glob(TVector2 x2)const{TVector3 x3(-x2.Y()+0.5*AliRICHParam::PcSizeY(),0,x2.X()-0.5*AliRICHParam::PcSizeX());x3.Transform(fRot); x3+=fPcX3;return x3;}
-  
-  TVector2    Glob2Loc(TLorentzVector x4)            const{return Glob2Loc(x4.Vect());}
-  
-  void        Print(Option_t *sOption)const;//virtual      
-   
-
-  inline void SetToZenith();
-  TRotMatrix *GetRotMatrix()       const{return fpRotMatrix;}  
+  void        RotX(Double_t a)       {a*=TMath::DegToRad();fRot.RotateX(a);fCenterX3.RotateX(a);fRadX3.RotateX(a);fPcX3.RotateX(a);}//degrees around X
+  void        RotY(Double_t a)       {a*=TMath::DegToRad();fRot.RotateY(a);fCenterX3.RotateY(a);fRadX3.RotateY(a);fPcX3.RotateY(a);}//degrees around Y
+  void        RotZ(Double_t a)       {a*=TMath::DegToRad();fRot.RotateZ(a);fCenterX3.RotateZ(a);fRadX3.RotateZ(a);fPcX3.RotateZ(a);}//degrees around Z
+  TVector3    Rad()               const{return fRadX3;}         //provides center of radiator position in MRS, cm   
+  TVector3    Pc()                const{return fPcX3;}          //provides center of photocathond position in MRS, cm
+  TVector3    Center()            const{return fCenterX3;}      //provides center of chamber position in MRS, cm
+  void        Print(Option_t *sOption)const;                    //virtual interface from TObject
+//Transformations for photcathode plane  
+  TVector2    Mrs2Pc(TVector3 x3)const{x3-=fPcX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());}
+  TVector3    Pc2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fPcX3;return x3;}  
+  TVector2    Mrs2Pc(TLorentzVector x4)            const{return Mrs2Pc(x4.Vect());}
+//Transformations for radiator plane  
+  TVector2    Mrs2Rad(TVector3 x3)const{x3-=fRadX3;x3.Transform(fRot.Inverse());return TVector2(-x3.X()+0.5*AliRICHParam::PcSizeX(),x3.Y()+0.5*AliRICHParam::PcSizeY());}
+  TVector3    Rad2Mrs(TVector2 x2)const{TVector3 x3(-x2.X()+0.5*AliRICHParam::PcSizeX(),x2.Y()-0.5*AliRICHParam::PcSizeY(),0);x3.Transform(fRot); x3+=fRadX3;return x3;}  
+  TVector3    PMrs2Loc(TVector3 p3)const{TVector3 ploc=Rot().Invert()*p3;ploc.SetXYZ(-ploc.Px(),ploc.Py(),ploc.Pz()); return ploc;}  
 protected:
-  TVector3      fCenterV3;        //chamber center position in MRS (cm) 
+  TVector3      fCenterX3;        //chamber center position in MRS (cm) 
+  TVector3      fRadX3;           //radiator entrance center position in MRS (cm)
   TVector3      fPcX3;            //PC center position in MRS (cm)
   TRotation     fRot;             //chamber rotation in MRS
   TRotMatrix   *fpRotMatrix;      //rotation matrix of the chamber with respect to MRS 
-  ClassDef(AliRICHChamber,6)      //single RICH chamber description
+  ClassDef(AliRICHChamber,7)      //single RICH chamber description
 };//class AliRICHChamber
-//__________________________________________________________________________________________________
-void AliRICHChamber::SetToZenith()
-{
-//Put the chamber to zenith. Position of PC is shifted in X-Z plane since the origin of chamber local system is in
-//left hand down coner.     
-  fCenterV3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2           ,0); 
-      fPcX3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2+5.276+0.25,0);   
-}
-//__________________________________________________________________________________________________
+
 #endif //AliRICHChamber_h
index c12ff0e..e22e715 100644 (file)
@@ -29,33 +29,33 @@ ClassImp(AliRICHClusterFinder)
 //__________________________________________________________________________________________________
 AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)   
 {//main ctor
-  fRICH = pRICH;//this goes first as GetDebug depends on it.
-  if(GetDebug())Info("main ctor","Start.");
+  fRICH = pRICH;
+  AliDebug(1,"main ctor Start.");
   
   fDigitMap = 0;
   fRawCluster.Reset();
   fResolvedCluster.Reset();
-  if(GetDebug())Info("main ctor","Stop.");
+  AliDebug(1,"main ctor Stop.");
 }//main ctor
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::Exec()
 {
 //Main method of cluster finder. Loops on  events and chambers, everything else is done in FindClusters()  
-  if(GetDebug()) Info("Exec","Start.");
+  AliDebug(1,"Exec Start.");
     
   R()->GetLoader()                ->LoadDigits();   
-  R()->GetLoader()->GetRunLoader()->LoadHeader();
-  R()->GetLoader()->GetRunLoader()->LoadKinematics();
+//  R()->GetLoader()->GetRunLoader()->LoadHeader(); 
+  R()->GetLoader()->GetRunLoader()->LoadKinematics(); //header is already loaded
 
   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
-    if(GetDebug()) Info("Exec","Processing event %i...",iEventN);
+    AliDebug(1,Form("Processing event %i...",iEventN));
     R()->GetLoader()->GetRunLoader()->GetEvent(iEventN);
     
     R()->GetLoader()->MakeTree("R");  R()->MakeBranch("R");
     R()->ResetDigits();               R()->ResetClusters();
     
     R()->GetLoader()->TreeD()->GetEntry(0);
-    for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop
+    for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){//chambers loop
       FindClusters(iChamber);
     }//chambers loop
     R()->GetLoader()->TreeR()->Fill();  R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event
@@ -65,17 +65,17 @@ void AliRICHClusterFinder::Exec()
   R()->ResetClusters();
   R()->GetLoader()                ->UnloadDigits(); 
   R()->GetLoader()                ->UnloadRecPoints();  
-  R()->GetLoader()->GetRunLoader()->UnloadHeader();
+//  R()->GetLoader()->GetRunLoader()->UnloadHeader();
   R()->GetLoader()->GetRunLoader()->UnloadKinematics();
 
-  if(GetDebug()) Info("Exec","Stop.");      
+  AliDebug(1,"Stop.");      
 }//Exec()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FindClusters(Int_t iChamber)
 {
 //Loops on digits for a given chamber, forms raw clusters, then tries to resolve them if requested
   Int_t iNdigits=R()->Digits(iChamber)->GetEntriesFast();
-  if(GetDebug())Info("FindClusters","Start for chamber %i with %i digits.",iChamber,iNdigits);  
+  AliDebug(1,Form("Start for chamber %i with %i digits.",iChamber,iNdigits));  
   
   if(iNdigits==0)return;//no digits for a given chamber, nothing to do
 
@@ -87,9 +87,9 @@ void AliRICHClusterFinder::FindClusters(Int_t iChamber)
     if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit
        
     FormRawCluster(i,j);//form raw cluster starting from (i,j) pad 
-    if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();}  
+    AliDebug(1,"After FormRawCluster:");ToAliDebug(1,fRawCluster.Print());  
     FindLocalMaxima();  //find number of local maxima and initial center of gravity
-    if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();}  
+    AliDebug(1,"After FindLocalMaxima:");ToAliDebug(1,fRawCluster.Print());  
     
     if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){
       FitCoG(); //serialization of resolved clusters will happen inside
@@ -101,13 +101,13 @@ void AliRICHClusterFinder::FindClusters(Int_t iChamber)
 
   delete fDigitMap;
   
-  if(GetDebug())Info("FindClusters","Stop.");  
+  AliDebug(1,"Stop.");  
 }//FindClusters()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
 {
 //Finds cerenkov-feedback-mip mixture for a given cluster
-  if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();}
+  AliDebug(1,"Start.");ToAliDebug(1,pCluster->Print());
   
   TObjArray *pDigits = pCluster->Digits();
   if(!pDigits) return; //??????????
@@ -124,13 +124,13 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
   }//loop on digits of a given cluster
   TMath::Sort(contribs.GetSize(),contribs.GetArray(),pindex);
   for(Int_t iDigN=0;iDigN<3*pCluster->Size()-1;iDigN++) {//loop on digits to sort tids
-    if(GetDebug()) Info("FindClusterContribs","%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN);
+    AliDebug(1,Form("%4i for digit n. %4i",contribs[pindex[iDigN]],iDigN));
     if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
       TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]);
       Int_t code   = particle->GetPdgCode();
       Double_t charge = 0;
       if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
-      if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge);
+      AliDebug(1,Form(" charge of particle %f",charge));
 
       if(code==50000050) iNckovs++;
       if(code==50000051) iNfeeds++;
@@ -144,7 +144,7 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
      Int_t code   = particle->GetPdgCode();
      Double_t charge = 0;
      if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
-     if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge);
+     AliDebug(1,Form(" charge of particle %f",charge));
      if(code==50000050) iNckovs++;
      if(code==50000051) iNfeeds++;
      if(charge!=0) iNmips++;
@@ -153,13 +153,14 @@ void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
   pCluster->CFM(iNckovs,iNfeeds,iNmips);
 //  
   delete [] pindex; 
-  if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");}
+  ToAliDebug(1,pCluster->Print());
+  AliDebug(1,"Stop.");
 }//FindClusterContribs()
 //__________________________________________________________________________________________________
 void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
 {
 //Builds the raw cluster (before deconvolution). Starts from the first pad (i,j) then calls itself recursevly  for all neighbours.
-  if(GetDebug()) Info("FormRawCluster","Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q());
+  AliDebug(1,Form("Start with digit(%i,%i) Q=%f",i,j,((AliRICHdigit*)fDigitMap->GetHit(i,j))->Q()));
   
   fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
   fDigitMap->FlagHit(i,j);//flag this pad as taken  
@@ -204,30 +205,30 @@ void AliRICHClusterFinder::FindLocalMaxima()
 void AliRICHClusterFinder::WriteRawCluster()
 {
 //Add the current raw cluster to the list of clusters
-  if(GetDebug()) Info("WriteRawCluster","Start.");
+  AliDebug(1,"Start.");
   
   FindClusterContribs(&fRawCluster);  
   R()->AddCluster(fRawCluster);
   
-  if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");}  
+  ToAliDebug(1,fRawCluster.Print()); AliDebug(1,"Stop."); 
 }//WriteRawCluster()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::WriteResolvedCluster()
 {
 //Add the current resolved cluster to the list of clusters
-  if(GetDebug()) Info("WriteResolvedCluster","Start.");
+  AliDebug(1,"Start.");
   
   FindClusterContribs(&fResolvedCluster);  
   R()->AddCluster(fResolvedCluster);
   
-  if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");}  
+  ToAliDebug(1,fResolvedCluster.Print()); AliDebug(1,"Stop.");  
 }//WriteResolvedCluster()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FitCoG()
 {
 //Fits cluster of size  by the corresponding number of Mathieson shapes.
 //This methode is only invoked in case everything is ok to start deconvolution  
-  if(GetDebug()) {Info("FitCoG","Start with:"); fRawCluster.Print();}
+  AliDebug(1,"Start with:"); ToAliDebug(1,fRawCluster.Print());
   
   Double_t arglist;
   Int_t ierflag = 0;
@@ -293,7 +294,7 @@ void AliRICHClusterFinder::FitCoG()
     fResolvedCluster.Fill(&fRawCluster,xCoG[i],yCoG[i],qfracCoG[i],fLocalC[i]);
     WriteResolvedCluster();
   }
-  if(GetDebug()) Info("CoG","Stop.");
+  AliDebug(1,"Stop.");
 }//FitCoG()
 //__________________________________________________________________________________________________
 void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, Int_t )
index 656872f..54a253f 100644 (file)
@@ -80,7 +80,7 @@ void AliRICHDisplFast::Exec(Option_t *)
           AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j);
           if(pHit->C()==iChamber){
             TVector3 hitGlobX3= pHit->OutX3();
-            TVector2 hitLocX2 = pRich->C(iChamber)->Glob2Loc(hitGlobX3);
+            TVector2 hitLocX2 = pRich->C(iChamber)->Mrs2Pc(hitGlobX3);
             pHitsH2->Fill(hitLocX2.X(),hitLocX2.Y(),200);
           }//if
         }//hits loop         
@@ -95,7 +95,7 @@ void AliRICHDisplFast::Exec(Option_t *)
       if(!isClusters) {l.SetTextColor(kBlue) ;l.DrawLatex(0.8,0.01,"No CLUSTERS");}
       Display->Update();
       Display->Modified();
-      getchar();             
+      gPad->WaitPrimitive();
 //deals with digits      
       if(isDigits){
         pRich->GetLoader()->TreeD()->GetEntry(0);
@@ -108,7 +108,7 @@ void AliRICHDisplFast::Exec(Option_t *)
         pDigitsH2->Draw("same");
         Display->Update();
         Display->Modified();       
-        getchar();
+        gPad->WaitPrimitive();
       }//if(isDigits)      
 //deals with clusters      
       if(isClusters){
@@ -121,7 +121,7 @@ void AliRICHDisplFast::Exec(Option_t *)
         pClustersH2->Draw("same");
         Display->Update();
         Display->Modified();
-        getchar();
+        gPad->WaitPrimitive();
       }//if(isClusters)
     }//chambers loop
     delete [] Hits;
diff --git a/RICH/AliRICHHelix.cxx b/RICH/AliRICHHelix.cxx
new file mode 100644 (file)
index 0000000..c5ad6d6
--- /dev/null
@@ -0,0 +1,12 @@
+#include "AliRICHHelix.h"
+#include <AliLog.h>
+
+ClassImp(AliRICHHelix)
+
+void  AliRICHHelix::Print(Option_t *opt) const
+{
+// Debug printout
+  fX0.Print(opt);fP0.Print(opt);
+  AliInfo("Point of interest:");
+  fX.Print(opt); fP.Print(opt);
+}//Print()
diff --git a/RICH/AliRICHHelix.h b/RICH/AliRICHHelix.h
new file mode 100644 (file)
index 0000000..70cee59
--- /dev/null
@@ -0,0 +1,99 @@
+#ifndef AliRICHHelix_h
+#define AliRICHHelix_h
+
+#include <TObject.h>
+#include <TVector3.h>
+#include "AliRICHParam.h"
+#include "AliRICHChamber.h"
+
+class AliRICHHelix: public TObject
+{
+public:
+           AliRICHHelix()                                              {;}
+           AliRICHHelix(TVector3 x0,TVector3 p0,Int_t q,Double_t b=0.4)  {fX0=x0;fP0=p0;fQ=q;fBz=b;} //takes position and momentum at parametrised point
+  virtual ~AliRICHHelix()                                              {;}        
+  inline void   Propagate(Double_t lenght);
+         Bool_t Intersection(TVector3 plane)          {return Intersection(plane,plane.Unit());} // special plane given by point only
+  inline Bool_t Intersection(TVector3 planePoint,TVector3 planeNorm); //intersection with plane given by point and normal vector
+  inline Int_t  RichIntersect(AliRICHParam *pParam);
+  TVector3 X()const {return fX;}
+  TVector3 P()const {return fP;}
+  TVector3 X0()const {return fX0;}
+  TVector3 P0()const {return fP0;}
+  TVector2 PosPc() const {return fPosPc;}           //returns position of intersection with PC
+  TVector2 PosRad()const {return fPosRad;}          //returns position of intersection with radiator 
+  TVector3 Ploc()const   {return fPloc;}            //returns momentum at the position of intersection with radiator 
+  
+         void   Print(Option_t *sOption)const; //virtual interface from TObject
+protected:
+  TVector3 fX0;            //helix position in parametrised point, cm    in MRS
+  TVector3 fP0;            //helix momentum in parametrised point, GeV/c in MRS
+  TVector3 fX;             //helix position in point of interest,  cm    in MRS
+  TVector3 fP;             //helix momentum in point of interest,  GeV/c in MRS
+  Double_t fLength;        //helix length in point of interest    
+  Int_t    fQ;             //sign of track charge (value not provided by current ESD)  
+  Double_t fBz;            //magnetic field along z value in Tesla under assumption of uniformity 
+  TVector2 fPosPc;         //position on PC in local system
+  TVector2 fPosRad;        //position on radiator in local system 
+  TVector3 fPloc;          //momentum in local system
+  ClassDef(AliRICHHelix,0) //General helix
+};//class AliRICHHelix
+//__________________________________________________________________________________________________
+void AliRICHHelix::Propagate(Double_t length)
+{
+// Propogates the helix to the position of interest defined by helix length s  
+// Assumes uniform magnetic field along z direction.  
+  const Double_t c = 0.00299792458;//this value provides that coordinates are in cm momentum in GeV/c
+  Double_t a = -c*fBz*fQ;
+  Double_t rho = a/fP0.Mag();
+  fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*length)/a-fP0.Y()*(1-TMath::Cos(rho*length))/a  );
+  fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*length)/a+fP0.X()*(1-TMath::Cos(rho*length))/a  ); 
+  fX.SetZ( fX0.Z()+fP0.Z()*length/fP0.Mag()                                               );
+  fP.SetX( fP0.X()*TMath::Cos(rho*length)-fP0.Y()*TMath::Sin(rho*length)                  );
+  fP.SetY( fP0.Y()*TMath::Cos(rho*length)+fP0.X()*TMath::Sin(rho*length)                  );
+  fP.SetZ( fP0.Z()                                                                        );
+  fLength=length;
+}
+//__________________________________________________________________________________________________
+Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm)
+{
+// Finds point of intersection (if exists) of the helix to the plane given by point and normal vector.
+// Returns kTrue if helix intersects the plane, kFALSE otherwise.
+// Stores result in current helix fields fX and fP.   
+  
+  Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist;
+
+  while(TMath::Abs(dist)>0.0001){
+    Propagate(s);                        //calculates helix at the distance s from x0 ALONG the helix
+    dist=(fX-planePoint)*planeNorm;      //distance between current helix position and plane
+    if(TMath::Abs(dist) > TMath::Abs(distPrev)) { return kFALSE;}
+    distPrev=dist;
+    s-=dist;
+  }
+  return kTRUE;
+}
+//__________________________________________________________________________________________________
+Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam)
+{
+// Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection
+// On exit fPosRad contain position of intersection in Local System with radiator     
+//         fPosPc  contains the same for photocathode  
+  for(Int_t iChamberN=1;iChamberN<=kNchambers;iChamberN++){//chamber loop
+    if(Intersection(pParam->C(iChamberN)->Rad())){//there is intersection with radiator plane
+      fPosRad=pParam->C(iChamberN)->Mrs2Rad(fX);//position on radiator plane
+      if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone)
+        if(Intersection(pParam->C(iChamberN)->Pc())){//there is intersection with photocathode
+          fPosPc=pParam->C(iChamberN)->Mrs2Pc(fX);//position on photcathode plane
+          if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone)
+            fPloc=pParam->C(iChamberN)->PMrs2Loc(fP);
+            return iChamberN;
+          }//if inside PC
+        }//if for PC
+      }//if inside radiator
+    }//if for radiator       
+  }//chamber loop
+  return 0;
+}
+//__________________________________________________________________________________________________    
+#endif
index ad1679f..0271244 100644 (file)
@@ -28,15 +28,20 @@ Float_t  AliRICHParam::fgSigmaThSpread        =0.035; //
 //__________________________________________________________________________________________________
 void AliRICHParam::Print(Option_t*)
 {
-  ::Info("","Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec());
-  fpChambers->Print();
-}
+  AliInfo(Form("Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec()));
+  ToAliInfo(fpChambers->Print());
+}//Print()
 //__________________________________________________________________________________________________
 void AliRICHParam::CreateChambers()
 {
 //Create all RICH Chambers on each call. Previous chambers deleted.
   if(fpChambers) delete fpChambers;
-  fpChambers=new TObjArray(kNchambers);
+  if(IsRadioSrc()){ 
+    fpChambers=new TObjArray(1);//test beam configuration 1 chamber
+    fpChambers->AddAt(new AliRICHChamber(0),0);  
+  }else{ 
+    fpChambers=new TObjArray(kNchambers);//normal configuration 7 chambers
+    for(int iChamberN=0;iChamberN<kNchambers;iChamberN++)  fpChambers->AddAt(new AliRICHChamber(iChamberN+1),iChamberN);  
+  }
   fpChambers->SetOwner();
-  for(int iChamberN=0;iChamberN<kNchambers;iChamberN++)  fpChambers->AddAt(new AliRICHChamber(iChamberN+1),iChamberN);  
-}//void AliRICH::CreateChambers()
+}//CreateChambers()
index 4557e9e..2a664c1 100644 (file)
@@ -5,19 +5,24 @@
 #include <TMath.h>
 #include <TObjArray.h>
 #include <TObject.h>
+#include <TMath.h>
 #include <TRandom.h>
 #include <TVector.h>
 #include <TVector2.h>
 #include <TVector3.h>
+#include <TRandom.h>
+#include <TError.h>
+#include <TObjArray.h>
+#include <AliLog.h>
+#include <TClass.h>
+
 
-static const int kNCH=7;           //number of RICH chambers ???
 static const int kNchambers=7;     //number of RICH chambers 
 static const int kNpadsX = 160;    //number of pads along X in single chamber
 static const int kNpadsY = 144;    //number of pads along Y in single chamber
 static const int kBad=-101;        //useful static const to mark initial (uninitalised) values
 static const int kNsectors=6;      //number of sectors per chamber
 
-static const int kadc_satm  =  4096;  //dynamic range (10 bits)
 static const int kCerenkov=50000050;  //??? go to something more general like TPDGCode
 static const int kFeedback=50000051;  //??? go to something more general like TPDGCode
 
@@ -29,7 +34,7 @@ public:
                   AliRICHParam():TObject(),fpChambers(0)  {CreateChambers();}
   virtual        ~AliRICHParam()                          {delete fpChambers;}
          void     CreateChambers();
-         AliRICHChamber* C(Int_t i)          {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);}      //returns pointer to chamber i
+  AliRICHChamber* C(Int_t i)                 {return (AliRICHChamber*)fpChambers->UncheckedAt(i-1);}      //returns pointer to chamber i
   static Int_t    NpadsX()                   {return kNpadsX;}                           //pads along X in chamber
   static Int_t    NpadsY()                   {return kNpadsY;}                           //pads along Y in chamber
   static Int_t    NpadsXsec()                {return NpadsX()/2;}                        //pads along X in sector
@@ -37,41 +42,26 @@ public:
   static Double_t DeadZone()                 {return 2.6;}                               //dead zone size in cm  
   static Double_t PadSizeX()                 {return 0.8;}                               //pad size x in cm 
   static Double_t PadSizeY()                 {return 0.84;}                              //pad size y in cm   
+  
   static Double_t SectorSizeX()              {return NpadsX()*PadSizeX()/2;}             //sector size x in cm
   static Double_t SectorSizeY()              {return NpadsY()*PadSizeY()/3;}             //sector size y in cm 
-  static Double_t PcSizeX()                  {return NpadsX()*PadSizeX()+DeadZone();}    //photocathode size x in cm
-  static Double_t PcSizeY()                  {return NpadsY()*PadSizeY()+2*DeadZone();}  //photocathode size y in cm 
-  static Double_t SizeX()                    {return 132.6;}
-  static Double_t SizeY()                    {return 26;}
-  static Double_t SizeZ()                    {return 136.7;}                             
-  static Double_t Offset()                   {return 490+1.267;}                         //distance from IP to center of chamber in cm 
-  static Double_t AngleYZ()                  {return 19.5*TMath::DegToRad();}            //angle between chambers in YZ plane, rad
-  static Double_t AngleXY()                  {return 20*TMath::DegToRad();}              //angle between chambers in XY plane, rad
-  static Double_t AngleRot()                 {return fgAngleRot*TMath::DegToRad();}      //RICH rotation around Z, rad
-  static Double_t FreonThickness()           {return 1.5;}   
-  static Double_t QuartzThickness()          {return 0.5;}   
-  
-  static Double_t GapProx()                  {return 8.0;}                               //cm between CsI PC and radiator quartz window
-  static Double_t GapColl()                  {return 7.0;}                               //cm between CsI PC and third wire grid (collection wires)     
-  static Double_t GapAnod()                  {return 0.204;}                             //cm between CsI PC and first wire grid (anod wires)     
-  static Double_t GapAmp()                   {return 0.445;}                             //cm between CsI PC and second wire grid (cathode wires)
+  static Double_t PcSizeX()                  {return NpadsX()*PadSizeX()+DeadZone();}    //PC size x, cm
+  static Double_t PcSizeY()                  {return NpadsY()*PadSizeY()+2*DeadZone();}  //PC size y, cm
+   
+  static Double_t Zfreon()                   {return 1.5;}                               //freon thinkness, cm
+  static Double_t Zwin()                     {return 0.5;}                               //radiator quartz window, cm   
+  static Double_t Pc2Win()                   {return 8.0;}                               //cm between CsI PC and radiator quartz window
+  static Double_t Pc2Coll()                  {return 7.0;}                               //cm between CsI PC and third wire grid (collection wires)     
+  static Double_t Pc2Anod()                  {return 0.204;}                             //cm between CsI PC and first wire grid (anod wires)     
+  static Double_t Pc2Cath()                  {return 0.445;}                             //cm between CsI PC and second wire grid (cathode wires)
+  static Double_t Freon2Pc()                 {return Zfreon()+Zwin()+Pc2Win();}          //cm between CsI PC and entrance to freon
   static Double_t PitchAnod()                {return PadSizeY()/2;}                      //cm between anode wires
   static Double_t PitchCath()                {return PadSizeY()/4;}                      //cm between cathode wires
-  static Double_t PitchColl()                {return 0.5;}                               //cm between collect wires
+  static Double_t PitchColl()                {return 0.5;}                               //cm between collection wires
   
-  static Double_t GapThickness()             {return 8.0;}      
-  static Double_t RadiatorToPads()           {return FreonThickness()+QuartzThickness()+GapThickness();}   
-  static Double_t AnodeCathodeGap()          {return 0.2;}                               //between CsI PC and first wire grid     
-  static Double_t QuartzLength()             {return 133;}   
-  static Double_t QuartzWidth()              {return 127.9;}
-  static Double_t OuterFreonLength()         {return 133;}   
-  static Double_t OuterFreonWidth()          {return 41.3;}   
-  static Double_t InnerFreonLength()         {return 133;}   
-  static Double_t InnerFreonWidth()          {return 41.3;}   
   static Double_t IonisationPotential()      {return 26.0e-9;}                            //for CH4 in GeV taken from ????
   static TVector2 MathiesonDelta()           {return TVector2(5*0.18,5*0.18);}            //area of 5 sigmas of Mathieson distribution (cm)
   static Int_t    MaxQdc()                   {return 4095;}                               //QDC number of channels          
-  static Double_t AlphaFeedback(Int_t )      {return 0.030;}                              //determines number of feedback photons
     
   static Bool_t   IsResolveClusters()         {return fgIsResolveClusters;}  //go after resolved clusters?
   static Bool_t   IsWireSag()                 {return fgIsWireSag;}          //take wire sagita in account?
@@ -89,6 +79,10 @@ public:
   static void     SetWireSag(Bool_t status)   {fgIsWireSag=status;}  
   static void     SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;}  
   static void     SetAngleRot(Double_t rot)   {fgAngleRot =rot;}
+  static Double_t IndOfRefC6F14(Double_t eV)  {return eV*0.0172+1.177;}          // eV = photon energy in eV
+  static Double_t IndOfRefSiO2(Double_t eV)   {Double_t e1=10.666,e2=18.125,f1=46.411,f2= 228.71;
+                                     return TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(eV,2))+f2/(e2*e2-TMath::Power(eV,2)));}//TDR p.35
+  static Double_t IndOfRefCH4()               {return 1.000444;}
 
   inline static TVector  Loc2Area(TVector2 x2);                                                    //return area affected by hit x2
   inline static TVector  Loc2Pad(TVector2 x2);                                                     //return pad containing given position
@@ -107,10 +101,11 @@ public:
          static Float_t  SigmaThMean()                 {return fgSigmaThMean;}                     //QDC electronic noise mean
          static Float_t  SigmaThSpread()               {return fgSigmaThSpread;}                   //QDC electronic noise width
                 void     Print(const Option_t *opt="");                                            //virtual
-  inline static void     PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p);                
                 
   inline static Int_t    Loc2Sec(TVector2 &x2);             //return sector, x2->Sector RS
-  inline static Int_t    Pad2Sec(TVector pad);              //return sector
+  inline static Int_t    Pad2Sec(const TVector &pad);              //return sector
+         static Bool_t   IsAccepted(const TVector2 &x2) {return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) 
+? kTRUE:kFALSE;}
 protected:
          TObjArray *fpChambers;                             //list of chambers    
   static Bool_t     fgIsWireSag;                            //wire sagitta ON/OFF flag
@@ -156,12 +151,12 @@ Int_t AliRICHParam::Loc2Sec(TVector2 &v2)
   Double_t x=v2.X(),y=v2.Y();  
   if     (v2.X() >= x0 && v2.X() <= x1 )  {sector=1;}
   else if(v2.X() >= x2 && v2.X() <= x3 )  {sector=2; x=v2.X()-x2;}
-  else                                    {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in X",v2.X(),v2.Y());return kBad;}
+  else                                    {return kBad;}
   
   if     (v2.Y() >= y0 && v2.Y() <= y1 )  {}                                  //sectors 1 or 2 
   else if(v2.Y() >= y2 && v2.Y() <= y3 )  {sector+=2; y=v2.Y()-y2;}           //sectors 3 or 4
   else if(v2.Y() >= y4 && v2.Y() <= y5 )  {sector+=4; y=v2.Y()-y4;}           //sectors 5 or 6
-  else                                    {sector=kBad; ::Error("Loc2Sec","Position %6.2f,%6.2f is out of chamber in Y",v2.X(),v2.Y());return kBad;}
+  else                                    {return kBad;}
   v2.Set(x,y);
   return sector;
 }//Loc2Sec(Double_t x, Double_t y)
@@ -188,27 +183,27 @@ TVector AliRICHParam::Loc2Pad(TVector2 x2)
   return pad;
 }
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Pad2Sec(TVector pad)
+Int_t AliRICHParam::Pad2Sec(const TVector &pad)
 {
 // Determines sector containing the given pad.
   Int_t sector=kBad;      
   if     (pad[0] >= 1           && pad[0] <=   NpadsXsec() )    {sector=1;}
   else if(pad[0] >  NpadsXsec() && pad[0] <=   NpadsX()    )    {sector=2;} 
-  else                                                         ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+  else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
     
   if     (pad[1] >= 1             && pad[1] <=   NpadsYsec() )    {}
   else if(pad[1] >  NpadsYsec()   && pad[1] <= 2*NpadsYsec() )    {sector+=2;}
   else if(pad[1] >  2*NpadsYsec() && pad[1] <=   NpadsY()    )    {sector+=4;}
-  else                                                         ::Error("Pad2Sec","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+  else                                                         AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
 
   return sector;
 }//Pad2Sec()
 //__________________________________________________________________________________________________
 TVector2 AliRICHParam::Pad2Loc(TVector pad)
 {
-// Returns position of the center of the given pad in local system of the chamber    
+// Returns position of the center of the given pad in local system of the chamber (cm)    
 // y ^  5 6
-//   |  3 4        chamber structure
+//   |  3 4        sector numbers
 //   |  1 2
 //    -------> x  
   Double_t x=kBad,y=kBad;
@@ -217,7 +212,7 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad)
   else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6
     x=(pad[0]-0.5)*PadSizeX()+DeadZone();
   else
-    ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+    AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
   
   if(pad[1] > 0 && pad[1] <= NpadsYsec())//it's 1 or 2
     y=(pad[1]-0.5)*PadSizeY();
@@ -226,7 +221,7 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad)
   else if(pad[1] > 2*NpadsYsec() && pad[1]<= NpadsY())//it's 5 or 6
     y=(pad[1]-0.5)*PadSizeY()+2*DeadZone();
   else
-    ::Error("Pad2Loc","Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]);
+    AliDebugClass(1,Form("Wrong pad (%3.0f,%3.0f)",pad[0],pad[1]));
     
   return TVector2(x,y);
 }
@@ -234,7 +229,8 @@ TVector2 AliRICHParam::Pad2Loc(TVector pad)
 Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
 {
 // Returns % of gain variation due to wire sagita.
-// All curves are parametrized as per sector basis, so y must be apriory transformed to the Sector RS.    
+// All curves are parametrized as per sector basis, so x must be apriory transformed to the Sector RS.    
+// Here x is a distance along wires.  
   x-=SectorSizeX()/2;
   if(x>SectorSizeX()) x-=SectorSizeX(); 
   switch(HV(sector)){
@@ -250,8 +246,7 @@ Int_t AliRICHParam::TotQdc(TVector2 x2,Double_t eloss)
 {
 // Calculates the total charge produced by the eloss in point x2 (Chamber RS).
 // Returns this change parametrised in QDC channels, or 0 if the hit in the dead zone.
-// eloss=0 means photons which provided for only 1 electron
-// eloss > 0 for Mip
+// eloss=0 means photon which produces 1 electron only eloss > 0 for Mip
   if(Loc2Sec(x2)==kBad) return 0; //hit in the dead zone     
   Int_t iNelectrons=Int_t(eloss/IonisationPotential()); if(iNelectrons==0) iNelectrons=1;
   Double_t qdc=0;
@@ -264,13 +259,13 @@ Double_t AliRICHParam::FracQdc(TVector2 x2,TVector pad)
 // Calculates the charge fraction induced to given pad by the hit from the given point.
 // Integrated Mathieson distribution is used.  
   TVector2 center2=Pad2Loc(pad);//gives center of requested pad
-  Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2)  /AnodeCathodeGap();
-  Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2)  /AnodeCathodeGap();
-  Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2)  /AnodeCathodeGap();
-  Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2)  /AnodeCathodeGap();
-  
-  if(Loc2Sec(x2)!=Pad2Sec(pad)) return 0;//requested pad does not belong to the sector of the given hit position  
-  else                          return Mathieson(normXmin, normYmin, normXmax, normYmax);
+  Double_t normXmin=(x2.X()-center2.X()-PadSizeX()/2)  /Pc2Cath();//parametrise for Mathienson
+  Double_t normXmax=(x2.X()-center2.X()+PadSizeX()/2)  /Pc2Cath();
+  Double_t normYmin=(x2.Y()-center2.Y()-PadSizeY()/2)  /Pc2Cath();
+  Double_t normYmax=(x2.Y()-center2.Y()+PadSizeY()/2)  /Pc2Cath();
+//requested pad might not belong to the sector of the given hit position, hence the check:
+  return (Loc2Sec(x2)!=Pad2Sec(pad)) ? 0:Mathieson(normXmin, normYmin, normXmax, normYmax);
 }
 //__________________________________________________________________________________________________
 Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax)
@@ -306,21 +301,4 @@ Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q)
 // Checks if the current q is over threshold and FEE will save this value to data concentrator.
   return (q>NsigmaTh()*(SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread()));
 }
-//__________________________________________________________________________________________________
-void AliRICHParam::PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p)
-{
-// Propogates the helix given by (x0,p0) in MRS to the position of interest defined by helix length s  
-  const Double_t c = 0.00299792458;
-  const Double_t Bz = 0.5;       //field in Tesla
-  const Double_t q = 1;          //charge in electron units
-  Double_t a = -c*Bz*q;
-  Double_t rho = a/p0.Mag();
-  p->SetX(p0.X()*TMath::Cos(rho*s)-p0.Y()*TMath::Sin(rho*s));
-  p->SetY(p0.Y()*TMath::Cos(rho*s)+p0.X()*TMath::Sin(rho*s));
-  p->SetZ(p0.Z());
-  x->SetX(x0.X()+p0.X()*TMath::Sin(rho*s)/a-p0.Y()*(1-TMath::Cos(rho*s))/a);
-  x->SetY(x0.Y()+p0.Y()*TMath::Sin(rho*s)/a+p0.X()*(1-TMath::Cos(rho*s))/a);
-  x->SetZ(x0.Z()+p0.Z()*s/p->Mag());
-}
 #endif //AliRICHParam_h
index 16e4bc2..5ceada6 100644 (file)
 //////////////////////////////////////////////////////////////////////////
 
 #include <Riostream.h>
-#include <TCanvas.h>
-#include <TFile.h>
-#include <TH2.h>
 #include <TMath.h>
-#include <TMath.h>
-#include <TMinuit.h>
-#include <TNtuple.h>
-#include <TParticle.h>
-#include <TRandom.h>
 #include <TRotation.h>
 #include <TVector3.h>
 
-#include "AliLoader.h"
 #include "AliRICH.h"
-#include "AliRICHChamber.h"
 #include "AliRICHParam.h"
 #include "AliRICHRecon.h"
-#include "AliRun.h"
-#include "AliStack.h"
+#include "AliRICHHelix.h"
+#include <AliLog.h>
 
 #define NPointsOfRing 201
 
-TMinuit *gAliRICHminuit ;
-
-void fcnrecon(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
 //__________________________________________________________________________________________________
-AliRICHRecon::AliRICHRecon(const char*name, const char*title)
-             :TTask(name,title)
+AliRICHRecon::AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId)
+             :TTask("RichRec","RichPat")
 {
-  // main ctor
+// main ctor
+  SetFreonScaleFactor(1);
+  fIsWEIGHT = kFALSE;
   fThetaBin=750; fThetaMin = 0.0; fThetaMax = 0.75; 
-  fDTheta       = 0.001;   fWindowWidth  = 0.060;       
-  fNpadX         = AliRICHParam::NpadsY();
-  fNpadY         = AliRICHParam::NpadsX();
-  fPadSizeX      = AliRICHParam::PadSizeY();
-  fPadSizeY      = AliRICHParam::PadSizeX();
-  fRadiatorWidth = AliRICHParam::FreonThickness();
-  fQuartzWidth   = AliRICHParam::QuartzThickness();
-  fGapWidth      = AliRICHParam::RadiatorToPads();
-  fXmin         = -AliRICHParam::PcSizeY()/2.;
-  fXmax         =  AliRICHParam::PcSizeY()/2.;
-  fYmin         = -AliRICHParam::PcSizeX()/2.;
-  fYmax         =  AliRICHParam::PcSizeX()/2.;  
-  fRich = (AliRICH*)gAlice->GetDetector("RICH");
-  fOutFile=new TFile("Anal.root","RECREATE","My Analysis histos"); 
-  if(fIsDISPLAY) fDisplay = new TCanvas("Display","RICH Display",0,0,1200,750);      
-  fNtuple=new TNtuple("hn","ntuple",
-"Run:Trig:VertZ:Pmod:Pt:Eta:TrackTheta:TrackPhi:TrackThetaFit:TrackPhiFit:Charge::NPhotons:NPhotonsFit:InRing:MassOfParticle:HoughArea:Multiplicity:TPCLastZ");
+  fDTheta       = 0.001;   fWindowWidth  = 0.060;
+  fRadiatorWidth = AliRICHParam::Zfreon();
+  fQuartzWidth   = AliRICHParam::Zwin();
+  fGapWidth      = AliRICHParam::Freon2Pc() - fRadiatorWidth - fQuartzWidth;
+  fXmin         = -AliRICHParam::PcSizeX()/2.;
+  fXmax         =  AliRICHParam::PcSizeX()/2.;
+  fYmin         = -AliRICHParam::PcSizeY()/2.;
+  fYmax         =  AliRICHParam::PcSizeY()/2.; 
+  SetTrackTheta(pHelix->Ploc().Theta());
+  SetTrackPhi(pHelix->Ploc().Phi());
+  SetMipIndex(iMipId);
+  SetShiftX(pHelix->PosRad().X());
+  SetShiftY(pHelix->PosRad().Y());
+  fpClusters = pClusters;
 }
 //__________________________________________________________________________________________________
-void AliRICHRecon::StartProcessEvent()
+Double_t AliRICHRecon::ThetaCerenkov()
 {
-  //start to process for pattern recognition
-  
-  Float_t trackThetaStored    = 0;
-  Float_t trackPhiStored      = 0;
-  Float_t thetaCerenkovStored = 0;
-  Int_t houghPhotonsStored    = 0;
-  
-  SetFreonScaleFactor(0.994);
-
-  if(fIsDISPLAY) 
-    {
-      DrawEvent(0);
-//      Waiting();
-    }
-
-  AliLoader * richLoader = Rich()->GetLoader();
-  AliRunLoader * runLoader = richLoader->GetRunLoader();
-
-  if (richLoader->TreeH() == 0x0) richLoader->LoadHits();
-  if (richLoader->TreeR() == 0x0) richLoader->LoadRecPoints();
-  if (richLoader->TreeD() == 0x0) richLoader->LoadDigits();
-  if (runLoader->TreeE() == 0x0)  runLoader->LoadHeader();
-  if (runLoader->TreeK() == 0x0)  runLoader->LoadKinematics();    
-  
-  richLoader->TreeR()->GetEntry(0);
-
-  Float_t clusX[7][500],clusY[7][500];
-  Int_t clusQ[7][500],clusMul[7][500];    
-  Int_t nClusters[7];
-    
-  for (Int_t ich=0;ich<7;ich++) {
-    nClusters[ich] = Rich()->Clusters(ich+1)->GetEntries();    
-    for(Int_t k=0;k<nClusters[ich];k++) {
-      AliRICHcluster *pCluster = (AliRICHcluster *)Rich()->Clusters(ich+1)->At(k);
-      clusX[ich][k] = pCluster->X();
-      clusY[ich][k] = pCluster->Y();
-      clusQ[ich][k] = pCluster->Q();
-      clusMul[ich][k] = pCluster->Size();
-      //        pCluster->Print();
-    }
-  }
-  
-  Int_t nPrimaries = (Int_t)richLoader->TreeH()->GetEntries();
-  
-  cout << " N. primaries " << nPrimaries << endl;
-  
-  for(Int_t i=0;i<nPrimaries;i++){
-    
-    richLoader->TreeH()->GetEntry(i);
-    
-    //      Rich()->Hits()->Print();
-    Int_t iPrim = 0;
-    
-    AliRICHhit* pHit=0;
-    
-    for(Int_t j=0;j<Rich()->Hits()->GetEntries();j++) {
-      
-      pHit = (AliRICHhit*)Rich()->Hits()->At(j);
-      if(pHit->GetTrack() < nPrimaries) break;
-      iPrim++;
-    }
-    
-    cout << " iPrim " << iPrim << " pHit " << pHit << endl;
-    
-    if (!pHit) return;
-    
-    //      pHit->Print();
-    
-    TParticle *pParticle = runLoader->Stack()->Particle(pHit->GetTrack());
-    Float_t pmod     = pParticle->P();
-    Float_t pt       = pParticle->Pt();
-    Float_t trackEta = pParticle->Eta();
-    Int_t q          = (Int_t)TMath::Sign(1.,pParticle->GetPDG()->Charge());        
-    
-    //      pParticle->Print();
-    
-    cout << " pmod " << pmod << " pt " << pt << " Eta " << trackEta << " charge " << q << endl;
-    
-    SetTrackMomentum(pmod); 
-    SetTrackPt(pt);
-    SetTrackEta(trackEta);
-    SetTrackCharge(q);
-    
-    TVector3 pLocal(0,0,0);//?????
-    
-    TVector2 primLocal =Rich()->C(pHit->C())->Glob2Loc(pHit->InX3());
-    
-    //      Float_t pmodFreo = pLocal.Mag();
-    Float_t trackTheta = pLocal.Theta();
-    Float_t trackPhi = pLocal.Phi();
-    
-    //      cout << " trackTheta " << trackTheta << " trackPhi " << trackPhi << endl;
-    
-    SetTrackTheta(trackTheta);
-    SetTrackPhi(trackPhi);
-    
-    Int_t maxInd = 0;
-    Float_t minDist =  999.;
-    
-    //      cout << " n Clusters " << nClusters[pHit->Chamber()-1] << " for chamber n. " << pHit->Chamber() << endl;
-    
-    for(Int_t j=0;j<nClusters[pHit->Chamber()-1];j++)
-      {
-       Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][j];
-       Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][j];
-       
-       
-       Float_t diff = sqrt(diffx*diffx + diffy*diffy);
-       
-       if(diff < minDist)
-         {
-           minDist = diff;
-           maxInd = j;
-         }
-       
-      }
-    
-    Float_t diffx = primLocal.X() - clusX[pHit->Chamber()-1][maxInd];
-    Float_t diffy = primLocal.Y() - clusY[pHit->Chamber()-1][maxInd];
-    
-    cout << " diffx " << diffx << " diffy " << diffy << endl;
-    
-    
-    SetMipIndex(maxInd);
-    SetTrackIndex(i);
-    
-    Float_t shiftX = 0;//primLocal.X()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.X(); ????? 
-    Float_t shiftY = 0;//primLocal.Y()/primLocal.Z()*(fRadiatorWidth+fQuartzWidth+fGapWidth) + primLocal.Y(); ?????
-    
-    SetShiftX(shiftX);
-    SetShiftY(shiftY);
-    
-    Float_t *pclusX = &clusX[pHit->Chamber()-1][0];
-    Float_t *pclusY = &clusY[pHit->Chamber()-1][0];
-    
-    SetCandidatePhotonX(pclusX);
-    SetCandidatePhotonY(pclusY);
-    SetCandidatePhotonsNumber(nClusters[pHit->Chamber()-1]);
-    
-    Int_t qch = clusQ[pHit->Chamber()-1][maxInd];
-    
-    
-    if(minDist < 3.0 && qch > 120 && maxInd !=0) 
-      {
-       
-       if(fIsBACKGROUND)
-         {
-           
-           Float_t xrndm = fXmin + (fXmax-fXmin)*gRandom->Rndm(280964);
-           Float_t yrndm = fYmin + (fYmax-fYmin)*gRandom->Rndm(280964);
-           SetShiftX(xrndm);
-           SetShiftY(yrndm);
-           
-         }
-       
-       PatRec();
-       
-       trackThetaStored = GetTrackTheta();
-       trackPhiStored = GetTrackPhi();
-       thetaCerenkovStored = GetThetaCerenkov();
-       houghPhotonsStored = GetHoughPhotons();
-       
-       Int_t diffNPhotons = 999;
-       Int_t nsteps = 0;
-       Float_t diffTrackTheta = 999.;
-       Float_t diffTrackPhi   = 999.;
-       
-       while(fIsMINIMIZER && GetHoughPhotons() > 2 
-             && diffNPhotons !=0 
-             && diffTrackTheta > 0.0001
-             && nsteps < 10)
-         {
-           
-           Int_t   houghPhotonsBefore  = GetHoughPhotons();
-           
-           Float_t trackThetaBefore = GetTrackTheta();
-           Float_t trackPhiBefore   = GetTrackPhi();
-           
-           Minimization(); 
-           
-           PatRec();
-           
-           diffNPhotons = TMath::Abs(houghPhotonsBefore - GetHoughPhotons()); 
-           
-           Float_t trackThetaAfter = GetTrackTheta();
-           Float_t trackPhiAfter   = GetTrackPhi();
-           
-           diffTrackTheta = TMath::Abs(trackThetaAfter - trackThetaBefore);
-           diffTrackPhi   = TMath::Abs(trackPhiAfter - trackPhiBefore);
-           
-           if(fDebug)
-              cout << " houghPhotonsBefore " << houghPhotonsBefore
-                   << " GetHoughPhotons()  " << GetHoughPhotons();
-           
-           nsteps++;
-         }
-       
-       SetFittedThetaCerenkov(GetThetaCerenkov());
-       SetFittedHoughPhotons(GetHoughPhotons());
-       
-       SetTrackTheta(trackThetaStored);
-       SetTrackPhi(trackPhiStored);
-       SetThetaCerenkov(thetaCerenkovStored);
-       SetHoughPhotons(houghPhotonsStored);
-       
-       SetMinDist(minDist);
-       
-       FillHistograms();
-       
-       if(fIsDISPLAY) DrawEvent(1);
-       
-       Waiting();
-       
-      }
-  }
-  if(fIsDISPLAY) fDisplay->Print("display.ps");
-}//StartProcessEvent()
-//__________________________________________________________________________________________________
-void AliRICHRecon::EndProcessEvent()
-{
-  // function called at the end of the event loop
-
-  fOutFile->Write();
-  fOutFile->Close();                                                     
-}
-//__________________________________________________________________________________________________
-void AliRICHRecon::PatRec()
-{
-  //pattern recognition method based on Hough transform
-
-  
-  Float_t trackTheta = GetTrackTheta();
-  Float_t trackPhi   = GetTrackPhi();
-  Float_t pmod       = GetTrackMomentum();
-  Int_t iMipIndex   = GetMipIndex();
+// Pattern recognition method based on Hough transform
+// Return theta Cerenkov for a given track and list of clusters which are set in ctor  
 
+  if(fpClusters->GetEntries()==0) return kBad;//no clusters at all for a given track
   Bool_t kPatRec = kFALSE;  
 
   Int_t candidatePhotons = 0;
 
-  Float_t shiftX = GetShiftX();
-  Float_t shiftY = GetShiftY();
-
-  Float_t* candidatePhotonX = GetCandidatePhotonX();
-  Float_t* candidatePhotonY = GetCandidatePhotonY();
-
-  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
-
-  if(fDebug) cout << " n " << candidatePhotonsNumber << endl;
-
   SetThetaCerenkov(999.);
   SetHoughPhotons(0);
   SetHoughPhotonsNorm(0);
-  SetHoughRMS(999.);
-
-  for (Int_t j=0; j < candidatePhotonsNumber; j++)
-    {
-
-      SetPhotonIndex(j);
-
-      SetPhotonFlag(0);
-      SetPhotonEta(-999.);
-      SetPhotonWeight(0.);
-
-      if (j == iMipIndex) continue;
-
-        
-      if(candidatePhotonX[j] < -64.) continue; /* avoid artificial clusters from edge uesd by Yale.... */
-
-      Float_t xtoentr = candidatePhotonX[j] - shiftX;
-      Float_t ytoentr = candidatePhotonY[j] - shiftY;
-
-      //      Float_t chargehit = fHits_charge[j]; 
-      //      if(chargehit > 150) continue;
-
-      SetEntranceX(xtoentr);
-      SetEntranceY(ytoentr);
-
-      FindPhiPoint();
-
-      Int_t photonStatus = PhotonInBand();
-      if(fDebug)
-         {
-            cout << " Photon n. " << j << " Status " << photonStatus << " accepted " << endl;
-            cout << " CandidatePhotonX[j] " << candidatePhotonX[j] << " CandidatePhotonY[j] " << candidatePhotonY[j] << endl;
-         }
-    
-      if(photonStatus == 0) continue;
-
-      SetPhotonFlag(1);
-
-      FindThetaPhotonCerenkov();
-
-      Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
-
-      if(fDebug) cout << " theta photon " << thetaPhotonCerenkov << endl;
 
-      SetPhotonEta(thetaPhotonCerenkov);
-
-      candidatePhotons++;
-
-      
-    }
+  for (Int_t j=0; j < fpClusters->GetEntries(); j++){//clusters loop
+    SetPhotonIndex(j);
+    SetPhotonFlag(0);
+    SetPhotonEta(-999.);
+    SetPhotonWeight(0.);
+    if (j == GetMipIndex()) continue; // do not consider MIP cluster as a candidate photon
+    Float_t xtoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->X() - GetShiftX();
+    Float_t ytoentr = ((AliRICHcluster*)fpClusters->UncheckedAt(j))->Y() - GetShiftY();
+    SetEntranceX(xtoentr);
+    SetEntranceY(ytoentr);
+    FindPhiPoint();
+//      Int_t photonStatus = PhotonInBand();
+//      if(photonStatus == 0) continue;
+    SetPhotonFlag(1);
+    FindThetaPhotonCerenkov();
+    Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
+    AliDebug(1,Form("THETA CERENKOV ---> %i",thetaPhotonCerenkov));
+    SetPhotonEta(thetaPhotonCerenkov);
+    candidatePhotons++;
+  }//clusters loop
 
   if(candidatePhotons >= 1) kPatRec = kTRUE;
 
-  if(!kPatRec) return;
-    {
-       SetThetaCerenkov(999.);
-       SetHoughPhotons(0);
-    }
-  SetPhotonsNumber(candidatePhotonsNumber);
+  if(!kPatRec) return kBad;
+
+  SetPhotonsNumber(fpClusters->GetEntries());
 
   HoughResponse();
   
@@ -399,21 +110,23 @@ void AliRICHRecon::PatRec()
     {
       SetThetaCerenkov(999.);
       SetHoughPhotonsNorm(0.);
-      return;
+      return kBad;
     }
 
   if(fIsWEIGHT) FindWeightThetaCerenkov();
 
   Float_t thetaCerenkov = GetThetaCerenkov();
 
+  AliDebug(1,Form("Number of clusters accepted --->  %i",nPhotonHough));
+  
   SetThetaOfRing(thetaCerenkov);
-  FindAreaAndPortionOfRing();
+//  FindAreaAndPortionOfRing();
 
-  Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
-  SetHoughPhotonsNorm(nPhotonHoughNorm);
+//  Float_t nPhotonHoughNorm = ((Float_t)nPhotonHough)/GetPortionOfRing();
+//  SetHoughPhotonsNorm(nPhotonHoughNorm);
 
   // Calculate the area where the photon are accepted...
-
+/*
   Float_t thetaInternal = thetaCerenkov - 0.5*fWindowWidth; 
   SetThetaOfRing(thetaInternal);
   FindAreaAndPortionOfRing();
@@ -427,63 +140,11 @@ void AliRICHRecon::PatRec()
   Float_t houghArea = externalArea - internalArea;
 
   SetHoughArea(houghArea);
+*/
+  return GetThetaCerenkov();
 
-  if(fDebug)
-    {
-      cout << " ----- SUMMARY OF RECONSTRUCTION ----- " << endl; 
-      cout << " Rings found " << fNrings << " with thetac " << thetaCerenkov << endl;
-      
-      
-      cout << " Nphotons " << GetPhotonsNumber() 
-          << " Hough    " << nPhotonHough 
-          << " norm     " << nPhotonHoughNorm << endl;
-      
-      cout << " In PatRec:p " << pmod << " theta " << trackTheta << " phi " << trackPhi << endl;
-      cout << " ------------------------------------- " << endl; 
-    }
-
-  Int_t nPhotons = GetPhotonsNumber();
-
-  Float_t xmean = 0.;
-  Float_t x2mean = 0.;
-  Int_t nev = 0;
-
-  for (Int_t j=0; j < nPhotons;j++)
-    {
-      SetPhotonIndex(j);
-
-      Float_t eta = GetPhotonEta();
-
-      if(eta != -999.) 
-       {
-         if(GetPhotonFlag() == 2) 
-           {
-
-
-             xmean += eta;
-             x2mean += eta*eta;
-             nev++;
-           }
-       }
-    }
-
-  if(nev > 0)
-    {
-      xmean /=(Float_t)nev;
-      x2mean /=(Float_t)nev;
-    } else {
-      xmean = 0.;
-      x2mean = 0.;
-    }
-
-  Float_t vRMS = sqrt(x2mean - xmean*xmean);
-
-  SetHoughRMS(vRMS);
-
-  if(fDebug) cout << " RMS " << vRMS << endl;
-
-}
-
+}//ThetaCerenkov()
+//__________________________________________________________________________________________________
 void AliRICHRecon::FindEmissionPoint()
 {
   //estimate the emission point in radiator
@@ -501,6 +162,7 @@ void AliRICHRecon::FindEmissionPoint()
   Float_t emissionPoint = fRadiatorWidth + photonLenghtMin - photonLenghtMax;
 
   SetEmissionPoint(emissionPoint);
+  SetEmissionPoint(fRadiatorWidth/2); // tune the emission point
 }
 
 
@@ -509,7 +171,6 @@ Int_t AliRICHRecon::PhotonInBand()
   //search band fro photon candidates
 
   //  Float_t massOfParticle;
-  Float_t beta;
   Float_t nfreon;
 
   Float_t thetacer;
@@ -522,26 +183,16 @@ Int_t AliRICHRecon::PhotonInBand()
 
   Float_t phpad = GetPhiPoint();
 
-  //  Float_t pmod = GetTrackMomentum();
-  //  Float_t trackTheta = GetTrackTheta();
-  //  Float_t trackPhi = GetTrackPhi();
 
   // inner radius //
   SetPhotonEnergy(5.6);
   SetEmissionPoint(fRadiatorWidth -0.0001);
-  SetMassHypotesis(0.93828);
-
-  SetBetaOfParticle();
   SetFreonRefractiveIndex();
 
-  beta   = GetBetaOfParticle();
   nfreon = GetFreonRefractiveIndex();
-
-  thetacer = Cerenkovangle(nfreon,beta);
-
   thetacer = 0.;
 
-  if(fDebug) cout << " thetacer in photoninband min " << thetacer << endl;
+  AliDebug(1,Form("thetacer in photoninband min %f",thetacer));
 
   FindThetaAtQuartz(thetacer);
 
@@ -569,19 +220,15 @@ Int_t AliRICHRecon::PhotonInBand()
   SetPhotonEnergy(7.7);
   SetEmissionPoint(0.);
 //  SetMassHypotesis(0.139567);
-  SetMassHypotesis(0.);
-
-  SetBetaOfParticle();
   SetFreonRefractiveIndex();
 
-  beta   = GetBetaOfParticle();
   nfreon = GetFreonRefractiveIndex();
 
-  thetacer = Cerenkovangle(nfreon,beta);
+  thetacer = Cerenkovangle(nfreon,1);
 
   //  thetacer = 0.75;
 
-  if(fDebug) cout << " thetacer in photoninband max " << thetacer << endl;
+  AliDebug(1,Form("thetacer in photoninband max %f",thetacer));
 
   FindThetaAtQuartz(thetacer);
 
@@ -606,7 +253,7 @@ Int_t AliRICHRecon::PhotonInBand()
 
   Float_t padradius = sqrt(TMath::Power(xtoentr,2)+TMath::Power(ytoentr,2));
   
-  if(fDebug) printf(" rmin %f r %f rmax %f \n",innerRadius,padradius,outerRadius);
+  AliDebug(1,Form("rmin %f r %f rmax %f",innerRadius,padradius,outerRadius));
 
   if(padradius>=innerRadius && padradius<=outerRadius) return 1;
   return 0;
@@ -628,8 +275,6 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
 
   if(trackTheta == 0) {
 
-    if(fDebug) cout << " Theta sol unique " << thetaCerenkov << endl;  
-
     thetaAtQuartz = thetaCerenkov;
     SetThetaAtQuartz(thetaAtQuartz);
     return;
@@ -649,17 +294,7 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
 
   Double_t underSqrt = 1 + b*b - c*c;
 
-  if(fDebug)
-    {
-      cout << " trackTheta    " << trackTheta    << endl;
-      cout << " TrackPhi      " << trackPhi      << endl;
-      cout << " PhiPoint      " << phiPoint      << endl;
-      cout << " ThetaCerenkov " << thetaCerenkov << endl;
-      cout << " den b c " << den << " b " << b << " c " << c << endl;
-    }
-
   if(underSqrt < 0) {
-    if(fDebug) cout << " sqrt negative !!!!" << underSqrt << endl;
     SetThetaAtQuartz(999.);
     return;
   }
@@ -670,12 +305,11 @@ void AliRICHRecon::FindThetaAtQuartz(Float_t thetaCerenkov)
   Double_t thetaSol1 = 2*TMath::ATan(sol1);
   Double_t thetaSol2 = 2*TMath::ATan(sol2);
 
-  if(fDebug) cout << " Theta sol 1 " << thetaSol1 
-                 << " Theta sol 2 " << thetaSol2 << endl;  
-
   if(thetaSol1>0 && thetaSol1 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol1;
   if(thetaSol2>0 && thetaSol2 < TMath::Pi()) thetaAtQuartz = (Float_t)thetaSol2;
 
+//  AliDebug(1,Form(" Theta @ quartz window %f ",thetaAtQuartz));
+
   SetThetaAtQuartz(thetaAtQuartz);
 }
 
@@ -692,9 +326,6 @@ void AliRICHRecon::FindThetaPhotonCerenkov()
 
   const Float_t kTollerance = 0.05;
 
-  //  Float_t pmod = GetTrackMomentum();
-  //  Float_t trackTheta = GetTrackTheta();
-  //  Float_t trackPhi = GetTrackPhi();
 
   Float_t phiPoint = GetPhiPoint();
 
@@ -703,9 +334,9 @@ void AliRICHRecon::FindThetaPhotonCerenkov()
 
   Float_t xPoint = GetEntranceX();
   Float_t yPoint = GetEntranceY();
-  Float_t distPoint = sqrt(xPoint*xPoint + yPoint*yPoint);
+  Float_t distPoint = TMath::Sqrt(xPoint*xPoint + yPoint*yPoint);
 
-  if(fDebug) cout << " DistPoint " << distPoint << endl;
+  AliDebug(1,Form(" DistPoint %f ",distPoint));
 
   // Star minimization...
 
@@ -756,8 +387,7 @@ void AliRICHRecon::FindThetaPhotonCerenkov()
       radiusMean = FromEmissionToCathode();
     }
 
-  if(fDebug) cout << " r1 " << radiusMin << " rmean " 
-                 << radiusMean << " r2 " << radiusMax << endl;
+  AliDebug(1,Form(" r1 %f rmean %f r2 %f",radiusMin,radiusMean,radiusMax));
 
   while (TMath::Abs(radiusMean-distPoint) > kTollerance)
     {
@@ -784,12 +414,13 @@ void AliRICHRecon::FindThetaPhotonCerenkov()
 
       nIteration++;
       if(nIteration>=50) {
-       if(fDebug) printf(" max iterations in FindPhotonCerenkov\n");
+       AliDebug(1,Form(" max iterations in FindPhotonCerenkov ",nIteration));
        SetThetaPhotonCerenkov(999.);
        return;
       }
     }
 
+  AliDebug(1,Form(" distpoint %f radius %f ",distPoint,radiusMean));
   SetThetaPhotonCerenkov(thetaCerMean);
 
 }
@@ -811,9 +442,6 @@ void AliRICHRecon::FindAreaAndPortionOfRing()
   Float_t x0 = xemiss + shiftX;
   Float_t y0 = yemiss + shiftY;
 
-  //  Float_t pmod = GetTrackMomentum();
-  //  Float_t trackTheta = GetTrackTheta();
-  //  Float_t trackPhi = GetTrackPhi();
 
   SetPhotonEnergy(6.85);
   SetFreonRefractiveIndex();
@@ -1001,31 +629,6 @@ Int_t AliRICHRecon::CheckDetectorAcceptance() const
   return 999;
 }
 //__________________________________________________________________________________________________
-Float_t AliRICHRecon::PhotonPositionOnCathode()
-{ 
-  // find the photon position on the CsI
-  //  Float_t massOfParticle;
-  Float_t beta;
-  Float_t nfreon;
-
-
-  SetPhotonEnergy(6.85);
-  SetEmissionPoint(fRadiatorWidth/2.);
-  SetMassHypotesis(0.139567);
-
-  SetBetaOfParticle();
-  SetFreonRefractiveIndex();
-
-  beta   = GetBetaOfParticle();   
-  nfreon = GetFreonRefractiveIndex();
-
-
-  Float_t radius = FromEmissionToCathode();
-  if (radius == 999.) return 999.;
-
-  return 0;
-}
-
 void AliRICHRecon::FindPhotonAnglesInDRS()
 {
   // Setup the rotation matrix of the track...
@@ -1143,9 +746,9 @@ void AliRICHRecon::FindPhiPoint()
 
   Float_t argY = ytoentr - emissionPoint*tan(trackTheta)*sin(trackPhi);
   Float_t argX = xtoentr - emissionPoint*tan(trackTheta)*cos(trackPhi);
-  Float_t phipad = atan2(argY,argX); 
+  Float_t phi = atan2(argY,argX); 
 
-  SetPhiPoint(phipad);
+  SetPhiPoint(phi);
 
 }
 
@@ -1406,7 +1009,7 @@ void AliRICHRecon::FlagPhotons()
   Int_t nPhotonHough = 0;
 
   Float_t thetaCerenkov = GetThetaCerenkov();
-  if(fDebug) cout << " fThetaCerenkov " << thetaCerenkov << endl;
+  AliDebug(1,Form(" fThetaCerenkov ",thetaCerenkov));
 
   Float_t thetaDist= thetaCerenkov - fThetaMin;
   Int_t steps = (Int_t)(thetaDist / fDTheta);
@@ -1418,9 +1021,6 @@ void AliRICHRecon::FlagPhotons()
   tmin = tavg - 0.5*fWindowWidth;
   tmax = tavg + 0.5*fWindowWidth;
 
-  if(fDebug) cout << " tmin " << tmin << " tmax " << tmax << endl;
-  if(fDebug) cout << " thetac " << thetaCerenkov << endl;
-
   //  Int_t candidatePhotonsNumber = GetCandidatePhotonsNumber();
 
   Int_t nPhotons = GetPhotonsNumber();
@@ -1444,305 +1044,3 @@ void AliRICHRecon::FlagPhotons()
   SetHoughPhotons(nPhotonHough);
 }
 
-void AliRICHRecon::DrawEvent(Int_t flag) const
-{
-  // draw event with rings
-
-  flag=1; // dummy to be removed...
-}
-
-Float_t  AliRICHRecon::FindMassOfParticle()
-{
-  // find mass of the particle from theta cerenkov
-
-  Float_t pmod = GetTrackMomentum();
-
-  SetPhotonEnergy(6.85);
-  SetFreonRefractiveIndex();
-
-  Float_t thetaCerenkov = GetThetaCerenkov();
-  FindBetaFromTheta(thetaCerenkov);
-
-  Double_t beta = (Double_t)(GetBetaOfParticle());
-  Double_t den = 1. - beta*beta;
-  if(den<=0.) return 999.;
-
-  Double_t gamma = 1./TMath::Sqrt(den);
-
-  Float_t mass = pmod/(beta*(Float_t)gamma);
-
-  return mass;
-}
-
-
-void AliRICHRecon::FillHistograms()
-{
-  // fill histograms..
-
-  Float_t fittedTrackTheta, fittedTrackPhi;
-
-  Float_t thetaCerenkov    = GetThetaCerenkov();
-  if(thetaCerenkov == 999.) return;
-
-  Float_t vertZ = GetEventVertexZ();
-
-  Float_t trackTheta = GetTrackTheta();
-  Float_t trackPhi   = GetTrackPhi();
-  Float_t pmod       = GetTrackMomentum();
-  Float_t pt         = GetTrackPt();
-  Float_t trackEta   = GetTrackEta();
-  Int_t q            = GetTrackCharge();
-  Float_t tPCLastZ   = GetTrackTPCLastZ(); 
-  Float_t minDist    = GetMinDist(); 
-
-  fittedTrackTheta = GetFittedTrackTheta();
-  fittedTrackPhi   = GetFittedTrackPhi();
-  Int_t fittednPhotonHough = GetFittedHoughPhotons();
-  
-  if(fDebug)
-    {
-      cout << " p " << pmod  << " ThetaC " << thetaCerenkov 
-          << " rings " << fNrings << endl;
-    }
-
-  Int_t nPhotonHough     = GetHoughPhotons();
-//  Float_t nPhotonHoughNorm = GetHoughPhotonsNorm();
-  Float_t inRing = GetPortionOfRing();
-
-  Float_t massOfParticle = FindMassOfParticle();
-
-  Float_t houghArea = GetHoughArea();
-  Float_t multiplicity = GetEventMultiplicity();
-
-
-  Float_t var[20];
-
-  var[0] = 0; 
-  var[1] = 0;
-  var[2] = vertZ;
-  var[3] = pmod;
-  var[4] = pt;
-  var[5] = trackEta;
-  var[6] = trackTheta;
-  var[7] = trackPhi;
-  var[8] = fittedTrackTheta;
-  var[9] = fittedTrackPhi;
-  var[10] = q;
-  var[11] = thetaCerenkov;
-  var[12] = (Float_t)nPhotonHough;
-  var[13] = (Float_t)fittednPhotonHough;
-  var[14] = inRing;
-  var[15] = massOfParticle;
-  var[16] = houghArea;
-  var[17] = multiplicity;
-  var[18] = tPCLastZ;
-  var[19] = minDist;
-
-  fNtuple->Fill(var);
-
-
-  fittedTrackTheta = GetFittedTrackTheta();
-  fittedTrackPhi = GetFittedTrackPhi();
-
-
-
-  if(thetaCerenkov > 0.505 && thetaCerenkov < 0.605) {
-      SetPhotonEnergy(6.85);
-      SetFreonRefractiveIndex();
-  }
-
-  Int_t nPhotons = GetPhotonsNumber();
-
-  for (Int_t j=0; j < nPhotons;j++)
-    SetPhotonIndex(j);
-}//FillHistograms()
-//__________________________________________________________________________________________________
-void AliRICHRecon::Minimization()
-{
-  // minimization to find the best theta and phi of the track
-
-  Double_t arglist;
-  Int_t ierflag = 0;
-
-  static Double_t vstart[2];
-  static Double_t lower[2], upper[2];
-  static Double_t step[2]={0.001,0.001};
-
-  Double_t trackThetaNew,trackPhiNew;
-  TString chname;
-  Double_t eps, b1, b2;
-  Int_t ierflg;
-
-  gAliRICHminuit = new TMinuit(2);
-  gAliRICHminuit->SetObjectFit((TObject *)this);
-  gAliRICHminuit->SetFCN(fcnrecon);
-  gAliRICHminuit->mninit(5,10,7);
-
-  vstart[0] = (Double_t)GetTrackTheta();
-  vstart[1] = (Double_t)GetTrackPhi();
-
-  lower[0] = vstart[0] - 0.03;
-  if(lower[0] < 0) lower[0] = 0.;
-  upper[0] = vstart[0] + 0.03;
-  lower[1] = vstart[1] - 0.03;
-  upper[1] = vstart[1] + 0.03;
-
-
-  gAliRICHminuit->mnparm(0,"theta",vstart[0],step[0],lower[0],upper[0],ierflag);
-  gAliRICHminuit->mnparm(1," phi ",vstart[1],step[1],lower[1],upper[1],ierflag);
-
-  arglist = -1;
-
-  //  gAliRICHminuit->FixParameter(0);
-
-  gAliRICHminuit->SetPrintLevel(-1);
-//  gAliRICHminuit->mnexcm("SET PRI",&arglist, 1, ierflag);
-  gAliRICHminuit->mnexcm("SET NOGR",&arglist, 1, ierflag);
-  gAliRICHminuit->mnexcm("SET NOW",&arglist, 1, ierflag);
-  arglist = 1;
-  gAliRICHminuit->mnexcm("SET ERR", &arglist, 1,ierflg);
-  arglist = -1;
-
-  //  gAliRICHminuit->mnscan();
-
-//  gAliRICHminuit->mnexcm("SIMPLEX",&arglist, 0, ierflag);
-  gAliRICHminuit->mnexcm("MIGRAD",&arglist, 0, ierflag);
-  gAliRICHminuit->mnexcm("EXIT" ,&arglist, 0, ierflag);
-  
-  gAliRICHminuit->mnpout(0,chname, trackThetaNew, eps , b1, b2, ierflg);
-  gAliRICHminuit->mnpout(1,chname, trackPhiNew, eps , b1, b2, ierflg);
-
-  //values after the fit...
-  SetFittedTrackTheta((Float_t)trackThetaNew);
-  SetFittedTrackPhi((Float_t)trackPhiNew);
-
-  delete gAliRICHminuit;
-
-}
-
-void AliRICHRecon::EstimationOfTheta()
-{
-  // theta estimate
-
-  Int_t nPhotons = 0;
-
-  Float_t shiftX = GetShiftX();
-  Float_t shiftY = GetShiftY();
-
-  Float_t *candidatePhotonX = GetCandidatePhotonX();
-  Float_t *candidatePhotonY = GetCandidatePhotonY();
-
-  Int_t nPhotonsCandidates = GetCandidatePhotonsNumber();
-
-  //  cout << "MINIM: Nphotons " << nPhotonsCandidates << endl;
-
-  for (Int_t j=0; j < nPhotonsCandidates; j++)
-    {
-
-      SetPhotonIndex(j);
-
-      if(!GetPhotonFlag()) continue;
-
-      Float_t xtoentr = candidatePhotonX[j] - shiftX;
-      Float_t ytoentr = candidatePhotonY[j] - shiftY;
-
-      SetEntranceX(xtoentr);
-      SetEntranceY(ytoentr);
-
-      FindPhiPoint();
-
-      FindThetaPhotonCerenkov();
-
-      Float_t thetaPhotonCerenkov = GetThetaPhotonCerenkov();
-
-      //      cout << " ACCEPTED!!! " << thetaPhotonCerenkov << endl;
-
-      SetPhotonEta(thetaPhotonCerenkov);
-
-      nPhotons++;
-
-    }
-
-  Float_t xmean = 0.;
-  Float_t x2mean = 0.;
-  Int_t nev = 0;
-
-  for (Int_t j=0; j < nPhotonsCandidates;j++)
-    {
-      SetPhotonIndex(j);
-
-      Float_t eta = GetPhotonEta();
-
-      if(eta != -999.) 
-       {
-         if(GetPhotonFlag() == 2) 
-           {
-             xmean += eta;
-             x2mean += eta*eta;
-             nev++;
-           }
-       }
-    }
-
-  if(nev > 0)
-    {
-      xmean /=(Float_t)nev;
-      x2mean /=(Float_t)nev;
-    } else {
-      xmean = 0.;
-      x2mean = 0.;
-    }
-
-  Float_t vRMS = sqrt(x2mean - xmean*xmean);
-
-  //  cout << " RMS " << vRMS;
-
-  SetEstimationOfTheta(xmean);
-  SetEstimationOfThetaRMS(vRMS);
-}
-
-void fcnrecon(Int_t& /*npar*/, Double_t* /*gin*/, Double_t &f, Double_t *par, Int_t)
-{
-  // function to be minimized
-  AliRICHRecon *gMyRecon = (AliRICHRecon*)gAliRICHminuit->GetObjectFit();
-
-  Float_t p0 = (Float_t)par[0];
-  Float_t p1 = (Float_t)par[1];
-
-  gMyRecon->SetTrackTheta(p0);
-  gMyRecon->SetTrackPhi(p1);
-
-  gMyRecon->EstimationOfTheta();
-  Float_t vRMS = gMyRecon->GetEstimationOfThetaRMS();
-
-  Int_t houghPhotons = gMyRecon->GetHoughPhotons();
-
-
-  f = (Double_t)(1000*vRMS/(Float_t)houghPhotons);
-
-//   if(fDebug) cout << "   f   " << f
-//               << " theta " << par[0] << " phi " << par[1] 
-//                   << " HoughPhotons " << houghPhotons << endl;
-//   
-//   if(fDebug&&iflag == 3)
-//     {
-//             cout << " --- end convergence...summary --- " << endl;
-//             cout << " theta " << par[0] << endl;
-//             cout << "  phi  " << par[1] << endl;
-//     }
-}
-
-void AliRICHRecon::Waiting()
-{
-  // wait, wait....
-  if(!fIsDISPLAY) return;
-  cout << " Press any key to continue...";
-
-//  gSystem->ProcessEvents();
-  getchar(); 
-
-  cout << endl;
-
-  return;
-}
-
index 03132c5..419ce81 100644 (file)
 
 #include <TTask.h>
 
-class AliRICH;
-class TFile;
-class TNtuple;
-class TCanvas;
+class AliRICHHelix;
 
 class AliRICHRecon : public TTask 
 {
 public : 
-    AliRICHRecon(const char*, const char*);
-   ~AliRICHRecon(){EndProcessEvent();}
+    AliRICHRecon(AliRICHHelix *pHelix,TClonesArray *pClusters,Int_t iMipId);
+    virtual ~AliRICHRecon(){;}
 
-  AliRICH* Rich() {return fRich;}             //main pointer to RICH
-  void StartProcessEvent();                   //
-  void EndProcessEvent();                     //
-  //PH  void InitRecon();                           //
-  void PatRec();                              //
-  void Minimization();                        //
-  void FillHistograms();                      //
+  Double_t ThetaCerenkov();                   // it returns reconstructed Theta Cerenkov
   void FindThetaPhotonCerenkov();             //
   void FindAreaAndPortionOfRing();            //
   void FindEmissionPoint();                   //
@@ -45,12 +36,8 @@ public :
   void FindWeightThetaCerenkov();             //
   void EstimationOfTheta();                   //
   void FindIntersectionWithDetector();        //
-  void Waiting();                             //
-  Float_t FindMassOfParticle();               //
   Float_t Cerenkovangle(Float_t n, Float_t b);//
-  Float_t PhotonPositionOnCathode();          //
   Int_t   PhotonInBand();                     //
-  void DrawEvent(Int_t flag)  const;          //
   Int_t   CheckDetectorAcceptance() const;    //
   Int_t   GetFittedHoughPhotons()                   const{ return fFittedHoughPhotons;}             //
   Int_t   GetPhotonFlag()                           const{ return fPhotonFlag[fPhotonIndex];}       //
@@ -61,13 +48,10 @@ public :
   Int_t   GetTrackIndex()                           const{ return fTrackIndex;}                     //
   Int_t   GetCandidatePhotonsNumber()               const{ return fCandidatePhotonsNumber;}         //
   Int_t   GetHoughPhotons()                         const{ return fHoughPhotons;}                   //
-  Float_t GetEventVertexZ()                         const{ return fEventVertZ;}                     //
-  Float_t GetEventMultiplicity()                    const{ return fEventMultiplicity;}              //
   Float_t GetPhotonEnergy()                         const{ return fPhotonEnergy;}                   //
   Float_t GetFreonRefractiveIndex()                 const{ return fFreonRefractiveIndex;}           //
   Float_t GetQuartzRefractiveIndex()                const{ return fQuartzRefractiveIndex;}          //
   Float_t GetGasRefractiveIndex()                   const{ return fGasRefractiveIndex;}             //
-  Float_t GetFreonScaleFactor()                     const{ return fFreonScaleFactor;}               //
   Float_t GetEmissionPoint()                        const{ return fLengthEmissionPoint;}            //
   Float_t GetMassHypotesis()                        const{ return fMassHypotesis;}                  //
   Float_t GetBetaOfParticle()                       const{ return fTrackBeta;}                      //
@@ -75,13 +59,8 @@ public :
   Float_t GetEntranceY()                            const{ return fYtoentr;}                        //
   Float_t GetThetaCerenkov()                        const{ return fThetaCerenkov;}                  //
   Float_t GetThetaPhotonCerenkov()                  const{ return fThetaPhotonCerenkov;}            //
-  Float_t GetTrackMomentum()                        const{ return fTrackMomentum;}                  //
-  Float_t GetTrackEta()                             const{ return fTrackEta;}                       //
   Float_t GetTrackTheta()                           const{ return fTrackTheta;}                     //
   Float_t GetTrackPhi()                             const{ return fTrackPhi;}                       //
-  Float_t GetTrackPt()                              const{ return fTrackPt;}                        //
-  Float_t GetTrackTPCLastZ()                        const{ return fTrackTPCLastZ;}                  //
-  Float_t GetMinDist()                              const{ return fMinDist;}                        //
   Float_t GetXPointOnCathode()                      const{ return fPhotonLimitX;}                   //
   Float_t GetYPointOnCathode()                      const{ return fPhotonLimitY;}                   //
   Float_t GetThetaPhotonInDRS()                     const{ return fThetaPhotonInDRS;}               //
@@ -111,24 +90,17 @@ public :
   Float_t GetPhotonEta()                            const{ return fPhotonEta[fPhotonIndex];}        //
   Float_t GetPhotonWeight()                         const{ return fPhotonWeight[fPhotonIndex];}     //
   Float_t GetHoughRMS()                             const{ return fHoughRMS;}                       //
-  Float_t* GetCandidatePhotonX()                    const{ return fCandidatePhotonX;}               //
-  Float_t* GetCandidatePhotonY()                    const{ return fCandidatePhotonY;}               //
-  Float_t GetHoughPhotonsNorm()                     const{ return fHoughPhotonsNorm;}               //
   Float_t GetFittedTrackTheta()                     const{ return fFittedTrackTheta;}               //
   Float_t GetFittedTrackPhi()                       const{ return fFittedTrackPhi;}                 //
   Float_t GetFittedThetaCerenkov()                  const{ return fFittedThetaCerenkov;}            //
   Float_t GetEstimationOfTheta()                    const{ return fEstimationOfTheta;}              //
   Float_t GetEstimationOfThetaRMS()                 const{ return fEstimationOfThetaRMS;}           //
-  void SetEventVertexZ(Float_t EventVertZ) { fEventVertZ = EventVertZ;}                             //
-  void SetEventMultiplicity(Float_t EventMultiplicity) { fEventMultiplicity = EventMultiplicity;}   //
   void SetPhotonEnergy(Float_t PhotonEnergy) { fPhotonEnergy = PhotonEnergy;}                       //
   void SetFreonRefractiveIndex() {fFreonRefractiveIndex = fFreonScaleFactor*(1.177+0.0172*fPhotonEnergy);}//
   void SetQuartzRefractiveIndex() {fQuartzRefractiveIndex = sqrt(1+(46.411/(113.763556-TMath::Power(fPhotonEnergy,2)))+(228.71/(328.51563-TMath::Power(fPhotonEnergy,2))));}//
   void SetGasRefractiveIndex() { fGasRefractiveIndex = 1.;}                                         //
   void SetFreonScaleFactor(Float_t FreonScaleFactor) {fFreonScaleFactor = FreonScaleFactor;}        //
   void SetEmissionPoint(Float_t LengthEmissionPoint) { fLengthEmissionPoint = LengthEmissionPoint;} //
-  void SetMassHypotesis(Float_t mass) {fMassHypotesis = mass;}                                      //
-  void SetBetaOfParticle() { fTrackBeta = fTrackMomentum/sqrt(TMath::Power(fTrackMomentum,2)+TMath::Power(fMassHypotesis,2));}//
   void SetEntranceX(Float_t Xtoentr) { fXtoentr = Xtoentr;}                                         //
   void SetEntranceY(Float_t Ytoentr) { fYtoentr = Ytoentr;}                                         //
   void SetThetaPhotonInTRS(Float_t Theta) {fThetaPhotonInTRS = Theta;}                              //
@@ -149,14 +121,9 @@ public :
   void SetRadiusOuterRing(Float_t OuterRadius) {fOuterRadius = OuterRadius;}                        //
   void SetThetaCerenkov(Float_t ThetaCer) {fThetaCerenkov = ThetaCer;}                              //
   void SetThetaPhotonCerenkov(Float_t ThetaPhotCer) {fThetaPhotonCerenkov = ThetaPhotCer;}          //
-  void SetTrackMomentum(Float_t TrackMomentum) {fTrackMomentum = TrackMomentum;}                    //
-  void SetTrackEta(Float_t TrackEta) {fTrackEta = TrackEta;}                                        //
   void SetTrackTheta(Float_t TrackTheta) { fTrackTheta = TrackTheta;}                               //
   void SetTrackPhi(Float_t TrackPhi) { fTrackPhi = TrackPhi;}                                       //
-  void SetTrackPt(Float_t TrackPt) { fTrackPt = TrackPt;}                                           //
   void SetTrackCharge(Int_t TrackCharge) { fTrackCharge = TrackCharge;}                             //
-  void SetTrackTPCLastZ(Float_t TrackTPCLastZ) { fTrackTPCLastZ = TrackTPCLastZ;}                   //
-  void SetMinDist(Float_t MinDist) { fMinDist = MinDist;}                                           //
   void SetShiftX(Float_t ShiftX) { fShiftX = ShiftX;}                                               //
   void SetShiftY(Float_t ShiftY) { fShiftY = ShiftY;}                                               //
   void SetDetectorWhereX(Float_t Xcoord) { fXcoord = Xcoord;}                                       //
@@ -175,9 +142,6 @@ public :
   void SetHoughRMS(Float_t HoughRMS) { fHoughRMS = HoughRMS;}                                       //
   void SetMipIndex(Int_t MipIndex) { fMipIndex = MipIndex;}                                         //
   void SetTrackIndex(Int_t TrackIndex) { fTrackIndex = TrackIndex;}                                 //
-  void SetCandidatePhotonX(Float_t *CandidatePhotonX) { fCandidatePhotonX = CandidatePhotonX;}      //
-  void SetCandidatePhotonY(Float_t *CandidatePhotonY) { fCandidatePhotonY = CandidatePhotonY;}      //
-  void SetCandidatePhotonsNumber(Int_t CandidatePhotonsNumber) { fCandidatePhotonsNumber = CandidatePhotonsNumber;}//
   void SetHoughPhotons(Int_t HoughPhotons) { fHoughPhotons = HoughPhotons;}                         //
   void SetHoughPhotonsNorm(Float_t HoughPhotonsNorm) { fHoughPhotonsNorm = HoughPhotonsNorm;}       //
   void SetFittedTrackTheta(Float_t FittedTrackTheta)    { fFittedTrackTheta = FittedTrackTheta;}    //
@@ -191,7 +155,7 @@ public :
   Float_t FromEmissionToCathode();                                                                  //
 
 protected:
-  AliRICH* fRich;                             // main poiter to RICH
+  TClonesArray *fpClusters;                   // poiter to clusters
   Int_t   fTrackCharge;                       // charge track
   Int_t fMipIndex;                            // index for Mip
   Int_t fTrackIndex;                          // index for track
@@ -201,14 +165,9 @@ protected:
   Int_t fCandidatePhotonsNumber;              // number of candidate photons
   Int_t fHoughPhotons;                        // n. photons after Hough
   Int_t   fFittedHoughPhotons;                // n. photons after Hough and after minimization
-  Float_t fEventVertZ;                        // z coord. of the primary vertex
-  Float_t fEventMultiplicity;                 // event primary multiplicity
+
   Float_t fTrackTheta;                        // Theta of track at RICH
   Float_t fTrackPhi;                          // Phi of track at RICH
-  Float_t fTrackMomentum;                     // track momentum
-  Float_t fTrackEta;                          // track pseudorapidity
-  Float_t fTrackPt;                           // pt of track
-  Float_t fTrackTPCLastZ;                     // z last point of the TPC
   Float_t fMinDist;                           // min distance between extrapolated track and MIP
   Float_t fTrackBeta;                         // beta of the track
   Float_t fXtoentr;                           // X entrance to RICH
@@ -236,7 +195,6 @@ protected:
   Float_t fPhotonLimitX;                      // X phys limit for photon
   Float_t fPhotonLimitY;                      // Y phys limit for photon 
   Float_t fDistanceFromCluster;               // distance from cluster
-  Float_t fMassHypotesis;                     // reconstructed mass
   Float_t fCerenkovAnglePad;                  // cherenkov angle of pad
   Float_t fThetaPhotonCerenkov;               // theta cerenkov for photon
   Float_t fShiftX;                            // x shift to entrance in radiator
@@ -245,6 +203,7 @@ protected:
   Float_t fYcoord;                            // ..
   Float_t fIntersectionX;                     // ..
   Float_t fIntersectionY;                     // ..
+  Float_t fMassHypotesis;                     //
   Float_t fThetaOfRing;                       // theta of ring
   Float_t fAreaOfRing;                        // area of the ring
   Float_t fPortionOfRing;                     // fraction of the accepted ring
@@ -263,19 +222,9 @@ protected:
   Int_t   fThetaBin;                          // bin in theta
   Float_t fThetaMin,fThetaMax;                // min max
   Float_t fXmin,fXmax,fYmin,fYmax;            // xy min max
-  TFile  *fOutFile;                           // histogram output
-  TNtuple *fNtuple;                           // output ntuple
-  TCanvas *fDisplay;                          // for display
   Int_t   fNrings;                            //current number of reconstructed rings
-  Bool_t  fDebug;                             // flag for debug
-  Bool_t  fIsDISPLAY;                         // flag for display
   Bool_t  fIsWEIGHT;                          // flag to consider weight procedure
   Bool_t  fIsBACKGROUND;                      // flag to simulate bkg
-  Bool_t  fIsMINIMIZER;                       // flag to start with (theta,phi) minimization
-  Int_t   fNpadX;                             // npadx
-  Int_t   fNpadY;                             // npady
-  Float_t fPadSizeX;                          // sizepad x
-  Float_t fPadSizeY;                          // sizepad y
   Float_t fRadiatorWidth;                     // radiator width
   Float_t fQuartzWidth;                       // quartz width
   Float_t fGapWidth;                          // gap width
index 581478b..47424a9 100644 (file)
@@ -13,7 +13,6 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/* $Id$ */
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
 
 
 #include "AliRICHReconstructor.h"
-#include "AliRunLoader.h"
-#include "AliRun.h"
 #include "AliRICHClusterFinder.h"
-
+#include "AliRICHHelix.h"
+#include <AliRunLoader.h>
+#include <AliRun.h>
+#include <AliESD.h>
 
 ClassImp(AliRICHReconstructor)
 
-
-//_____________________________________________________________________________
-void AliRICHReconstructor::Reconstruct(AliRunLoader* runLoader) const
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::Reconstruct(AliRunLoader* pAL) const
 {
-// reconstruct clusters
-
-  AliRICH* rich = GetRICH(runLoader);
-  if (!rich) return;
-  AliRICHClusterFinder clusterer(rich);
-  clusterer.Exec();
+//Finds clusters out of digits
+  AliDebug(1,"Start cluster finder.");AliRICHClusterFinder clus(GetRICH(pAL));  clus.Exec();
 }
-
-//_____________________________________________________________________________
-void AliRICHReconstructor::FillESD(AliRunLoader* /*runLoader*/, 
-                                  AliESD* /*esd*/) const
+//__________________________________________________________________________________________________
+void AliRICHReconstructor::FillESD(AliRunLoader* /*pAL*/, AliESD* pESD) const
 {
-// nothing to be done
-
-}
-
+//This methode fills AliESDtrack with information from RICH  
+  AliDebug(1,Form("Start with %i tracks",pESD->GetNumberOfTracks()));
+/*  const Double_t masses[5]={0.000511,0.105658,0.139567,0.493677,0.93828};//electron,muon,pion,kaon,proton
+  const Double_t refIndex = 1.29052;
 
-//_____________________________________________________________________________
-AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* runLoader) const
+  Double_t thetaTh[5];
+  Double_t sinThetaThNorm;
+  Double_t sigmaThetaTh[5];
+  Double_t height[5];
+  Double_t totalHeight=0;
+  Double_t x[3],p[3]; //tmp storage for track parameters
+  for(Int_t iTrackN=0;iTrackN<pESD->GetNumberOfTracks();iTrackN++){//ESD tracks loop
+    AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
+//    if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
+    pTrack->GetXYZ(x); pTrack->GetPxPyPz(p);          Double_t pmod=pTrack->GetP();//get running track parameters
+    continue;  
+//  AliRICHHelix helix(x[0],x[1],x[2],p[0],p[1],p[2]);      //construct helix from running track parameters
+    
+//    TVector rad(1,5); TVector pc(1,5);
+//    helix.RichIntersect(GetRICH(pAL)->P(),rad,pc); //returns cross point of track with RICH PC in LRS
+    
+    for(Int_t iPart=4;iPart>=0;iPart--){
+      Double_t cosThetaTh = TMath::Sqrt(masses[iPart]*masses[iPart]+pmod*pmod)/(refIndex*pmod);
+      if(cosThetaTh>=1) {break;}
+      thetaTh[iPart] = TMath::ACos(cosThetaTh);
+      sinThetaThNorm = TMath::Sin(thetaTh[iPart])/TMath::Sqrt(1-1/(refIndex*refIndex));
+      sigmaThetaTh[iPart] = (0.014*(1/sinThetaThNorm-1) + 0.0043)*1.25;
+      height[iPart] = TMath::Gaus(thetaExp,thetaTh[iPart],sigmaThetaTh[iPart]);
+      totalHeight +=height[iPart];
+    }
+    
+    Double_t richPID[5];
+    for(Int_t iPart=0;iPart<5;iPart++) richPID[iPart] = height[iPart]/totalHeight;    
+    pTrack->SetRICHpid(richPID); 
+  }//ESD tracks loop
+  */
+}//FillESD
+//__________________________________________________________________________________________________
+AliRICH* AliRICHReconstructor::GetRICH(AliRunLoader* pAL) const
 {
 // get the RICH detector
 
-  if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
-  if (!runLoader->GetAliRun()) {
-    Error("GetRICH", "couldn't get AliRun object");
-    return NULL;
-  }
-  AliRICH* rich = (AliRICH*) runLoader->GetAliRun()->GetDetector("RICH");
-  if (!rich) {
-    Error("GetRICH", "couldn't get RICH detector");
-    return NULL;
-  }
-  return rich;
+  if (!pAL->GetAliRun()) pAL->LoadgAlice();
+  if (!pAL->GetAliRun()) {AliError("couldn't get AliRun object"); return NULL;  }
+  AliRICH* pRich = (AliRICH*) pAL->GetAliRun()->GetDetector("RICH");
+  if (!pRich) {AliError("couldn't get RICH detector");    return NULL;  }
+  return pRich;
 }
+//__________________________________________________________________________________________________
index a17c42d..6bd045d 100644 (file)
@@ -1,27 +1,26 @@
-#ifndef ALIRICHRECONSTRUCTOR_H
-#define ALIRICHRECONSTRUCTOR_H
+#ifndef AliRICHReconstructor_h
+#define AliRICHReconstructor_h
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id$ */
-
-#include "AliReconstructor.h"
+#include <AliReconstructor.h>
+#include <AliLog.h>
+#include "AliRICHTracker.h"
 
 class AliRICH;
 
-
-class AliRICHReconstructor: public AliReconstructor {
+class AliRICHReconstructor: public AliReconstructor 
+{
 public:
-  AliRICHReconstructor(): AliReconstructor() {};
-  virtual ~AliRICHReconstructor() {};
-
-  virtual void         Reconstruct(AliRunLoader* runLoader) const;
-  virtual void         FillESD(AliRunLoader* runLoader, AliESD* esd) const;
-
+           AliRICHReconstructor(): AliReconstructor()          {AliDebug(1,"Start.");}
+  virtual ~AliRICHReconstructor()                              {AliDebug(1,"Start.");}  
+  virtual AliTracker*  CreateTracker(AliRunLoader*)const       {AliDebug(1,"Start.");return new AliRICHTracker;}    //virtual from AliReconstructor
+  virtual void         Reconstruct(AliRunLoader* pAL) const;              //virtual from AliReconstructor
+  virtual void         FillESD(AliRunLoader* pAL, AliESD* pESD) const;    //virtual from AliReconstructor
 private:
-  AliRICH*             GetRICH(AliRunLoader* runLoader) const;
+  AliRICH*             GetRICH(AliRunLoader* pAL) const;
 
-  ClassDef(AliRICHReconstructor, 0)   // class for the RICH reconstruction
+  ClassDef(AliRICHReconstructor, 0)   //class for the RICH reconstruction
 };
 
 #endif
diff --git a/RICH/AliRICHTracker.cxx b/RICH/AliRICHTracker.cxx
new file mode 100644 (file)
index 0000000..0208b58
--- /dev/null
@@ -0,0 +1,72 @@
+#include "AliRICHTracker.h"
+#include <AliESD.h>
+#include <TVector3.h>
+#include <TTree.h>
+#include "AliRICH.h"
+#include "AliRICHHelix.h"
+#include <AliMagF.h>
+#include "AliRICHRecon.h"
+#include <AliStack.h>
+#include <TParticle.h>
+
+ClassImp(AliRICHTracker)
+
+
+Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
+{
+  Int_t iNtracks=pESD->GetNumberOfTracks();
+  Double_t b=GetFieldMap()->SolenoidField()/10;// magnetic field in Tesla
+  AliDebug(1,Form("Start with %i tracks in %f Tesla field",iNtracks,b));
+  Double_t xb[3],pb[3];//tmp storage for track parameters
+  TVector3 x0(0,0,0); TVector3 p0(0,0,0);//tmp storage for AliRICHHelix
+  
+  AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));
+  
+//  pRich->GetLoader()->GetRunLoader()->LoadHeader();pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+//  AliStack *pStack =   pRich->GetLoader()->GetRunLoader()->Stack();
+//  TParticle *pParticle = pStack->Particle(0);
+//  p0.SetMagThetaPhi(pParticle->P(),pParticle->Theta(),pParticle->Phi());
+
+  for(Int_t iTrackN=0;iTrackN<iNtracks;iTrackN++){//ESD tracks loop
+    AliESDtrack *pTrack = pESD->GetTrack(iTrackN);// get next reconstructed track
+//  if((pTrack->GetStatus()&AliESDtrack::kTOFout)==0) continue; //ignore tracks not recontructed by TOF
+    pTrack->GetXYZ(xb); 
+      pTrack->GetPxPyPz(pb); 
+          Int_t status=pTrack->GetStatus()&AliESDtrack::kTOFout;//get running track parameters
+    AliDebug(1,Form("Track %i pmod=%f mass=%f stat=%i",iTrackN,pTrack->GetP(),pTrack->GetMass(),status));
+//      x.Print();p.Print();
+    x0.SetXYZ(xb[0],xb[1],xb[2]); p0.SetXYZ(xb[0],xb[1],xb[2]);
+    AliRICHHelix helix(x0,p0,pTrack->GetSign(),b);   
+    Int_t iChamber=helix.RichIntersect(pRich->P());        
+    AliDebug(1,Form("intersection with %i chamber found",iChamber));
+    if(!iChamber) continue;//intersection with no chamber found
+    
+    Double_t distMip=9999; //min distance between clusters and track position on PC 
+    Int_t iMipId=0; //index of that min distance cluster 
+    for(Int_t iClusN=0;iClusN<pRich->Clusters(iChamber)->GetEntries();iClusN++){//clusters loop for intersected chamber
+      AliRICHcluster *pClus=(AliRICHcluster*)pRich->Clusters(iChamber)->UncheckedAt(iClusN);//get pointer to current cluster
+      Double_t distCurrent=pClus->DistTo(helix.PosPc());//ditance between current cluster and helix intersection with PC
+      if(distCurrent<distMip){distMip=distCurrent;iMipId=iClusN;}//find cluster nearest to the track 
+      
+      AliDebug(1,Form("Ploc (%f,%f,%f) dist= %f",helix.Ploc().Mag(),helix.Ploc().Theta()*TMath::RadToDeg(),
+                                                                    helix.Ploc().Phi()*TMath::RadToDeg(),pClus->DistTo(helix.PosPc())));
+    }////clusters loop for intersected chamber
+    
+    AliDebug(1,Form("Min distance cluster: %i dist is %f",iMipId,distMip));
+    
+    AliRICHRecon recon(&helix,pRich->Clusters(iChamber),iMipId);
+    Double_t thetaCerenkov=recon.ThetaCerenkov(); //search for mean Cerenkov angle for this track
+    AliDebug(1,Form("FINAL Theta Cerenkov=%f",thetaCerenkov));
+    pTrack->SetRICHsignal(thetaCerenkov);
+        
+    
+  }//ESD tracks loop
+  AliDebug(1,"Stop.");  
+  return 0;
+} //pure virtual from AliTracker
+
+Int_t AliRICHTracker::LoadClusters(TTree *pTree)
+{
+// Load clusters for RICH
+  AliDebug(1,"Start.");  pTree->GetEntry(0);  AliDebug(1,"Stop."); return 0;
+}
diff --git a/RICH/AliRICHTracker.h b/RICH/AliRICHTracker.h
new file mode 100644 (file)
index 0000000..7fd064d
--- /dev/null
@@ -0,0 +1,28 @@
+#ifndef AliRICHTracker_h
+#define AliRICHTracker_h
+
+#include <AliTracker.h>
+#include <AliLog.h>
+
+class AliCluster;
+class AliESD;
+class TTree;
+
+class AliRICHTracker : public AliTracker
+{
+public:
+           AliRICHTracker() :AliTracker()      {AliDebug(1,"Start.");}
+  virtual ~AliRICHTracker()                    {AliDebug(1,"Start.");}
+  Int_t Clusters2Tracks(AliESD *)              {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker 
+  Int_t RefitInward(AliESD *)                  {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker 
+  void UnloadClusters()                        {AliDebug(1,"Start.");}          //pure virtual from AliTracker 
+  AliCluster *GetCluster(Int_t )const          {AliDebug(1,"Start.");return 0;} //pure virtual from AliTracker 
+  Int_t PropagateBack(AliESD *);                                                //pure virtual from AliTracker 
+  Int_t LoadClusters(TTree *);                                                  //pure virtual from AliTracker 
+
+protected:
+
+  ClassDef(AliRICHTracker,0)
+};//class AliRICHTracker
+
+#endif//AliRICHTracker_h
index 9bfebf4..06adbc8 100644 (file)
@@ -29,25 +29,21 @@ ClassImp(AliRICHv0)
 void AliRICHv0::StepManager()
 {
 //This StepManager is a provision for different test-learn activities on the current MC layer  
+  static Int_t iStepN;
   const char *sParticle;
   switch(gMC->TrackPid()){
-    case kProton:
-      sParticle="proton";break;
-    case kNeutron:
-      sParticle="neutron";break;
-    case kGamma:
-      sParticle="gamma";break;
-    case kCerenkov:
-      sParticle="photon";break;
-    case kPi0:
-      sParticle="Pi0";break;  
-    case kElectron:
-      sParticle="electron";break;  
-    default:
-      sParticle="not known";break;
+    case kProton:      sParticle="proton"    ;break;
+    case kNeutron:     sParticle="neutron"   ;break;
+    case kGamma:       sParticle="gamma"     ;break;
+    case kCerenkov:    sParticle="photon"    ;break;
+    case kPi0:         sParticle="Pi0"       ;break;  
+    case kPiPlus:      sParticle="Pi+"       ;break;  
+    case kPiMinus:     sParticle="Pi-"       ;break;  
+    case kElectron:    sParticle="electron"  ;break;  
+    default:           sParticle="not known" ;break;
   }
 
-  Info("","event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f",
+  Info(Form("Step %i",iStepN),"event=%i hunt=%i tid=%i pid=%i(%s) m=%f q=%3.1f dEdX=%9.3f",
                             gMC->CurrentEvent(),
                             fIshunt,
                             gAlice->GetMCApp()->GetCurrentTrackNumber(),
@@ -55,8 +51,8 @@ void AliRICHv0::StepManager()
                             sParticle,
                             gMC->TrackMass(),
                             gMC->TrackCharge(),
-                            gMC->Edep());
-  Info("","Flags:alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)",
+                            gMC->Edep()*1e9);
+  Info("Flags","alive(%i) disap(%i) enter(%i) exit(%i) inside(%i) out(%i) stop(%i) new(%i)",
                             gMC->IsTrackAlive(),
                             gMC->IsTrackDisappeared(),
                             gMC->IsTrackEntering(),
@@ -65,12 +61,12 @@ void AliRICHv0::StepManager()
                             gMC->IsTrackOut(),
                             gMC->IsTrackStop(),
                             gMC->IsNewTrack());
-  Int_t copy0,copy1,copy2,copy3;
+  Int_t copy0=0,copy1=0,copy2=0,copy3=0;
   Int_t vid0=gMC->CurrentVolID(copy0);
   Int_t vid1=gMC->CurrentVolOffID(1,copy1);
   Int_t vid2=gMC->CurrentVolOffID(2,copy2);
   Int_t vid3=gMC->CurrentVolOffID(3,copy3);
-  Info("","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i   %s-%s-%s-%s",
+  Info("Volumes","vid0=%i(%s)c%i vid1=%i(%s)c%i vid2=%i(%s)c%i vid3=%i(%s)c%i   %s-%s-%s-%s",
                       vid0,gMC->VolName(vid0),copy0, 
                       vid1,gMC->VolName(vid1),copy1, 
                       vid2,gMC->VolName(vid2),copy2, 
@@ -82,72 +78,44 @@ void AliRICHv0::StepManager()
   
   Float_t a,z,den,rad,abs; a=z=den=rad=abs=kBad;
   Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
-  Info("","mid=%i a=%7.2f z=%7.2f den=%7.2f rad=%7.2f abs=%7.2f",mid,a,z,den,rad,abs);
+  Info("Material","id=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f",mid,a,z,den,rad,abs);
+  
+  Int_t iTmedId=gMC->GetMedium();
+  const char *sTmed;
+  switch(iTmedId){
+    case kAir:        sTmed="Air"      ;break;
+    case kRoha:       sTmed="Rohacell" ;break;
+    case kSiO2:       sTmed="SiO2"     ;break;
+    case kC6F14:      sTmed="C6F14"    ;break;
+    case kGridCu:     sTmed="GridCu"   ;break;
+    case kOpSiO2:     sTmed="OpSiO2"   ;break;
+    case kGap:        sTmed="Gap"      ;break;
+    case kGlass:      sTmed="Glass"    ;break;
+    case kCu:         sTmed="Cu"       ;break;
+    case kSteel:      sTmed="Steel"    ;break;
+    case kPerpex:     sTmed="Perpex"   ;break;
+    case kCH4:        sTmed="CH4"      ;break;
+    case kCsI:        sTmed="CsI"      ;break;
+    case kAl:         sTmed="Al"       ;break;
+    case kW:          sTmed="W"        ;break;
+    default:          sTmed="not known";break;
+  }
+  Info("Medium","id=%i (%s)",iTmedId,sTmed);
   
   TLorentzVector x4;
   gMC->TrackPosition(x4);
   Float_t glo[3],loc[3];
   glo[0]=x4.X();glo[1]=x4.Y();glo[2]=x4.Z();  
   gMC->Gmtod(glo,loc,1);
-  Info("","glo(%+8.3f,%+8.3f,%+8.3f) r=%8.3f theta=%8.3f phi=%8.3f",
+  Info("X4 glo","(x,y,z)=(%+8.3f,%+8.3f,%+8.3f) (r,theta,phi)=(%8.3f,%8.3f,%8.3f)",
                       glo[0],glo[1],glo[2],x4.Rho(),x4.Theta()*TMath::RadToDeg(),x4.Phi()*TMath::RadToDeg());  
-  Info("","loc(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()",         loc[0],loc[1],loc[2]);  
-  if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy0)){
+  Info("X4 loc","(x,y,z)=(%+8.3f,%+8.3f,%8.3f) by gMC->Gmtod()",loc[0],loc[1],loc[2]);  
+  if(gMC->VolId("GAP ")==gMC->CurrentVolID(copy0)){
     Int_t iChamber;
     gMC->CurrentVolOffID(2,iChamber);
-    TVector2 x2=C(iChamber)->Glob2Loc(x4);
-    Info("","loc(%+8.3f,%+8.3f) by Glob2Loc",      x2.X(),x2.Y());  
+    TVector2 x2=C(iChamber)->Mrs2Pc(x4);
+    Info("X4 loc","(%+8.3f,%+8.3f) by Mrs2Pc",      x2.X(),x2.Y());  
   }
-  Info("","end of current step\n");
+  Info(Form("Step %i",iStepN++),"end of current step\n");
 }//StepManager()
 
-void AliRICHv0::CreateGeometry()
-{
-  if(GetDebug())Info("CreateGeometry","Start v0.");  
-  
-  Double_t cm=1,mm=0.1;
-  //place radioactive source compartment if needed
-  Double_t containerLen=21.8*mm/2                        ,containerR=38*mm/2;  
-  Double_t screwLen=15*mm/2                              ,screwR=5*mm/2;  
-  Double_t srcLen=10*mm/2                                ,srcR=2*mm/2;  
-  Double_t perpexLen=20*mm/2                             ,perpexR=34*mm/2;  
-  Double_t perpexWholeLen=10*mm/2                        ,perpexWholeR=4*mm/2;  
-  Double_t alBottomLen=containerLen-perpexLen            ,alWholeR=5*mm/2;      
-//volumes    
-  Double_t par[3];
-
-//  Double_t anodWireD=20*mkm, cathWireD=100*mkm, collWireD=100*mkm;
-        
-  Double_t pcZ2=0.5*mm;
-  par[0]=68.8*cm;            par[1]=70.86*cm;                par[2]=13*cm;           gMC->Gsvolu("RICH","BOX ",(*fIdtmed)[kCH4],par,3); //RICH    
-  par[0]=P()->PcSizeX()/2;   par[1]=P()->PcSizeY()/2;        par[2]=pcZ2;            gMC->Gsvolu("CSI ","BOX ",(*fIdtmed)[kCSI],par,3); //CSI
-  par[0]=P()->PcSizeX()/2;   par[1]=P()->PcSizeY()/2;        par[2]=P()->GapAmp();   gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kCH4],par,3); //GAP
-//  par[0]=0;                  par[1]=;                        par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WIAN","TUBE",(*fIdtmed)[kW],par,3);//Anod wire
-//  par[0]=0;                  par[1]=srcR;                    par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICA","TUBE",(*fIdtmed)[kCu],par,3);//Cathod wire
-//  par[0]=0;                  par[1]=srcR;                    par[2]=P()->PcSizeX()/2;gMC->Gsvolu("WICO","TUBE",(*fIdtmed)[kCu],par,3);//Collect wire
-    
-  par[0]=0      ;par[1]=containerR  ;par[2]=containerLen;  gMC->Gsvolu("CONT","TUBE",(*fIdtmed)[kCH4]   ,par,3); //container 
-  par[0]=perpexR;par[1]=containerR  ;par[2]=containerLen;  gMC->Gsvolu("ALWA","TUBE",(*fIdtmed)[kAl]    ,par,3); //Al cylindric wall
-  par[0]=0      ;par[1]=perpexR     ;par[2]=alBottomLen;   gMC->Gsvolu("ALBO","TUBE",(*fIdtmed)[kAl]    ,par,3); //Al bottom 
-  par[0]=0      ;par[1]=alWholeR    ;par[2]=alBottomLen;   gMC->Gsvolu("ALWH","TUBE",(*fIdtmed)[kCH4]   ,par,3); //Whole in Al bottom
-  par[0]=0      ;par[1]=perpexR     ;par[2]=perpexLen;     gMC->Gsvolu("PEPL","TUBE",(*fIdtmed)[kPerpex],par,3); //Perpex plug
-  par[0]=0      ;par[1]=perpexWholeR;par[2]=perpexWholeLen;gMC->Gsvolu("PEWH","TUBE",(*fIdtmed)[kCH4]   ,par,3); //Whole in Perpex
-  par[0]=0      ;par[1]=screwR      ;par[2]=screwLen;      gMC->Gsvolu("SCRE","TUBE",(*fIdtmed)[kSteel] ,par,3); //Screw
-  par[0]=0      ;par[1]=srcR        ;par[2]=srcLen;        gMC->Gsvolu("SOUR","TUBE",(*fIdtmed)[kSteel] ,par,3); //Source itself
-//nodes    
-  gMC->Gspos("RICH",1,"ALIC",0,0,-9                                ,0,"ONLY");     //RICH in ALIC
-  gMC->Gspos("CSI ",1,"RICH",0,0,-(P()->GapProx()+pcZ2)           ,0,"ONLY");     //CsI  in ALIC
-  gMC->Gspos("GAP ",1,"RICH",0,0,-(P()->GapProx()-P()->GapAmp()/2),0,"ONLY");     //GAP in ALIC
-    
-  gMC->Gspos("CONT",1,"RICH",0,0,1*cm                             ,0,"ONLY");     //Sr90 container in RICH
-  gMC->Gspos("ALWA",1,"CONT",0,0,0                                ,0,"ONLY");     //Al wall
-  gMC->Gspos("ALBO",1,"CONT",0,0,-containerLen+alBottomLen        ,0,"ONLY");     //Al bottom
-  gMC->Gspos("PEPL",1,"CONT",0,0,containerLen-perpexLen           ,0,"ONLY");     //Perpex plug
-   
-  gMC->Gspos("ALWH",1,"ALBO",6*mm,0,0                             ,0,"ONLY");     //Whole in Al bottom
-    
-  gMC->Gspos("PEWH",1,"PEPL",6*mm,0,-perpexLen+perpexWholeLen     ,0,"ONLY");     //Whole in Perpex plug
-  gMC->Gspos("SOUR",1,"PEPL",6*mm,0, perpexLen-srcLen             ,0,"ONLY");     //Source in Perpex plug
-  gMC->Gspos("SCRE",1,"PEPL",0,0,perpexLen-screwLen               ,0,"ONLY");     //Screw in Perpex plug
-  if(GetDebug())Info("CreateGeometry","Stop v0.");  
-}//CreateGeometry()
index 3594ff8..bf5ec4b 100644 (file)
@@ -15,7 +15,6 @@ public:
   virtual void   Init()                                              {;}                      //interface from AliRICH
   virtual Int_t  IsVersion()                                    const{return 0;}              //interface from AliRICH
   virtual void   StepManager();                                                               //interface from AliRICH
-  virtual void   CreateGeometry();                                                            //interface from AliRICH
 protected:
   ClassDef(AliRICHv0,1)  //RICH coarse version for material budget study and debuging
 };
index 8d9efda..cca5906 100644 (file)
@@ -32,42 +32,42 @@ ClassImp(AliRICHv1)
 void AliRICHv1::StepManager()
 {
 // Full Step Manager.
-// 3- Ckovs absorbed in Collection electrods
-// 5- Ckovs absorbed in Cathode wires
-// 6- Ckovs absorbed in Anod wires
+// 3- Ckovs absorbed on Collection electrods
+// 5- Ckovs absorbed on Cathode wires
+// 6- Ckovs absorbed on Anod wires
          
   Int_t          copy;
   static Int_t   iCurrentChamber;
 //history of Cerenkovs
   if(gMC->TrackPid()==kCerenkov){
-    if(gMC->IsNewTrack() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(0)++;// 0- Ckovs produced in Freon
-    if(!gMC->IsTrackAlive() && (gMC->CurrentVolID(copy)==gMC->VolId("FRE1")||gMC->CurrentVolID(copy)==gMC->VolId("FRE2"))) fCounters(1)++;// 1- Ckovs absorbed in Freon
-    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("QUAR")) fCounters(2)++;//2- Ckovs absorbed in Quartz
-    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("META")) fCounters(4)++;//4- Ckovs absorbed in CH4
+    if( gMC->IsNewTrack()   && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(0)++;// 0- Ckovs produced in radiator
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRAD")) fCounters(1)++;// 1- Ckovs absorbed in radiator
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RRWI")) fCounters(2)++;// 2- Ckovs absorbed in radiator window
+    if(!gMC->IsTrackAlive() && gMC->CurrentVolID(copy)==gMC->VolId("RICH")) fCounters(4)++;// 4- Ckovs absorbed in CH4
   }
           
 //Treat photons    
   static TLorentzVector cerX4;
-  if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("CSI ")){//photon in CSI
-    if(gMC->Edep()>0){//CF+CSI+DE
+  if((gMC->TrackPid()==kCerenkov||gMC->TrackPid()==kFeedback)&&gMC->CurrentVolID(copy)==gMC->VolId("RPC ")){//photon in PC
+    if(gMC->Edep()>0){//photon in PC +DE
       if(IsLostByFresnel()){ 
-        if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI
+        if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected from CsI
         gMC->StopTrack();
         return;
       }        
-      gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);
+      gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);//RICH-RPPF-RPC
        
       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
-      fCounters(8)++;//4- Ckovs converted in CsI
+      fCounters(8)++;//4- Ckovs converted to electron on CsI
       GenerateFeedbacks(iCurrentChamber);
-    }//CF+CSI+DE
-  }//CF in CSI
+    }//photon in PC and DE >0 
+  }//photon in PC
   
 //Treat charged particles  
   static Float_t eloss;
   static TLorentzVector mipInX4,mipOutX4;
-  if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("GAP ")){//MIP in GAP
-    gMC->CurrentVolOffID(3,iCurrentChamber);
+  if(gMC->TrackCharge() && gMC->CurrentVolID(copy)==gMC->VolId("RGAP")){//MIP in GAP
+    gMC->CurrentVolOffID(1,iCurrentChamber);//RICH-RGAP
     if(gMC->IsTrackEntering()||gMC->IsNewTrack()) {//MIP in GAP entering or newly created
       eloss=0;                                                           
       gMC->TrackPosition(mipInX4);
@@ -93,7 +93,7 @@ Bool_t AliRICHv1::IsLostByFresnel()
   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc),localMom[1]);
   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
   if(gMC->GetRandom()->Rndm() < Fresnel(p4.E()*1e9,cotheta,1)){
-    if(GetDebug()) Info("IsLostByFresnel","Photon lost");
+    AliDebug(1,"Photon lost");
     return kTRUE;
   }else
     return kFALSE;
@@ -106,11 +106,12 @@ void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
   
   TLorentzVector x4;
   gMC->TrackPosition(x4);  
-  TVector2 x2=C(iChamber)->Glob2Loc(x4);
-  Int_t sector=P()->Loc2Sec(x2);  if(sector==kBad) return; //hit in dead zone, nothing to produce
+  TVector2 x2=C(iChamber)->Mrs2Pc(x4);//hit position on photocathode plane
+  TVector2 xspe=x2;
+  Int_t sector=P()->Loc2Sec(xspe);  if(sector==kBad) return; //hit in dead zone, nothing to produce
   Int_t iTotQdc=P()->TotQdc(x2,eloss);
-  Int_t iNphotons=gMC->GetRandom()->Poisson(P()->AlphaFeedback(sector)*iTotQdc);    
-  if(GetDebug())Info("GenerateFeedbacks","N photons=%i",iNphotons);
+  Int_t iNphotons=gMC->GetRandom()->Poisson(P()->C(iChamber)->AlphaFeedback(sector)*iTotQdc);    
+  AliDebug(1,Form("N photons=%i",iNphotons));
   Int_t j;
   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
 //Generate photons
@@ -172,5 +173,5 @@ void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
                      outputNtracksStored,
                      1.0);    
   }//feedbacks loop
-  if(GetDebug()) Info("GenerateFeedbacks","Stop.");
+  AliDebug(1,"Stop.");
 }//GenerateFeedbacks()
index 1634afc..a1b5e53 100644 (file)
@@ -6,7 +6,7 @@ public:
                    kPMD=8192,kDIPO=16384,kEMCAL=32768,kVZERO=65536,kMUON=131072,kZDC=262144,kSHILD=524288};
   enum EProcesses {kDCAY=1,kPAIR=2,kCOMP=4,kPHOT=8,kPFIS=16,kDRAY=32,kANNI=64,kBREM=128,kMUNU=256,kCKOV=512,kHADR=1024,kLOSS=2048,kMULS=4096,
                    kRAYL=8192};
-  enum EGenTypes  {kGun1,kGun7,kPythia7,kHijing,kHijingPara};
+  enum EGenTypes  {kGun1=100,kGun7,kPythia7,kHijing,kHijingPara};
   
           KirConfig(const char*sFileName);
          ~KirConfig()                    {Info("ctor","");fMain->Cleanup(); delete fMain;}
@@ -51,9 +51,12 @@ KirConfig::KirConfig(const char *sFileName)
 //Generator  
   pVerFrame->AddFrame(pGenGrpFrm=new TGGroupFrame(pHorFrame,"Generator"));
   pGenGrpFrm->AddFrame(fGenTypeCombo=new TGComboBox(pGenGrpFrm,100));
-  fGenTypeCombo->AddEntry("gun to single chamber",kGun1);  fGenTypeCombo->AddEntry("gun to all chambers",kGun7);
+  fGenTypeCombo->AddEntry("gun to single chamber",kGun1);  
+  fGenTypeCombo->AddEntry("gun to all chambers",kGun7);
   fGenTypeCombo->AddEntry("7 guns on top of Pythia",kPythia7);
-  fGenTypeCombo->AddEntry("HIJING",kHijing);                fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara);
+  fGenTypeCombo->AddEntry("HIJING",kHijing);
+  fGenTypeCombo->AddEntry("parametrized HIJING",kHijingPara);
+  fGenTypeCombo->AddEntry("Sr90 source",kSr90);
   fGenTypeCombo->Select(kHijingPara);
   fGenTypeCombo->Resize(160,20);
   
@@ -86,7 +89,7 @@ KirConfig::KirConfig(const char *sFileName)
 // Magnetic Field
   pVerFrame->AddFrame(pFldGrpFrm=new TGGroupFrame(pHorFrame,"Magnetic Field"));
   pFldGrpFrm->AddFrame(fMagFldChkBtn=new TGCheckButton(pFldGrpFrm,"On/Off"));
-  fMagFldChkBtn->SetState(kButtonUp);
+  fMagFldChkBtn->SetState(kButtonDown);
 //Detectors
   pHorFrame->AddFrame(fDetButGrp=new TGButtonGroup(pHorFrame,"Detectors"));
   fDetButGrp->Connect("Pressed(Int_t)","KirConfig",this,"AddDetector(Int_t)");
@@ -206,14 +209,15 @@ void KirConfig::CreateConfigFile()
   switch(fGenTypeCombo->GetSelected()){
     case kHijingPara: 
       fprintf(fp,"  AliGenHIJINGpara *pGen=new AliGenHIJINGpara(%i);\n",(int)fGenNprimEntry->GetNumber());
-      fprintf(fp,"  pGen->SetMomentumRange(0,999); pGen->SetPhiRange(0,360); pGen->SetThetaRange(%f,%f);\n",Eta2Theta(8),Eta2Theta(-8));
+      fprintf(fp,"  pGen->SetMomentumRange(0,999); pGen->SetThetaRange(%f,%f); pGen->SetPhiRange(0,360);\n",Eta2Theta(8),Eta2Theta(-8));
       fprintf(fp,"  pGen->SetOrigin(0,0,0);  pGen->SetSigma(0,0,0);\n");
       fprintf(fp,"  pGen->Init();\n");
     break;
     case kGun1:     
       fprintf(fp,"  AliGenFixed *pGen=new AliGenFixed(1);\n");  
-      fprintf(fp,"  pGen->SetPart(%i); pGen->SetMomentum(%3.1f); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected(),float(fGenMomCombo->GetSelected())/10);  
-      fprintf(fp,"  pGen->SetPhiRange(pRICH->C(%i)->PhiD()); pGen->SetThetaRange(pRICH->C(%i)->ThetaD()-2);\n",fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());            
+      fprintf(fp,"  pGen->SetPart(%i); pGen->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());  
+      fprintf(fp,"  pGen->SetMomentum(%3.1f); pGen->SetTheta(pRICH->C(%i)->ThetaD()-2); pGen->SetPhi(pRICH->C(%i)->PhiD()); \n",
+                    float(fGenMomCombo->GetSelected())/10,fGenChamberCombo->GetSelected(),fGenChamberCombo->GetSelected());            
       fprintf(fp,"  pGen->Init();\n");
     break;    
     case kGun7:   
@@ -221,7 +225,7 @@ void KirConfig::CreateConfigFile()
       fprintf(fp,"  for(int i=1;i<=7;i++){\n");
       fprintf(fp,"    AliGenFixed *pFixed=new AliGenFixed(1);\n");
       fprintf(fp,"    pFixed->SetPart(%i); pFixed->SetMomentum(2.5+i*0.4); pFixed->SetOrigin(0,0,0);\n",fGenPartIdCombo->GetSelected());
-      fprintf(fp,"    pFixed->SetPhiRange(pRICH->C(i)->PhiD()); pFixed->SetThetaRange(pRICH->C(i)->ThetaD()-2);\n");                             
+      fprintf(fp,"    pFixed->SetPhi(pRICH->C(i)->PhiD()); pFixed->SetTheta(pRICH->C(i)->ThetaD()-2);\n");                             
       fprintf(fp,"    pCocktail->AddGenerator(pFixed,Form(\"Fixed %i\",i),1);\n  }\n");  
       fprintf(fp,"  pCocktail->Init();\n");
     break;
@@ -238,6 +242,18 @@ void KirConfig::CreateConfigFile()
       fprintf(fp,"  pPythia->SetProcess(kPyMb);  pPythia->SetEnergyCMS(14000);\n");      
       fprintf(fp,"  pCocktail->AddGenerator(pPythia,\"Pythia\",1);\n");  
       fprintf(fp,"  pCocktail->Init();\n");
+    case kSr90:  
+      fprintf(fp,"  AliGenRadioactive *pGen=new AliGenRadioactive(kSr90,%i);\n",(int)fGenNprimEntry->GetNumber());
+      fprintf(fp,"  pGen->SetOrigin(0,0,0);  pGen->SetSigma(0.1,0.1,0.05);\n");      
+      fprintf(fp,"  pGen->Init();\n");
+    break;  
+    case kHijing:  
+      fprintf(fp,"  AliGenHijing *pGen=new AliGenHijing(-1); pGen->SetEnergyCMS(5500); pGen->SetReferenceFrame(\"CMS\");\n");
+      fprintf(fp,"  pGen->SetProjectile(\"A\", 208, 82); pGen->SetTarget(\"A\", 208, 82);\n");      
+      fprintf(fp,"  pGen->SetJetQuenching(0); pGen->SetShadowing(0);\n");
+      fprintf(fp,"  pGen->KeepFullEvent(); pGen->SetSelectAll(0); \n");
+      fprintf(fp,"  pGen->SetImpactParameterRange(0., 5.);\n");
+      fprintf(fp,"  pGen->Init();\n");
     break;  
   }
 //central before RICH detectors                  
index d831c11..d0f46e4 100644 (file)
@@ -11,17 +11,20 @@ void Geom()
   Double_t dx,dy,dz,r1,r2,a=0,z=0,den=0,radlen=0,intlen=0;
   Double_t cm=1,m=100*cm,mm=0.1*cm,mkm=0.001*cm;  
 //Medium  
-                                                                                      pAir=new TGeoMedium("Air"      ,1, g->GetMaterial("Air"));
-  new TGeoMaterial("CH4"      ,a=26.98 ,z=13,den=0.4224                       );      pCH4=new TGeoMedium("CH4"      ,2, g->GetMaterial("CH4"));
-  new TGeoMaterial("CsI"      ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6);      pCsI=new TGeoMedium("CsI"      ,3, g->GetMaterial("CsI"));
-  new TGeoMaterial("Al"       ,a=26.98 ,z=13,den=2.7 ,radlen=24.01,intlen=70.6);      pAl =new TGeoMedium("Al"       ,4, g->GetMaterial("Al")); 
-  new TGeoMaterial("Fe"       ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8);      pFe =new TGeoMedium("Fe"       ,5, g->GetMaterial("Fe"));
-  new TGeoMaterial("W"        ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm);pW  =new TGeoMedium("W"        ,6, g->GetMaterial("W"));
-  new TGeoMaterial("Cu"       ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8);      pCu =new TGeoMedium("Cu"       ,7, g->GetMaterial("Cu"));
-  new TGeoMaterial("Perpex"   ,a=55.845,z=26,den=7.87,radlen=13.84,intlen=82.8);       new TGeoMedium("Perpex"   ,8, g->GetMaterial("Perpex")); 
-  new TGeoMaterial("Rohacell" ,a=12.01,z=6,den=0.1,radlen=18.8,intlen=0);              new TGeoMedium("Rohacell" ,9, g->GetMaterial("Rohacell"));
-  new TGeoMaterial("SiO2"     ,a=0,z=0,den=0);                                         new TGeoMedium("SiO2"     ,10, g->GetMaterial("SiO2"));
-  new TGeoMaterial("C6F14"    ,a=0,z=0,den=0);                                         new TGeoMedium("C6F14"    ,11, g->GetMaterial("C6F14"));
+                                                                                       pAir =new TGeoMedium("Air"      ,1, g->GetMaterial("Air"));
+  new TGeoMaterial("CH4"      ,a=26.98 ,z=13,den=0.4224                             ); pCH4 =new TGeoMedium("CH4"      ,2, g->GetMaterial("CH4"));
+  new TGeoMaterial("CsI"      ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pCsI =new TGeoMedium("CsI"      ,3, g->GetMaterial("CsI"));
+  new TGeoMaterial("Al"       ,a=26.98 ,z=13,den=2.7 ,radlen=24.01*cm,intlen=70.6*cm); pAl  =new TGeoMedium("Al"       ,4, g->GetMaterial("Al")); 
+  new TGeoMaterial("W"        ,a=183.84,z=27,den=19.3,radlen= 9.59*cm,intlen=0.35*cm); pW   =new TGeoMedium("W"      ,5, g->GetMaterial("W"));
+  new TGeoMaterial("Cu"       ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pCu  =new TGeoMedium("Cu"       ,6, g->GetMaterial("Cu"));
+  new TGeoMaterial("Rohacell",a=12.01 ,z=6 ,den=0.1,radlen=18.8,intlen=0);            pRoha =new TGeoMedium("Rohacell" ,7, g->GetMaterial("Rohacell"));
+  new TGeoMaterial("SiO2"     ,a=0     ,z=0 ,den=0);                                  pSiO2= new TGeoMedium("SiO2"     ,8,g->GetMaterial("SiO2"));
+  new TGeoMaterial("C6F14"    ,a=0     ,z=0 ,den=0);                                  pC6F14=new TGeoMedium("C6F14"    ,9,g->GetMaterial("C6F14"));
+//for Sr90 source  
+  new TGeoMaterial("Perpex"  ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pPerpex=new TGeoMedium("Perpex",10, g->GetMaterial("Perpex"));
+  new TGeoMaterial("Steel"   ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSteel=new TGeoMedium("Steel"  ,11, g->GetMaterial("Steel"));
+  new TGeoMaterial("Mylar"   ,a=55.845,z=26,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pMylar=new TGeoMedium("Mylar"  ,12, g->GetMaterial("Mylar"));
+  new TGeoMaterial("Sr"      ,a=87.62 ,z=38,den=7.87,radlen=13.84*cm,intlen=82.8*cm); pSr   =new TGeoMedium("Sr"     ,13, g->GetMaterial("Sr"));
 //Mother
   TGeoVolume *pMother=g->MakeBox("Mother",pAir,dx=10*m/2,dy=10*m/2,dz=30*m/2); //arbitrary values  
 //Rich  
@@ -43,7 +46,7 @@ void Geom()
   pRich->AddNode(pPpf,copy=5,new TGeoTranslation(-335*mm,+433*mm,8*cm+20*mm));
   pRich->AddNode(pPpf,copy=6,new TGeoTranslation(+335*mm,+433*mm,8*cm+20*mm));
   pPpf->AddNode( pPc ,copy=1,new TGeoTranslation(   0*mm,   0*mm,-19.15*mm));//PPF 2001P2 
-  pPpf->AddNode(g->GetVolume("PPFlarge"),copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm,  0.85*mm));
+  pPpf->AddNode(pPpfLarge,copy=1,new TGeoTranslation(-224.5*mm,-151.875*mm,  0.85*mm));
   pPpf->AddNode(g->GetVolume("PPFlarge"),copy=2,new TGeoTranslation(-224.5*mm,- 50.625*mm,  0.85*mm));
   pPpf->AddNode(g->GetVolume("PPFlarge"),copy=3,new TGeoTranslation(-224.5*mm,+ 50.625*mm,  0.85*mm));
   pPpf->AddNode(g->GetVolume("PPFlarge"),copy=4,new TGeoTranslation(-224.5*mm,+151.875*mm,  0.85*mm));
@@ -99,52 +102,49 @@ void Geom()
     TGeoVolume *pRadShort =g->MakeBox( "RadShort" ,g->GetMedium("Neoceram") ,dx=  10*mm/2 ,dy= 403*mm/2  ,dz=  15*mm/2);    
     TGeoVolume *pRadSpacer=g->MakeTube("RadSpacer",g->GetMedium("SiO2")     ,r1= 0        ,r2=10*mm/2  ,dz=  15*mm/2);         
     
-    pRich->AddNode(pRad,copy=1,new TGeoTranslation(   0*mm,-434*mm,   -12*mm)); 
-    pRich->AddNode(pRad,copy=2,new TGeoTranslation(   0*mm,   0*mm,   -12*mm)); 
-    pRich->AddNode(pRad,copy=3,new TGeoTranslation(   0*mm,+434*mm,   -12*mm)); 
-    pRad->AddNode(g->GetVolume("RadFront") ,copy=1,new TGeoTranslation(   0*mm,   0*mm, -10.0*mm));
-    pRad->AddNode(g->GetVolume("RadWin")   ,copy=1,new TGeoTranslation(   0*mm,   0*mm,   9.5*mm));
-    pRad->AddNode(g->GetVolume("RadLong")  ,copy=1,new TGeoTranslation(   0*mm,-204*mm,  -0.5*mm));
-    pRad->AddNode(g->GetVolume("RadLong")  ,copy=2,new TGeoTranslation(   0*mm,+204*mm,  -0.5*mm));
-    pRad->AddNode(g->GetVolume("RadShort") ,copy=1,new TGeoTranslation(-660*mm,   0*mm,  -0.5*mm));
-    pRad->AddNode(g->GetVolume("RadShort") ,copy=2,new TGeoTranslation(+660*mm,   0*mm,  -0.5*mm));
+    pRich->AddNode(pRad      ,copy=1,new TGeoTranslation(   0*mm,-434*mm,   -12*mm)); 
+    pRich->AddNode(pRad      ,copy=2,new TGeoTranslation(   0*mm,   0*mm,   -12*mm)); 
+    pRich->AddNode(pRad      ,copy=3,new TGeoTranslation(   0*mm,+434*mm,   -12*mm)); 
+    
+    pRad ->AddNode(pRadFront ,copy=1,new TGeoTranslation(   0*mm,   0*mm, -10.0*mm));
+    pRad ->AddNode(pRadWin   ,copy=1,new TGeoTranslation(   0*mm,   0*mm,   9.5*mm));
+    pRad ->AddNode(pRadLong  ,copy=1,new TGeoTranslation(   0*mm,-204*mm,  -0.5*mm));
+    pRad ->AddNode(pRadLong  ,copy=2,new TGeoTranslation(   0*mm,+204*mm,  -0.5*mm));
+    pRad ->AddNode(pRadShort ,copy=1,new TGeoTranslation(-660*mm,   0*mm,  -0.5*mm));
+    pRad ->AddNode(pRadShort ,copy=2,new TGeoTranslation(+660*mm,   0*mm,  -0.5*mm));
     for(int i=0;i<3;i++)
       for(int j=0;j<10;j++)
         pRad->AddNode(pRadSpacer,copy=10*i+j,new TGeoTranslation(-1330*mm/2+116*mm+j*122*mm,(i-1)*105*mm,-0.5*mm));
   }else{//Sr90
-    Double_t containerZ=21.8*mm/2,containerR=38*mm/2;  
-    Double_t screwLen=15*mm/2      ,screwR=5*mm/2;  
-    Double_t srcLen=10*mm/2        ,srcR=2*mm/2;  
-    Double_t perpexZ=20*mm/2       ,perpexR=34*mm/2;  
-    Double_t perpexWholeLen=10*mm/2,perpexWholeR=4*mm/2;  
-    Double_t alBottomLen=containerZ-perpexZ,alWholeR=5*mm/2;  
-    g->MakeTube("Sr90"           ,pAir   ,0      ,containerR  ,containerZ);             
-    g->MakeTube("Sr90AlWall"     ,pAl    ,perpexR,containerR  ,containerZ);             
-    g->MakeTube("Sr90AlBottom"   ,pAl    ,0      ,perpexR     ,alBottomLen);            
-    g->MakeTube("Sr90AlWhole"    ,pAir   ,0      ,alWholeR    ,alBottomLen);            
-    g->MakeTube("Sr90PerpexPlug" ,g->GetMedium("Perpex"),0      ,perpexR     ,perpexZ);                
-    g->MakeTube("Sr90PerpexWhole",pAir   ,0      ,perpexWholeR,perpexWholeLen);         
-    g->MakeTube("Sr90Screw"      ,g->GetMedium("Steel") ,0      ,screwR      ,screwLen);               
-    g->MakeTube("Sr90Source"     ,g->GetMedium("Steel") ,0      ,srcR        ,srcLen);     
+    pSrc               =g->MakeTube("Src"              ,pCH4    , 0 , 70*mm/2 ,  30*mm/2);       //top container
+      pAlGlass         =g->MakeTube("SrcAlGlass"        ,pAl    , 0 , 38*mm/2 ,21.8*mm/2);       //Al glass wall        
+        pPerpexPlug    =g->MakeTube("SrcPerpex"        ,pPerpex , 0 , 34*mm/2 ,  20*mm/2);       //Perpex plug         
+          pSteelScrew  =g->MakeTube("SrcSteelScrew"    ,pSteel  , 0 ,  5*mm/2 ,  15*mm/2);       //Steel screw in the center        
+          pSr90Screw   =g->MakeTube("SrcSr90Screw"     ,pSteel  , 0 ,  2*mm/2 ,  10*mm/2);       //Steel screw to support Sr90 
+            pSr90      =g->MakeTube("SrcSr90"          ,pSr     , 0 ,  1*mm/2 ,   1*mm/2);       //Sr90 source
+          pHolePerpex  =g->MakeTube("SrcHolePerpex"    ,pAir    , 0 ,  4*mm/2 ,  10*mm/2);       //Air hole in perpex plug      
+        pHoleAl        =g->MakeTube("SrcHoleAl"        ,pAir    , 0 ,  4*mm/2 , 1.8*mm/2);       //Air hole in Al glass bottom
+    pMylarFoil         =g->MakeTube("SrcMylarFoil"     ,pMylar  , 0 , 30*mm/2 , 50*mkm/2);       //Mylar foil                
                 
-    pRich->AddNode(pContainer,1,new TGeoTranslation(30*cm,0,containerZ));
-    pContainer ->AddNode(pAlWall,1);
-    pContainer ->AddNode(pAlBottom,1,new TGeoTranslation(0,0,-containerZ+alBottomLen));
-    pContainer ->AddNode(pPerpexPlug,1,new TGeoTranslation(0,0,containerZ-perpexZ));
-    pAlBottom  ->AddNode(pAlWhole,1,new TGeoTranslation(6*mm,0,0));
-    pPerpexPlug->AddNode(pPerpexWhole,1,new TGeoTranslation(6*mm,0,-perpexZ+perpexWholeLen));  
-    pPerpexPlug->AddNode(pSource,1     ,new TGeoTranslation(6*mm,0, perpexZ-srcLen));  
-    pPerpexPlug->AddNode(pScrew,1,new TGeoTranslation(0,0,perpexZ-screwLen));  
+    pRich->AddNode(pSrc,1,new TGeoTranslation(30*cm,0,alVesselZ));
+      pSrc ->AddNode(pMylarFoil,1,new TGeoTranslation(0,0,21.8*mm/2+50*mkm/2));
+      pSrc ->AddNode(pAlGlass,1,new TGeoTranslation(0,0,0));//Al glass to fake Src volume
+        pAlGlass->AddNode(pHoleAl,1,new TGeoTranslation(6*mm, 0*mm, 10.9*mm));
+        pAlGlass->AddNode(pPerpexPlug,1,new TGeoTranslation(0*mm,0*mm,-0.9*mm));
+          pPerpexPlug->AddNode(pSteelScrew,1,new TGeoTranslation(shift,0, perpexPlugZ-steelScrewZ));  
+          pPerpexPlug->AddNode(pHolePerpex,1,new TGeoTranslation(shift,0,-perpexPlugZ+airWholeZ));      
+          pPerpexPlug->AddNode(pSr90Screw ,1,new TGeoTranslation(shift,0, perpexPlugZ-sr90ScrewZ));  
+          pSr90Screw->AddNode(pSr90 ,1,new TGeoTranslation(0,0, -(sr90ScrewZ-sr90Z)));  
   }//if(!isSr90)
 //Sandbox  
-  TGeoVolume *pSandBox=g->MakeBox( "SandBox"  ,g->GetMedium("Air")      ,dx=1419*mm/2 ,dy=1378*mm/2   ,dz=50.5*mm/2);  //2072P1   
-  g->MakeBox( "SandCover",g->GetMedium("Al")       ,dx=1419*mm/2 ,dy=1378*mm/2   ,dz= 0.5*mm/2);  
-  g->MakeBox( "SandComb" ,g->GetMedium("Rohacell") ,dx=1359*mm/2 ,dy=1318*mm/2   ,dz=49.5*mm/2);  
+  TGeoVolume *pSandBox  =g->MakeBox( "SandBox"  ,pAir  ,dx=1419*mm/2 ,dy=1378*mm/2   ,dz=50.5*mm/2);  //2072P1   
+  TGeoVolume *pSandCover=g->MakeBox( "SandCover",pAl   ,dx=1419*mm/2 ,dy=1378*mm/2   ,dz= 0.5*mm/2);  
+  TGeoVolume *pSandComb =g->MakeBox( "SandComb" ,pRoha ,dx=1359*mm/2 ,dy=1318*mm/2   ,dz=49.5*mm/2);  
   
   pRich->AddNode(pSandBox,copy=1,new TGeoTranslation(  0*mm,0*mm, -73.75*mm)); //p.84 TDR
-  pSandBox->AddNode(g->GetVolume("SandComb")  ,copy=1,new TGeoTranslation(  0*mm,0*mm,      0*mm)); //2072P1
-  pSandBox->AddNode(g->GetVolume("SandCover") ,copy=1,new TGeoTranslation(  0*mm,0*mm,    +25*mm)); 
-  pSandBox->AddNode(g->GetVolume("SandCover") ,copy=2,new TGeoTranslation(  0*mm,0*mm,    -25*mm)); 
+    pSandBox->AddNode(pSandComb  ,copy=1,new TGeoTranslation(  0*mm,0*mm,      0*mm)); //2072P1
+    pSandBox->AddNode(pSandCover ,copy=1,new TGeoTranslation(  0*mm,0*mm,    +25*mm)); 
+    pSandBox->AddNode(pSandCover ,copy=2,new TGeoTranslation(  0*mm,0*mm,    -25*mm)); 
 //Colors  
   g->GetVolume("Pc")       ->SetLineColor(kGreen);
   g->GetVolume("PPF")      ->SetLineColor(33);
index a5143b5..669fad6 100644 (file)
@@ -15,12 +15,11 @@ void Opticals()
   Float_t aIdxSiO2[kNbins];
   Float_t aIdxOpSiO2[kNbins];
         
-  Float_t  e1=10.666,e2=18.125,f1=46.411,f2= 228.71; //RICH TDR page 35 
   for (i=0;i<kNbins;i++){
-    aIdxC6F14[i] = aPckov[i] *0.0172*1e9+1.177;
-    aIdxSiO2[i]  = TMath::Sqrt(1.+f1/(e1*e1-TMath::Power(aPckov[i]*1e9,2))+f2/(e2*e2-TMath::Power(aPckov[i]*1e9,2)));//TDR p.35
-    aIdxCH4[i]   =1.000444;
-    aIdxGrid[i]  =1;
+    aIdxC6F14[i] = (Float_t)AliRICHParam::IndOfRefC6F14(aPckov[i]*1e9);
+    aIdxSiO2[i]  = (Float_t)AliRICHParam::IndOfRefSiO2(aPckov[i]*1e9);
+    aIdxCH4[i]   = (Float_t)AliRICHParam::IndOfRefCH4();
+    aIdxGrid[i]    =1;
     aIdxOpSiO2[i]  =1;
   } 
     
index 87321e3..aaeed71 100644 (file)
@@ -16,4 +16,8 @@
 #pragma link C++ class  AliRICHDigitizer+;
 #pragma link C++ class  AliRICHDisplFast+;
 #pragma link C++ class  AliRICHReconstructor+;
+#pragma link C++ class  AliRICHHelix+;
+#pragma link C++ class  AliGenRadioactive+;
+#pragma link C++ class  AliRICHTracker+;
+#pragma link C++ enum   RadioNuclei_t;
 #endif
index 364e726..b7c808c 100644 (file)
@@ -1,20 +1,28 @@
 void RichBatch(const Int_t iNevents,const Bool_t isDebug,const char *sConfigFileName)
 {
-  if(isDebug) gAlice->SetDebug(1);
+  gSystem->Exec("rm -rf *.root hlt hough ZZZ*");
+  if(isDebug)   AliLog::SetGlobalDebugLevel(AliLog::kDebug);
 
-  Info("my/RichBatch.C","%i event(s) requested, debug %i,config file %s",iNevents,isDebug,sConfigFileName);  
   TStopwatch sw;TDatime time;  
-
+  
+  
   AliSimulation     *pSim=new AliSimulation;     pSim->Run(iNevents); delete pSim;
-  AliReconstruction *pRec=new AliReconstruction; pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH");  pRec->Run();         delete pRec;
+  AliReconstruction *pRec=new AliReconstruction; 
+  pRec->SetRunLocalReconstruction("ITS TPC TRD TOF RICH");
+  pRec->SetFillESD("ITS TPC TRD TOF RICH");
+  pRec->Run();         delete pRec;
+  
+  //pRec->SetRunLocalReconstruction("RICH"); pRec->SetFillESD("RICH");    
+  
+  cout<<"\n!!!!!!!!!!!!Info in <my/RichBatch.C>: Start creating Control Plots \n";
   
-  AliRunLoader* pAL = AliRunLoader::Open();
-  pAL->LoadgAlice();
+  AliRunLoader* pAL = AliRunLoader::Open();  pAL->LoadgAlice();
   AliRICH *pRICH=(AliRICH*)pAL->GetAliRun()->GetDetector("RICH");
   if(pRICH) pRICH->ControlPlots();
  
-  cout<<"\nInfo in <my/RichBatch.C>: Start time: ";time.Print();
-    cout<<"Info in <my/RichBatch.C>: Stop  time: ";time.Set();  time.Print();
-    cout<<"Info in <my/RichBatch.C>: Time  used: ";sw.Print();
+  cout<<"\n!!!!!!!!!!!!Info in <my/RichBatch.C>: Start time: ";time.Print();
+    cout<<"!!!!!!!!!!!!Info in <my/RichBatch.C>: Stop  time: ";time.Set();  time.Print();
+    cout<<"!!!!!!!!!!!!Info in <my/RichBatch.C>: Time  used: ";sw.Print();
   gSystem->Exec("touch ZZZ______finished_______ZZZ");
+//  gSystem->Exec("playwave .kde/my/end.wav");
 }
index 95c4933..0e059c7 100644 (file)
@@ -3,7 +3,7 @@ How to open session?
 How to retrive pointer to alice run loader:
         use pRICH->GetLoader()->GetRunLoader() (all detector classes inherit from AliDetector wich has GetLoader())
        use methode AliRun::GetRunLoader for gAlice (depricated)
-How to get pointers to deifferent root trees:
+How to get pointers to different root trees:
        TreeE belongs to AliRunLoader, available after AliRunLoader::LoadHeader()
        TreeK belongs to AliRunLoader, available after AliRunLoader::LoadKinematics()
        TreeH belongs to AliLoader   , available after AliLoader::LoadHits()
@@ -13,19 +13,20 @@ How to get pointers to deifferent root trees:
 
 
 How to work with the stack of particles?
-       - pointer to the stack is returned by gAlice->Stack() (global gAlice of type AliRun) or AliRunLoader::Stack() but
-        before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0.
-        Moreover loading header gives the information about number of particles only. 
-        To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics()        
-       - total amount of particles in stack for a given event:         AliRunLoader::Stack()->GetNtrack() or AliRun::GetEvent() (after LoadHeader())
-       - total amount of primiry particles in stack for a given event: AliRunLoader::Stack()->GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader())
+        - first of all, the stack includes primary as well as secondary particles
+       - pointer to the stack is returned by AliRun::Stack() (global gAlice of type AliRun - depricated way to do)
+          or by AliRunLoader::Stack() but before one needs to load event header by AliRunLoader::LoadHeader() otherwise both methods return 0.
+          Moreover loading header gives the information about number of particles only. 
+          To retrive the list of particle one also needs to load kinematics by AliRunLoader::LoadKinematics()        
+       - total amount of particles in stack for a given event:         AliStack::GetNtrack() or AliRun::GetEvent() (after LoadHeader())
+       - total amount of primiry particles in stack for a given event: AliStack::GetNprimary() or AliLoader::TreeH()->GetEntries() (after LoadHeader())
 
 
 
 How to retrive hits:
-       Hits a stored on primiry by primiry basis. To retrieve all hits one needs to do:
+       Hits a stored on primary by primary basis. To retrieve all hits one needs to do:
        initialise the root tree and containers:    AliLoader::LoadHits() 
-       read number of primiries in current event:
+       read number of primaries in current event:
        loop on the list of primiries:
 
 
@@ -36,13 +37,8 @@ How to retrive sdigits?
        One needs to say:
        pRich->GetLoader()->LoadSDigits(); this one open file, get the tree and invoke AliRICH::SetTreeAddress()    
 
-
-
-gAlice->GetMCApp()->GetCurrentTrackNumber()
-
-
 What are the debug methodes avail:
-       AliModule::GetDebug() 
+       AliModule::GetDebug() now depricated, use AliDebug printout instead
        AliModule::SetDebug()
        AliRun::GetDebug()
        AliRun::SetDebug()
@@ -61,5 +57,3 @@ How to deal with AliRunDigitizer?
 How to avoid using gAlice?
        Rich()->GetLoader()->GetRunLoader()->GetAliRun() returns gAlice global pointer.
          
-
-
index b3a7f34..4b1a095 100644 (file)
@@ -3,7 +3,7 @@ SRCS   =  AliRICH.cxx AliRICHv0.cxx AliRICHv1.cxx\
         AliRICHClusterFinder.cxx AliRICHMap.cxx\
         AliRICHChamber.cxx AliRICHRecon.cxx\
          AliRICHDisplFast.cxx\
-        AliRICHDigitizer.cxx AliRICHReconstructor.cxx
+        AliRICHDigitizer.cxx AliGenRadioactive.cxx AliRICHHelix.cxx AliRICHTracker.cxx AliRICHReconstructor.cxx
 
 HDRS =  $(SRCS:.cxx=.h)
 DHDR= RICHLinkDef.h
index c312014..10e3d18 100644 (file)
@@ -3,6 +3,7 @@ void ph(Int_t event=0) {R()->PrintHits(event);}    //utility print hits for 'eve
 void ps(Int_t event=0) {R()->PrintSDigits(event);} //utility print sdigits
 void pd(Int_t event=0) {R()->PrintDigits(event);}  //utility print digits
 void pc(Int_t event=0) {R()->PrintClusters(event);}//utility print clusters
+void pt(Int_t event=0) {R()->PrintTracks(event);}  //utility print tracks
 
 //__________________________________________________________________________________________________
 void pp(int tid)
@@ -19,6 +20,7 @@ void pp(int tid)
 //__________________________________________________________________________________________________
 void PrintParticleInfo(int tid)
 {
+// Prints particle info for a given TID
   TParticle *p=al->Stack()->Particle(tid);
   cout<<p->GetName()<<"("<<tid<<")";
   if(p->GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));}
@@ -27,6 +29,7 @@ void PrintParticleInfo(int tid)
 //__________________________________________________________________________________________________
 Int_t prim(Int_t tid)
 {
+// Provides mother TID for the given TID
   R()->GetLoader()->GetRunLoader()->LoadHeader();  R()->GetLoader()->GetRunLoader()->LoadKinematics();
   
   if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack())
@@ -147,13 +150,11 @@ AliRICH *r;
 Bool_t ReadAlice()
 {
   Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
-  AliLoader::SetDebug(0);
   if(gAlice) delete gAlice;      
   if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){
     gSystem->Exec("rm -rf *.root *.dat");
     Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one");
     new AliRun("gAlice","Alice experiment system");
-    gAlice->SetDebug(-1);
     gAlice->Init("Config.C");
     r=(AliRICH*)gAlice->GetDetector("RICH");
     return kFALSE;
@@ -164,7 +165,6 @@ Bool_t ReadAlice()
   a->SetDebug(0);    
 //RICH      
   if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
-  r->SetDebug(0);
   if(!(rl=al->GetLoader("RICHLoader")))          Warning("RICH/menu.C::ReadAlice","No RICH loader in file");        
         
   Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());      
@@ -353,22 +353,22 @@ void TestSeg()
   TPolyLine3D *pXaxis=new TPolyLine3D(2,X);pXaxis->SetLineColor(kRed);   pXaxis->Draw();
   TPolyLine3D *pYaxis=new TPolyLine3D(2,Y);pYaxis->SetLineColor(kGreen); pYaxis->Draw();
   TPolyLine3D *pZaxis=new TPolyLine3D(2,Z);pZaxis->SetLineColor(kBlue);  pZaxis->Draw();  
-//chambers  
+//Draw PC for all chambers by trasfering Pc plane using Pc2Mrs methode  
   Int_t iNpointsX=50,iNpointsY=50;  
   for(Int_t iChamberN=1;iChamberN<=7;iChamberN++){//chamber loop
     TPolyMarker3D *pChamber=new TPolyMarker3D(iNpointsX*iNpointsY);
     Int_t i=0;
     for(Double_t x=0;x<p.PcSizeX();x+=p.PcSizeX()/iNpointsX)
       for(Double_t y=0;y<p.PcSizeY();y+=p.PcSizeY()/iNpointsY){//step loop
-        TVector3 v3=p.C(iChamberN)->Loc2Glob(TVector2(x,y));
-        pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());
-        pChamber->SetMarkerSize(1);
-        pChamber->SetMarkerColor(iChamberN);
+        TVector3 v3=p.C(iChamberN)->Pc2Mrs(TVector2(x,y));//from regular grid of local PC points to MRS presentation
+        pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());//Pc plane poing in MRS
       }//step loop
+    pChamber->SetMarkerSize(1);
+    pChamber->SetMarkerColor(iChamberN);
     pChamber->Draw();  
     t.SetNDC();t.SetTextColor(iChamberN); t.DrawText(0.1,iChamberN*0.1,Form("Chamber %i",iChamberN));    
   }//chamber loop   
-  gPad->GetView()->RotateView(94,45);
+//  gPad->GetView()->RotateView(94,45);
 }//void TestSeg()
 //__________________________________________________________________________________________________
 void TestMenu()
@@ -423,3 +423,92 @@ void ReadRaw()
   }
   fclose(pRawFile);
 }
+//__________________________________________________________________________________________________
+void FillContribs(Int_t part,Int_t C, Bool_t print)
+{
+  static Int_t contrib[10][7][2];
+  static Bool_t kEnter=kFALSE;
+
+  C--; // chamber numbering from 1 to 7
+  if(!kEnter) {
+    for(Int_t i=0;i<10;i++) {
+      for(Int_t j=0;j<7;j++) {
+        for(Int_t k=0;k<2;k++) contrib[i][j][k]=0;
+      }
+    }
+    kEnter=kTRUE;
+  }
+
+  if(print) {
+    for(Int_t k=0;k<2;k++) {cout << "----------------Positives  to RICH ---------------" << endl;
+                            cout << " chamber    1    2    3     4     5     6    7    " << endl;
+                            cout << " -------------------------------------------------" << endl;
+      for(Int_t i=0;i<4;i++) {
+        if(i==0) cout << " e";
+        if(i==1) cout << "pi";
+        if(i==2) cout << " K";
+        if(i==3) cout << " p";
+        if(k==0) cout << "+         ";
+        if(k==1) cout << "-         ";
+        for(Int_t j=0;j<7;j++) {
+          cout << contrib[i][j][k] << "    ";
+        }
+          cout << endl;
+      }
+    }
+  }
+
+  // +ves
+  if(part==kPositron)    contrib[0][C][0]++;
+  if(part==kPiPlus)      contrib[1][C][0]++;
+  if(part==kKPlus)       contrib[2][C][0]++;
+  if(part==kProton)      contrib[3][C][0]++;
+
+  // -ves
+  if(part==kElectron)    contrib[0][C][1]++;
+  if(part==kPiMinus)     contrib[1][C][1]++;
+  if(part==kKMinus)      contrib[2][C][1]++;
+  if(part==kProtonBar)   contrib[3][C][1]++;
+}
+//__________________________________________________________________________________________________
+void ParticleContribs()
+{
+
+  TH1F *distH1 = new TH1F("distH1","distH1",100,0.,5.);
+  AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH");
+  Bool_t isHits    =!pRich->GetLoader()->LoadHits();
+
+  pRich->GetLoader()->GetRunLoader()->LoadHeader();  pRich->GetLoader()->GetRunLoader()->LoadKinematics();
+
+
+  if(!isHits){Error("Exec","No hits. No contribs to calculate.");return;}
+
+  for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events Loop
+    pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN);
+    cout << " event n. " << iEventN << endl;
+    Int_t nPrimaries = (Int_t)pRich->GetLoader()->TreeH()->GetEntries();
+    TObjArray * Hits = new TObjArray[nPrimaries];
+
+    for(Int_t i=0;i<nPrimaries;i++) {
+      pRich->GetLoader()->TreeH()->GetEntry(i);
+      Int_t nHits = pRich->Hits()->GetEntries();
+      for(Int_t k=0;k<nHits;k++)         Hits[i].Add(pRich->Hits()->At(k));
+
+    }
+//deals with hits
+      for(Int_t i=0;i<nPrimaries;i++){//prims loop
+        pRich->GetLoader()->TreeH()->GetEntry(i);
+        Int_t nHits = pRich->Hits()->GetEntries();
+        for(Int_t j=0;j<nHits;j++){//hits loop
+          AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j);
+          TParticle *pParticle = pRich->GetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(pHit->GetTrack());
+//          if(pParticle->P()>1) FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
+          FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
+          if(pParticle->GetPDG()->Charge()!=0) distH1->Fill(pHit->Length());
+        }//hits loop
+      }//prims loop
+    }// event loop
+    FillContribs(0,0,kTRUE);
+    distH1->Draw();
+
+}//ParticleContribs()