New convention on numbering and LRS
authorkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 May 2004 14:51:30 +0000 (14:51 +0000)
committerkir <kir@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 31 May 2004 14:51:30 +0000 (14:51 +0000)
17 files changed:
RICH/AliRICH.cxx
RICH/AliRICH.h
RICH/AliRICHChamber.cxx
RICH/AliRICHChamber.h
RICH/AliRICHClusterFinder.cxx
RICH/AliRICHClusterFinder.h
RICH/AliRICHDigitizer.cxx
RICH/AliRICHDigitizer.h
RICH/AliRICHDisplFast.cxx
RICH/AliRICHDisplFast.h
RICH/AliRICHParam.cxx
RICH/AliRICHParam.h
RICH/AliRICHv0.cxx
RICH/AliRICHv0.h
RICH/AliRICHv1.cxx
RICH/AliRICHv1.h
RICH/RICHLinkDef.h

index 2a8acba..d57feb0 100644 (file)
@@ -16,6 +16,7 @@
 #include "AliRICH.h"
 #include "AliRICHParam.h"
 #include "AliRICHChamber.h"
+#include "AliRICHClusterFinder.h"
 #include <TArrayF.h>
 #include <TGeometry.h>
 #include <TBRIK.h>
 #include <AliRun.h>
 #include <AliRunDigitizer.h>
 #include <AliMC.h>
+#include <AliESD.h>
 #include <TVirtualMC.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TStopwatch.h>
  
 ClassImp(AliRICHhit)
 //__________________________________________________________________________________________________
@@ -43,23 +48,27 @@ ClassImp(AliRICHdigit)
 //__________________________________________________________________________________________________
 void AliRICHdigit::Print(Option_t*)const
 {
-  ::Info("digit","csxy=%6i, cfm=%9i, c=%2i, x=%3i, y=%3i, q=%8.3f, TID1=%5i, TID2=%5i, TID3=%5i",
-                  Id(),fChFbMip,fChamber,fPadX,fPadY,fQdc,fTracks[0],fTracks[1],fTracks[2]);
+  ::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]);
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICHcluster)
 //__________________________________________________________________________________________________
 void AliRICHcluster::Print(Option_t*)const
 {
-  ::Info("cluster","CombiPid=%10i, c=%2i, size=%6i, dim=%5i, x=%7.3f, y=%7.3f, Q=%6i, st=%i",
-           fCombiPid,fChamber,fSize,fDimXY,fX,fY,fQdc,fStatus);
-}
-//__________________________________________________________________________________________________
-ClassImp(AliRICHreco)
-//__________________________________________________________________________________________________
-void AliRICHreco::Print(Option_t*)const
-{
-  ::Info("reco","ThetaCherenkov=%9.6f, Nphotons=%4i, TID=%9i",fThetaCherenkov,fNphotons,fTid);
+  char *status=0;
+  switch(fStatus){
+    case      kRaw: status="raw"     ;break;
+    case kResolved: status="resolved";break;
+    case    kEmpty: status="empty"   ;break;
+  }
+  if(fDigits)    
+    ::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);
+    
 }
 //__________________________________________________________________________________________________
 ClassImp(AliRICH)    
@@ -71,32 +80,23 @@ ClassImp(AliRICH)
 */
 //END_HTML
 //__________________________________________________________________________________________________
-AliRICH::AliRICH()
-        :AliDetector() 
+AliRICH::AliRICH():AliDetector(),fpParam(0),  fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0) 
 {
 //Default ctor should not contain any new operators
-  fpParam     =0;
-  fChambers   =0;   
 //AliDetector ctor deals with Hits and Digits  
-  fSdigits    =0; fNsdigits   =0;
-  fDigitsNew  =0; for(int i=0;i<kNCH;i++) fNdigitsNew[i]  =0;
-  fClusters   =0; for(int i=0;i<kNCH;i++) fNclusters[i]=0;
-  fRecos      =0; fNrecos     =0;
+  for(int i=0;i<kNchambers;i++) fNdigitsNew[i]  =0;
+  for(int i=0;i<kNchambers;i++) fNclusters[i]=0;
+//  fCounters.ResizeTo(20); fCounters.Zero();
 }//AliRICH::AliRICH()
 //__________________________________________________________________________________________________
 AliRICH::AliRICH(const char *name, const char *title)
-        :AliDetector(name,title)
+        :AliDetector(name,title),fpParam(new AliRICHParam),fSdigits(0),fNsdigits(0),fDigitsNew(0),fClusters(0)
 {
 //Named ctor
   if(GetDebug())Info("named ctor","Start.");
-  fpParam     =   new AliRICHParam;
-  fChambers = 0;  CreateChambers();
 //AliDetector ctor deals with Hits and Digits (reset them to 0, does not create them)
-  fHits=       0;     CreateHits();          gAlice->GetMCApp()->AddHitList(fHits);
-  fSdigits=    0;
-  fDigitsNew=  0;
-  fClusters=   0;
-  fRecos      =0;
+  CreateHits();          gAlice->GetMCApp()->AddHitList(fHits);
+  fCounters.ResizeTo(20); fCounters.Zero();
   if(GetDebug())Info("named ctor","Stop.");
 }//AliRICH::AliRICH(const char *name, const char *title)
 //__________________________________________________________________________________________________
@@ -106,59 +106,49 @@ AliRICH::~AliRICH()
   if(GetDebug()) Info("dtor","Start.");
 
   if(fpParam)    delete fpParam;
-  if(fChambers)  delete fChambers;
   
   if(fHits)      delete fHits;
   if(fSdigits)   delete fSdigits;
   if(fDigits)    delete fDigits;
   if(fDigitsNew) {fDigitsNew->Delete();   delete fDigitsNew;}
   if(fClusters)  {fClusters->Delete();    delete fClusters;}
-  if(fRecos)     delete fRecos;
   if(GetDebug()) Info("dtor","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.");
-
-  AliLoader * richLoader = GetLoader();
-  AliRunLoader * runLoader = GetLoader()->GetRunLoader();
-
   for(Int_t iEventN=0;iEventN<GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEventN++){//events loop
-    runLoader->GetEvent(iEventN);
+    GetLoader()->GetRunLoader()->GetEvent(iEventN);//get next event
   
-    if (!richLoader->TreeH()) richLoader->LoadHits();
-    if (!runLoader->TreeE()) runLoader->LoadHeader(); 
-    if (!runLoader->TreeK()) runLoader->LoadKinematics();//from
-    if (!richLoader->TreeS()) richLoader->MakeTree("S"); MakeBranch("S");//to
+    if(!GetLoader()->TreeH()) GetLoader()->LoadHits();    GetLoader()->GetRunLoader()->LoadHeader(); 
+                                                          GetLoader()->GetRunLoader()->LoadKinematics();//from
+    if(!GetLoader()->TreeS()) GetLoader()->MakeTree("S"); MakeBranch("S");//to
           
     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
-      richLoader->TreeH()->GetEntry(iPrimN);
+      GetLoader()->TreeH()->GetEntry(iPrimN);
       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){//hits loop 
-        AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);                
-        TVector2 x2 = P()->ShiftToWirePos(C(pHit->C())->Glob2Loc(pHit->OutX3()));                
-        Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());
+        AliRICHhit *pHit=(AliRICHhit*)Hits()->At(iHitN);//get current hit                
+        TVector2 x2 = C(pHit->C())->Glob2Loc(pHit->OutX3());//hit position in the chamber local system
+        Int_t iTotQdc=P()->TotQdc(x2,pHit->Eloss());//total charge produced by hit, 0 if hit in dead zone
         if(iTotQdc==0) continue;
-        Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax;
-        P()->Loc2Area(x2,iPadXmin,iPadYmin,iPadXmax,iPadYmax);//determine affected pads
-        if(GetDebug()) Info("Hits2SDigits","left-down=(%i,%i) right-up=(%i,%i)",iPadXmin,iPadYmin,iPadXmax,iPadYmax);
-        for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)//affected pads loop
-          for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++){
-            Double_t padQdc=iTotQdc*P()->FracQdc(x2,iPadX,iPadY);
-            if(padQdc>0.1) AddSDigit(pHit->C(),iPadX,iPadY,padQdc,
-              runLoader->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
+        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);
+        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);
+            if(padQdc>0.1) AddSDigit(pHit->C(),pad,padQdc,GetLoader()->GetRunLoader()->Stack()->Particle(pHit->GetTrack())->GetPdgCode(),pHit->GetTrack());
           }//affected pads loop 
       }//hits loop
     }//prims loop
-    richLoader->TreeS()->Fill();
-    richLoader->WriteSDigits("OVERWRITE");
+    GetLoader()->TreeS()->Fill();
+    GetLoader()->WriteSDigits("OVERWRITE");
     ResetSDigits();
   }//events loop  
-  richLoader->UnloadHits();
-  runLoader->UnloadHeader();
-  runLoader->UnloadKinematics();
+  GetLoader()->UnloadHits(); GetLoader()->GetRunLoader()->UnloadHeader(); GetLoader()->GetRunLoader()->UnloadKinematics();
   GetLoader()->UnloadSDigits();  
   if(GetDebug()) Info("Hit2SDigits","Stop.");
 }//Hits2SDigits()
@@ -177,7 +167,7 @@ void AliRICH::BuildGeometry()
   Float_t len=P()->SectorSizeY();
   new TBRIK("PHOTO","PHOTO","void",wid/2,0.1,len/2);
   
-  for(int i=1;i<=kNCH;i++){
+  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->SetLineColor(kRed);
@@ -208,18 +198,17 @@ void AliRICH::BuildGeometry()
 //______________________________________________________________________________
 void AliRICH::CreateMaterials()
 {
-    //
-    // *** DEFINITION OF AVAILABLE RICH MATERIALS *** 
-  
+// Definition of available RICH materials  
 #include "Opticals.h"
         
   Float_t a=0,z=0,den=0,radl=0,absl=0;
   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( 1, "Air     $",a=14.61,z=7.3, den=0.001205,radl=30420.0,absl=67500);//(Air)
-  AliMedium(1, "DEFAULT MEDIUM AIR$", 1, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  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);
   
   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);
@@ -230,11 +219,11 @@ void AliRICH::CreateMaterials()
   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);
   
-  AliMaterial(50, "ALUM",     a=26.98,z=13.0,den=2.7,     radl=8.9,    absl=0);    //aluminium sheet (Al)
+  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(31, "COPPER$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);    //(Cu)
-  AliMedium(12, "PCB_COPPER", 31, 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);
   
   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 
@@ -244,21 +233,30 @@ void AliRICH::CreateMaterials()
   AliMedium(8, "QUARTZO$", 21, 1, 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 (30, "FRE",aFreon,zFreon,den=1.7,-2,wmatFreon);//Freon (C6F14) 
-  AliMedium(4, "FREON$", 30, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  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);
   
   Float_t aMethane[2]={12.01,1}; Float_t zMethane[2]={6,1}; Float_t wmatMethane[2]={1,4};
-  AliMixture (40, "MET", aMethane, zMethane, den=7.17e-4,-2, wmatMethane);//methane (CH4)     
-  AliMedium(5, "METHANE$", 40, 1, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  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
+    AliMedium(kSteel, "STEEL", material, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
   
-  AliMixture (41, "METG", aMethane, zMethane, den=7.17e-4, -2, wmatMethane);
-  AliMedium(kGAP, "GAP$", 41, 1, isxfld, sxmgmx,tmaxfd, 0.1, -deemax, epsil, -stmin);
+    AliMaterial(material=46, "PERPEX$",  a=63.54,z=29.0,den=8.96,    radl=1.4,    absl=0);    //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
+    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 (32, "GLASS",aGlass, zGlass, den=1.74, 5, wGlass);//Glass 50%C+10.5%Si+35.5%O+3% + 1%
-  AliMedium(11, "GLASS", 32, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin);
+  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);
@@ -387,14 +385,14 @@ void AliRICH::MakeBranch(Option_t* option)
    
   if(cD&&fLoader->TreeD()){//D
     CreateDigits();
-    for(Int_t i=0;i<kNCH;i++){ 
+    for(Int_t i=0;i<kNchambers;i++){ 
       MakeBranchInTree(fLoader->TreeD(),Form("%s%d",GetName(),i+1),&((*fDigitsNew)[i]),kBufferSize,0);
     }
   }//D
   
   if(cR&&fLoader->TreeR()){//R
     CreateClusters();
-    for(Int_t i=0;i<kNCH;i++)
+    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.");   
@@ -420,7 +418,7 @@ void AliRICH::SetTreeAddress()
     
   if(fLoader->TreeD()){//D    
     if(GetDebug())Info("SetTreeAddress","tree D is requested.");
-    for(int i=0;i<kNCH;i++){      
+    for(int i=0;i<kNchambers;i++){      
       branch=fLoader->TreeD()->GetBranch(Form("%s%d",GetName(),i+1)); 
       if(branch){CreateDigits(); branch->SetAddress(&((*fDigitsNew)[i]));}
     }
@@ -428,7 +426,7 @@ void AliRICH::SetTreeAddress()
     
   if(fLoader->TreeR()){//R
     if(GetDebug())Info("SetTreeAddress","tree R is requested.");
-    for(int i=0;i<kNCH;i++){         
+    for(int i=0;i<kNchambers;i++){         
       branch=fLoader->TreeR()->GetBranch(Form("%sClusters%d" ,GetName(),i+1));
       if(branch){CreateClusters(); branch->SetAddress(&((*fClusters)[i]));}
     }
@@ -440,14 +438,14 @@ void AliRICH::Print(Option_t *option)const
 {
 //Debug printout
   TObject::Print(option);
-  P()->Dump();
-  fChambers->Print(option);  
+  P()->Print();
+  fCounters.Print();
 }//void AliRICH::Print(Option_t *option)const
 //__________________________________________________________________________________________________
 void AliRICH::CreateGeometry()
 {
 //Creates detailed geometry simulation (currently GEANT volumes tree)         
-  if(GetDebug())Info("CreateGeometry","Start.");
+  if(GetDebug())Info("CreateGeometry","Start main.");
 //Opaque quartz thickness
   Float_t oquaThickness = .5;
 //CsI dimensions
@@ -460,19 +458,16 @@ void AliRICH::CreateGeometry()
   Float_t zs;
   Int_t idrotm[1099];
   Float_t par[3];
-    
-//External aluminium box 
-  par[0]=68.8;par[1]=13;par[2]=70.86;  gMC->Gsvolu("RICH", "BOX ", idtmed[1009], par, 3);
-//Air 
-  par[0]=66.3;   par[1] = 13; par[2] = 68.35;      gMC->Gsvolu("SRIC", "BOX ", idtmed[1000], 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);
+//  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);
+//  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);
-//Aluminium sheet 
-  par[0]=66.3;par[1]=0.025;par[2]=68.35;           gMC->Gsvolu("ALUM", "BOX ", idtmed[1009], 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;
@@ -502,7 +497,7 @@ void AliRICH::CreateGeometry()
 //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()->ProximityGap()/2;par[2]=pcY/2;gMC->Gsvolu("GAP ","BOX ",(*fIdtmed)[kGAP],par,3);
+  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 
@@ -561,10 +556,10 @@ void AliRICH::CreateGeometry()
   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("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");
@@ -621,11 +616,11 @@ void AliRICH::CreateGeometry()
   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()->ProximityGap()/2 - 0.0001, 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()->ProximityGap()/2 - .1, 0., 0, "ONLY");
+  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
@@ -638,100 +633,244 @@ void AliRICH::CreateGeometry()
 //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");
-
+  
 //place chambers into mother volume ALIC
-  for(int i=1;i<=kNCH;i++){
+  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");
   }
 
-  if(GetDebug())Info("CreateGeometry","Stop.");  
+         
+  if(GetDebug())Info("CreateGeometry","Stop main.");  
 }//void AliRICH::CreateGeometry()
 //__________________________________________________________________________________________________
-void AliRICH::CreateChambers()
+void AliRICH::Reconstruct()const
 {
-//create all RICH Chambers on each call. Previous chambers deleted
-  if(fChambers) delete fChambers;
-  if(GetDebug())Info("CreateChambers","Creating RICH chambers.");
-  fChambers=new TObjArray(kNCH);
-  fChambers->SetOwner();
-  for(int i=0;i<kNCH;i++)  fChambers->AddAt(new AliRICHChamber(i+1,P()),i);  
-}//void AliRICH::CreateChambers()
+  AliRICHClusterFinder finder(const_cast<AliRICH*>(this));
+  finder.Exec();
+}
 //__________________________________________________________________________________________________
-void AliRICH::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
-{
-// Generate FeedBack photons 
-// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
-  
-  TLorentzVector x4;
-  gMC->TrackPosition(x4);  
-  TVector2 x2=C(iChamber)->Glob2Loc(x4);
-  Int_t sector=P()->Sector(x2);  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 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
-  for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
-    Double_t ranf[2];
-    gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
-    cthf=ranf[0]*2-1.0;
-    if(cthf<0) continue;
-    sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
-    phif = ranf[1] * 2 * TMath::Pi();
+void AliRICH::ControlPlots()
+{ 
+//Creates a set of hists to control the results of simulation
+     
+  TH1F *pHxD=0,*pHyD=0,*pNumClusH1=0,
+                   *pQdcH1=0,       *pSizeH1=0,
+                   *pPureMipQdcH1=0,*pPureMipSizeH1=0,
+                   *pPureCerQdcH1=0,*pPureCerSizeH1=0,
+                   *pPureFeeQdcH1=0,*pPureFeeSizeH1=0,
+                   *pMipQdcH1=0,    *pPhotQdcH1=0;  
+  TH2F *pMapH2=0,*pPureMipMapH2=0,*pPureCerMapH2=0,*pPureFeeMapH2=0;
+  
+  Bool_t isDig =!GetLoader()->LoadDigits();
+  Bool_t isClus=!GetLoader()->LoadRecPoints();
+
+  if(!isDig && !isClus){Error("ControlPlots","No digits and clusters! Nothing to do.");return;}
+  
+  TStopwatch sw;TDatime time;
     
-    if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
-      enfp = 7.5e-9;
-    else if(randomNumber<=0.7)
-      enfp = 6.4e-9;
-    else
-      enfp = 7.9e-9;
+  TFile *pFile = new TFile("$(HOME)/RCP.root","RECREATE");   
+  
+  if(isDig){
+    cout<<"Digits available\n";
+    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
+  
+  if(isClus){ 
+    cout<<"Clusters available\n";
+    pNumClusH1=new TH1F("NumClusPerEvent","Number of clusters per event;number",50,0,49);
     
-
-    dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
-    gMC->Gdtom(dir, mom, 2);
-    mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
-    mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
+    pQdcH1        =new TH1F("ClusQdc",   "Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+    pSizeH1       =new TH1F("ClusSize",  "Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+    pMapH2        =new TH2F("ClusMap",   "Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY());
+  
+    pMipQdcH1     =new TH1F("QdcMip"      ,"MIP Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+    pPhotQdcH1    =new TH1F("QdcPhot"     ,"Cer+Fee Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+        
+    pPureMipQdcH1 =new TH1F("QdcPureMip"  ,"MIP only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+    pPureMipSizeH1=new TH1F("SizePureMip" ,"MIP only Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+    pPureMipMapH2 =new TH2F("MapPureMip"  ,"MIP only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY());
+  
+    pPureCerQdcH1 =new TH1F("QdcPureCer"  ,"Cerenkov only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+    pPureCerSizeH1=new TH1F("SizePureCer" ,"Cernekov only Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+    pPureCerMapH2 =new TH2F("MapPureCer"  ,"Cerenkov only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY());
     
-    // Polarisation
-    e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
-    e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
-    e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
+    pPureFeeQdcH1 =new TH1F("QdcPureFee"  ,"Feedback only Cluster Charge all chambers;q [QDC]",R()->P()->MaxQdc(),0,R()->P()->MaxQdc());
+    pPureFeeSizeH1=new TH1F("SizePureFee" ,"Feedback only Cluster size all chambers;size [number of pads in cluster]",100,0,100);
+    pPureFeeMapH2 =new TH2F("MapPureFee"  ,"Feedback only Cluster map;x [cm];y [cm]",1000,0,R()->P()->PcSizeX(),1000,0,R()->P()->PcSizeY());
+  }//isClus
+  
+  for(Int_t iEvtN=0;iEvtN < GetLoader()->GetRunLoader()->GetAliRun()->GetEventsPerRun();iEvtN++){//events loop
+    GetLoader()->GetRunLoader()->GetEvent(iEvtN);    //gets current event
     
-    vmod=0;
-    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
-    if (!vmod) for(j=0;j<3;j++) {
-      uswop=e1[j];
-      e1[j]=e3[j];
-      e3[j]=uswop;
-    }
-    vmod=0;
-    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
-    if (!vmod) for(j=0;j<3;j++) {
-      uswop=e2[j];
-      e2[j]=e3[j];
-      e3[j]=uswop;
+    if(!GetLoader()->TreeH()) GetLoader()->LoadHits();
+    for(Int_t iPrimN=0;iPrimN < GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
+       GetLoader()->TreeH()->GetEntry(iPrimN);      
     }
     
-    vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;    
-    vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
+    if(isClus) GetLoader()->TreeR()->GetEntry(0);
+    if(isDig)  GetLoader()->TreeD()->GetEntry(0);  
     
-    phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
-    for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
-    gMC->Gdtom(pol, pol, 2);
-    Int_t outputNtracksStored;    
-    gAlice->GetMCApp()->PushTrack(1,                 //do not transport
-                     gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track 
-                     kFeedback,                      //PID
-                    mom[0],mom[1],mom[2],mom[3],    //track momentum  
-                     x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin 
-                     pol[0],pol[1],pol[2],           //polarization
-                    kPFeedBackPhoton,
-                     outputNtracksStored,
-                     1.0);    
-  }//feedbacks loop
-}//GenerateFeedbacks()
+    for(Int_t iChamN=1;iChamN<=7;iChamN++){//chambers loop
+      if(isClus){
+        Int_t iNclusCham=Clusters(iChamN)->GetEntries(); if(iNclusCham) pNumClusH1->Fill(iNclusCham);//number of clusters per event
+        for(Int_t iClusN=0;iClusN<iNclusCham;iClusN++){//clusters loop
+          AliRICHcluster *pClus=(AliRICHcluster*)Clusters(iChamN)->At(iClusN);
+                                       pQdcH1        ->Fill(pClus->Q());   
+                                       pSizeH1       ->Fill(pClus->Size());  
+                                       pMapH2        ->Fill(pClus->X(),pClus->Y()); //common
+                                       
+           if(pClus->IsSingleMip())     {pPureMipQdcH1 ->Fill(pClus->Q());
+                                       pPureMipSizeH1->Fill(pClus->Size());
+                                       pPureMipMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Mips
+                                       
+           if(pClus->IsSingleCerenkov()){pPureCerQdcH1 ->Fill(pClus->Q());
+                                       pPureCerSizeH1->Fill(pClus->Size());
+                                       pPureCerMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Cerenkovs
+                                       
+           if(pClus->IsSingleFeedback()){pPureFeeQdcH1 ->Fill(pClus->Q());
+                                       pPureFeeSizeH1->Fill(pClus->Size());
+                                       pPureFeeMapH2 ->Fill(pClus->X(),pClus->Y());}//Pure Feedbacks
+           
+           if(pClus->IsMip()) {pMipQdcH1 ->Fill(pClus->Q());} //MIP+ other contributions
+           if(!pClus->IsPureMip())     pPhotQdcH1->Fill(pClus->Q());  //not MIP
+        }//clusters loop
+      }//isClus
+      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());
+          pHxD->Fill(hitV2.X()-digV2.X()); pHyD->Fill(hitV2.Y()-digV2.Y());
+        }//digits loop
+      }//isDig
+    }//chambers loop
+    Info("ControlPlots","Event %i processed.",iEvtN);
+  }//events loop 
+  
+  if(isDig)  GetLoader()->UnloadDigits();
+  if(isClus) GetLoader()->UnloadRecPoints();
+  
+  pFile->Write(); delete pFile;
+  sw.Print();time.Print();
+}//ControlPlots()
+//__________________________________________________________________________________________________
+AliRICHhit* AliRICH::Hit(Int_t tid)
+{
+//defines which hit provided by given tid for the currently loaded event
+  R()->GetLoader()->LoadHits();
+  for(Int_t iPrimN=0;iPrimN<R()->GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop      
+    R()->GetLoader()->TreeH()->GetEntry(iPrimN);
+    for(Int_t iHitN=0;iHitN<R()->Hits()->GetEntries();iHitN++){
+      AliRICHhit *pHit=(AliRICHhit*)R()->Hits()->At(iHitN);
+      if(tid==pHit->Track()) {R()->GetLoader()->UnloadHits();return pHit;}
+    }//hits
+  }//prims loop
+  R()->GetLoader()->UnloadHits();
+  return 0;
+}
+//__________________________________________________________________________________________________
+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);
+  R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
+  if(R()->GetLoader()->LoadHits()) return;
+  
+  Int_t iTotalHits=0;
+  for(Int_t iPrimN=0;iPrimN<R()->GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
+    R()->GetLoader()->TreeH()->GetEntry(iPrimN);      
+    R()->Hits()->Print();
+    iTotalHits+=R()->Hits()->GetEntries();
+  }
+  R()->GetLoader()->UnloadHits();
+  Info("PrintHits","totally %i hits",iTotalHits);
+}
+//__________________________________________________________________________________________________
+void AliRICH::PrintSDigits(Int_t iEvtN)
+{
+//prints a list of RICH sdigits  for a given event
+  Info("PrintSDigits","List of RICH sdigits for event %i",iEvtN);
+  R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
+  if(R()->GetLoader()->LoadSDigits()) return;
+  
+  R()->GetLoader()->TreeS()->GetEntry(0);
+  R()->SDigits()->Print();
+  R()->GetLoader()->UnloadSDigits();
+  Info("PrintSDigits","totally %i sdigits",R()->SDigits()->GetEntries());
+}
+//__________________________________________________________________________________________________
+void AliRICH::PrintDigits(Int_t iEvtN)
+{
+//prints a list of RICH digits  for a given event
+  Info("PrintDigits","List of RICH digits for event %i",iEvtN);
+  R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
+  if(R()->GetLoader()->LoadDigits()) return;
+  
+  Int_t iTotalDigits=0;
+  R()->GetLoader()->TreeD()->GetEntry(0);
+  for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){
+    R()->Digits(iChamber)->Print();
+    iTotalDigits+=R()->Digits(iChamber)->GetEntries();
+  }
+  R()->GetLoader()->UnloadDigits();
+  Info("PrintDigits","totally %i Digits",iTotalDigits);
+}
+//__________________________________________________________________________________________________
+void AliRICH::PrintClusters(Int_t iEvtN)
+{
+//prints a list of RICH clusters  for a given event
+  Info("PrintClusters","List of RICH clusters for event %i",iEvtN);
+  R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
+  if(R()->GetLoader()->LoadRecPoints()) return;
+  
+  Int_t iTotalClusters=0;
+  R()->GetLoader()->TreeR()->GetEntry(0);
+  for(Int_t iChamber=1;iChamber<=kNchambers;iChamber++){
+    R()->Clusters(iChamber)->Print();
+    iTotalClusters+=R()->Clusters(iChamber)->GetEntries();
+  }
+  R()->GetLoader()->UnloadRecPoints();
+  Info("PrintClusters","totally %i clusters",iTotalClusters);
+}
+//__________________________________________________________________________________________________
+void AliRICH::FillESD(AliESD *pESD)const
+{
+//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) {pTrack->ResetRICH(); 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
+}
index 88a51e8..fa06ed9 100644 (file)
@@ -18,7 +18,7 @@
 class AliRICHhit : public AliHit
 {
 public:
-           AliRICHhit():AliHit()          {fChamber=kBad;  fEloss=kBad;  fInX3.SetXYZ(0,0,0);fOutX3.SetXYZ(0,0,0);}
+           AliRICHhit():AliHit(),fChamber(kBad),fEloss(kBad) {fInX3.SetXYZ(0,0,0);fOutX3.SetXYZ(0,0,0);}
            AliRICHhit(Int_t c,Int_t tid,TVector3 in,TVector3 out,Double_t e):AliHit(0,tid)
                                           {fChamber=c;fInX3=in; fOutX3=out;fEloss=e; fX=out.X();fY=out.Y();fZ=out.Z();}
   virtual ~AliRICHhit()                   {;}
@@ -42,26 +42,26 @@ protected:
 class AliRICHdigit :public AliDigit
 {
 public:
-  AliRICHdigit() {fChFbMip=fChamber=fPadX=fPadY=fTracks[0]=fTracks[1]=fTracks[2]=kBad;fQdc=kBad;}
-  AliRICHdigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t cpid,Int_t tid0,Int_t tid1,Int_t tid2)  
-                 {fPadX=x;fPadY=y;fQdc=q;fChamber=10*c+AliRICHParam::Sector(x,y);fChFbMip=cpid;fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;}
+  AliRICHdigit():AliDigit(),fCFM(0),fChamber(0),fPadX(0),fPadY(0),fQdc(kBad){fTracks[0]=fTracks[1]=fTracks[2]=kBad;}
+  AliRICHdigit(Int_t c,TVector pad,Double_t q,Int_t cfm,Int_t tid0,Int_t tid1,Int_t tid2)  
+       {fPadX=(Int_t)pad[0];fPadY=(Int_t)pad[1];fQdc=q;fChamber=10*c+AliRICHParam::Pad2Sec(pad);fCFM=cfm;fTracks[0]=tid0;fTracks[1]=tid1;fTracks[2]=tid2;}
   virtual ~AliRICHdigit() {;}  
   Int_t    Compare(const TObject *pObj) const
                  {if(Id()==((AliRICHdigit*)pObj)->Id())return 0;else if(Id()>((AliRICHdigit*)pObj)->Id())return 1;else return -1;}  //virtual      
-  Bool_t   IsSortable()                 const{return kTRUE;}                              //virtual
-  void     Print(Option_t *option="")   const;                                            //virtual
-  Int_t    ChFbMi()                     const{return fChFbMip;}                           //particle mixture for this digit
+  virtual Bool_t   IsSortable()                 const{return kTRUE;}                              //sort interface
+  virtual void     Print(Option_t *option="")   const;                                            //virtual
+  Int_t    ChFbMi()                     const{return fCFM;}                               //particle mixture for this digit
   Int_t    C()                          const{return fChamber/10;}                        //chamber number
   Int_t    S()                          const{return fChamber-(fChamber/10)*10;}          //sector number
   Int_t    X()                          const{return fPadX;}                              //x position of the pad
   Int_t    Y()                          const{return fPadY;}                              //y postion of the pad
+  TVector  Pad()                        const{Float_t v[2]={fPadX,fPadY}; return TVector(2,v);}
   Int_t    Id()                         const{return fChamber*10000000+fPadX*1000+fPadY;} //absolute id of this pad
   Double_t Q()                          const{return fQdc;}                               //charge in terms of ADC channels
-  Int_t    Tid(Int_t i)                 const{return fTracks[i];}                         //track reference produced this digit
   void     AddTidOffset(Int_t offset) 
     {for (Int_t i=0; i<3; i++) if (fTracks[i]>0) fTracks[i]+=offset;};
 protected:
-  Int_t    fChFbMip;  //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
+  Int_t    fCFM;  //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
   Int_t    fChamber;  //10*chamber number+ sector number 
   Int_t    fPadX;     //pad number along X
   Int_t    fPadY;     //pad number along Y
@@ -73,42 +73,48 @@ protected:
 class AliRICHcluster :public TObject
 {
 public:
-  enum ClusterStatus {kEdge,kShape,kSize,kRaw,kResolved};
-                    AliRICHcluster()                                 {fSize=fQdc=fStatus=fChamber=fDimXY=0;fX=fY=kBad;fDigits=0;}
-  virtual          ~AliRICHcluster()                                 {delete fDigits;}  
-  AliRICHcluster&   operator=(const AliRICHcluster&)                 {return *this;}
-         Int_t      Nlocals()                                   const{return fSize - 10000*(fSize/10000);} //
-         Int_t      Size()                                      const{return fSize/10000;}                 //
-         Int_t      Fsize()                                     const{return fSize;}                       //
-         Int_t      DimXY()                                     const{return fDimXY;}                      //
-         Int_t      C()                                         const{return fChamber/10;}                 //
-         Int_t      S()                                         const{return fChamber-(fChamber/10)*10;}   //
-         Int_t      Fchamber()                                  const{return fChamber;}                    //
-         Int_t      Q()                                         const{return fQdc;}                        // 
-         Double_t   X()                                         const{return fX;}                          //
-         Double_t   Y()                                         const{return fY;}                          // 
-         Int_t      Status()                                    const{return fStatus;}                     //
-         void       SetStatus(Int_t status)                          {fStatus=status;}                     //
-         Int_t      Nmips()                                     const{return fCombiPid-1000000*Ncerenkovs()-1000*Nfeedbacks();}  //
-         Int_t      Ncerenkovs()                                const{return fCombiPid/1000000;}           //
-         Int_t      Nfeedbacks()                                const{return (fCombiPid-1000000*Ncerenkovs())/1000;}      //
-         Bool_t     IsPureMip()                                 const{return fCombiPid<1000;}   
-         Bool_t     IsPureCerenkov()                            const{return Nmips()==0&&Nfeedbacks()==0;} //
-         Bool_t     IsPureFeedback()                            const{return Nmips()==0&&Ncerenkovs()==0;} //
-         Int_t      CombiPid()                                  const{return fCombiPid;} //
-         void       SetCombiPid(Int_t ckov,Int_t feeds,Int_t mips)   {fCombiPid=1000000*ckov+1000*feeds+mips;}            //
-         void       Fill(AliRICHcluster *pRaw,Double_t x,Double_t y, Double_t q, Int_t combipid)
-             {fCombiPid=combipid;fChamber=pRaw->Fchamber();fSize=pRaw->Fsize();
-              fQdc=(Int_t)(q*pRaw->Q());fX=x;fY=y;fStatus=kResolved;} //
-         TObjArray* Digits()                                    const{return fDigits;}                     //  
-         void       Print(Option_t *option="")const;                                                       //virtual
-  inline void       AddDigit(AliRICHdigit *pDig);                                                          //
-  inline void       CoG(Int_t nLocals);                                                                                 // 
-         void       Reset() {fSize=fQdc=fStatus=fChamber=fDimXY=kBad;fX=fY=kBad;delete fDigits;fDigits=0;} //
+  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();}  
+         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  
+  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
+         Int_t      Fsize()                     const{return fSize;}                                    //
+         Int_t      Shape()                     const{return fShape;}                                   //cluster shape rectangulare
+         Int_t      C()                         const{return fChamber/10;}                              //chamber number
+         Int_t      S()                         const{return fChamber-(fChamber/10)*10;}                //sector number
+         Int_t      Fchamber()                  const{return fChamber;}                                 //
+         Int_t      Q()                         const{return fQdc;}                                     //cluster charge in QDC channels 
+         Double_t   X()                         const{return fX;}                                       //cluster x position in LRS
+         Double_t   Y()                         const{return fY;}                                       //cluster y position in LRS 
+         Int_t      Status()                    const{return fStatus;}                                  //
+         void       SetStatus(Int_t status)         {fStatus=status;}                                     //
+         Int_t      Nmips()                     const{return fCFM-1000000*Ncerenkovs()-1000*Nfeedbacks();} //
+         Int_t      Ncerenkovs()                const{return fCFM/1000000;}                                //
+         Int_t      Nfeedbacks()                const{return (fCFM-1000000*Ncerenkovs())/1000;}            //
+         Bool_t     IsPureMip()                 const{return fCFM<1000;}                                   //
+         Bool_t     IsPureCerenkov()            const{return Nmips()==0&&Nfeedbacks()==0;}                 //
+         Bool_t     IsPureFeedback()            const{return Nmips()==0&&Ncerenkovs()==0;}                 //
+         Bool_t     IsSingleMip()               const{return Nmips()==1&&Ncerenkovs()==0&&Nfeedbacks()==0;}  //
+         Bool_t     IsSingleCerenkov()          const{return Nmips()==0&&Ncerenkovs()==1&&Nfeedbacks()==0;}  //
+         Bool_t     IsSingleFeedback()          const{return Nmips()==0&&Ncerenkovs()==0&&Nfeedbacks()==1;}  //
+         Bool_t     IsMip()                     const{return Nmips()!=0;}                                  //
+         Bool_t     IsCerenkov()                const{return Ncerenkovs()!=0;}                             //
+         Bool_t     IsFeedback()                const{return Nfeedbacks()!=0;}                             //
+         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;} //
 protected:
-  Int_t         fCombiPid;    //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
+  Int_t         fCFM;         //1000000*Ncerenkovs+1000*Nfeedbacks+Nmips  
   Int_t         fSize;        //10000*(how many digits belong to this cluster) + nLocalMaxima     
-  Int_t         fDimXY;       //100*xdim+ydim box containing the cluster
+  Int_t         fShape;       //100*xdim+ydim box containing the cluster
   Int_t         fQdc;         //QDC value
   Int_t         fChamber;     //10*module number+sector number 
   Double_t      fX;           //local x postion 
@@ -119,53 +125,35 @@ protected:
 };//class AliRICHcluster
 //__________________________________________________________________________________________________
 void AliRICHcluster::AddDigit(AliRICHdigit *pDig)
-{//    
-  if(!fDigits) {fQdc=fSize=fCombiPid=0;fDigits = new TObjArray;}
+{
+// Adds a given digit to the list of digits belonging to this cluster    
+  if(!fDigits) {fQdc=fSize=fCFM=0;fDigits = new TObjArray;}
   fQdc+=(Int_t)pDig->Q(); fDigits->Add(pDig);
   fChamber=10*pDig->C()+pDig->S();
   fSize+=10000;
+  fStatus=kRaw;
 }
 //__________________________________________________________________________________________________
 void AliRICHcluster::CoG(Int_t nLocals)
-{//
-  Int_t xmin=999,ymin=999,xmax=0,ymax=0;   
+{
+// Calculates naive cluster position as a center of gravity of its digits.
+  Float_t xmin=999,ymin=999,xmax=0,ymax=0;   
   fX=fY=0;
   for(Int_t iDig=0;iDig<Size();iDig++) {
     AliRICHdigit *pDig=(AliRICHdigit*)fDigits->At(iDig);
-    Int_t padX = pDig->X();Int_t padY = pDig->Y();Double_t q=pDig->Q();
-    TVector2 x2=AliRICHParam::Pad2Loc(padX,padY);
+    TVector pad=pDig->Pad(); Double_t q=pDig->Q();
+    TVector2 x2=AliRICHParam::Pad2Loc(pad);
     fX += x2.X()*q;fY +=x2.Y()*q;
-    if(padX<xmin)xmin=padX;if(padX>xmax)xmax=padX;if(padY<ymin)ymin=padY;if(padY>ymax)ymax=padY;
+    if(pad[0]<xmin)xmin=pad[0];if(pad[0]>xmax)xmax=pad[0];if(pad[1]<ymin)ymin=pad[1];if(pad[1]>ymax)ymax=pad[1];
    }
    fX/=fQdc;fY/=fQdc;//Center of Gravity
-   fDimXY = 100*(xmax-xmin+1)+ymax-ymin+1;//find box containing cluster
+   fShape=Int_t(100*(xmax-xmin+1)+ymax-ymin+1);//find box containing cluster
    fSize+=nLocals;
    fStatus=kRaw;
 }//CoG()
 
-
-class AliRICHreco: public TObject
-{
-public:
-            AliRICHreco() {fTid=fNphotons=kBad; fThetaCherenkov=kBad;}
-            AliRICHreco(Int_t tid,Double_t thetaCherenkov,Int_t nPhotons) {fTid=tid;fThetaCherenkov=thetaCherenkov;fNphotons=nPhotons;}
-
-  virtual  ~AliRICHreco() {;}
-
-  void    Print(Option_t *option="")const;      //virtual print
-  
-protected:
-   Int_t fTid;                 // track Id reference
-   Int_t fNphotons;            // number of photons contributed to the recontruction
-   Double_t fThetaCherenkov;   // reconstructed Theta Cerenkov for a given charged track
-
-   ClassDef(AliRICHreco,1)  //RICH reco class
-
-};//class AliRICHreco
-
 //__________________AliRICH_________________________________________________________________________
-class AliRICHParam;
-class AliRICHChamber;
+class AliESD;
 
 class AliRICH : public AliDetector 
 {
@@ -177,69 +165,66 @@ public:
           
   AliRICH&  operator=(const AliRICH&)                 {return *this;}
 //framework part  
-  virtual Int_t   IsVersion()                                            const =0;                                //virtual         
-  virtual void    StepManager()                                                =0;                                //virtual
-          void    Hits2SDigits();                                                                                 //virtual
-  AliDigitizer*   CreateDigitizer(AliRunDigitizer* man) const {return new AliRICHDigitizer(man);}                 //virtual
-          void    Print(Option_t *option)const;                                                                   //virtual
-          void    SetTreeAddress();                                                                               //virtual
-          void    MakeBranch(Option_t *opt=" ");                                                                  //virtual
-          void    CreateMaterials();                                                                              //virtual
-  virtual void    BuildGeometry();                                                                                //virtual
-  virtual void    CreateGeometry();                                                                               //virtual
+  virtual Int_t         IsVersion()                           const =0;                                  //interface from         
+  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          SetTreeAddress();                                                                //interface from AliLoader
+  virtual void          MakeBranch(Option_t *opt=" ");                                                   //interface from AliLoader
+  virtual void          CreateMaterials();                                                               //interface from AliMC
+  virtual void          CreateGeometry();                                                                //interface from AliMC
+  virtual void          BuildGeometry();                                                                 //interface 
 //private part  
           Float_t AbsoCH4(Float_t x)const;                               //calculates absorption length for methane
           Float_t Fresnel(Float_t ene,Float_t pdoti, Bool_t pola)const;  //deals with Fresnel absorption
-          void    GenerateFeedbacks(Int_t iChamber,Float_t eloss=0);     //generates feedback photons; eloss=0 for photon
-          void    CreateChambers();                                      //creates set of chambers
   inline  void    CreateHits();                                          //create hits container as a simple list
   inline  void    CreateSDigits();                                       //create sdigits container as a simple list
   inline  void    CreateDigits();                                        //create digits container as 7 lists, one per chamber
   inline  void    CreateClusters();                                      //create clusters container  as 7 lists, one per chamber
-  inline  void    CreateRecos();                                         //create recos container
 //        void    ResetHits()                {AliDetector::ResetHits();}  //virtual  
           void    ResetSDigits()             {fNsdigits=0;  if(fSdigits)  fSdigits ->Clear();}                                 
-          void    ResetDigits()              {if(fDigitsNew)for(int i=0;i<kNCH;i++){fDigitsNew->At(i)->Clear();fNdigitsNew[i]=0;}} //virtual
-          void    ResetClusters()            {if(fClusters) for(int i=0;i<kNCH;i++){fClusters ->At(i)->Clear();fNclusters[i]=0;}}
-          void    ResetRecos()               {if(fRecos) fRecos->Clear();fNrecos=0;}
+          void    ResetDigits()              {if(fDigitsNew)for(int i=0;i<kNchambers;i++){fDigitsNew->At(i)->Clear();fNdigitsNew[i]=0;}} //virtual
+          void    ResetClusters()            {if(fClusters) for(int i=0;i<kNchambers;i++){fClusters ->At(i)->Clear();fNclusters[i]=0;}}
   TClonesArray*   SDigits()             const{return fSdigits;}
   TClonesArray*   Digits(Int_t iC)      const{if(fDigitsNew) return (TClonesArray *)fDigitsNew->At(iC-1);else return 0;}
   TClonesArray*   Clusters(Int_t iC)    const{if(fClusters)  return (TClonesArray *)fClusters->At(iC-1);else return 0;}
-  TClonesArray*   Recos()               const{return fRecos;}
-  AliRICHChamber* C(Int_t iC)           const{return (AliRICHChamber*)fChambers->At(iC-1);}
-  AliRICHParam*   P()                   const{return fpParam;}
-          void    PrintDigits()         const{for(Int_t i=0;i<kNCH;i++) fDigitsNew->At(i)->Print();}
-          void    PrintClusters()       const{for(Int_t i=0;i<kNCH;i++) fClusters->At(i)->Print();}
+  AliRICHChamber* C(Int_t iC)           const{return fpParam->C(iC);}                       //provides pointer to a given chamber
+  AliRICHParam*   P()                   const{return fpParam;}                              //provides pointer to a RICH params
+  AliRICH*        R()                        {return this;}                                 //provides pointer to RICH main object
+  TVector         Counters()            const{return fCounters;}                            //provides a set of counters
+  void            ControlPlots();                                                           //utility
+  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 AddHit(Int_t chamber,Int_t tid,TVector3 iX3,TVector3 oX3,Double_t eloss=0)
-       {TClonesArray &tmp=*fHits;new(tmp[fNhits++])AliRICHhit(chamber,tid,iX3,oX3,eloss);} 
-  inline void AddSDigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t pid,Int_t tid); 
-  void AddDigit(int c,int x,int y,int q,int cfm,int *tid)
-       {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(c-1));new(tmp[fNdigitsNew[c-1]++])AliRICHdigit(c,x,y,q,cfm,tid[0],tid[1],tid[2]);}  
+  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); 
+  void AddDigit(int c,TVector pad,int q,int cfm,int *tid)//Add simulated digit
+       {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(c-1));new(tmp[fNdigitsNew[c-1]++])AliRICHdigit(c,pad,q,cfm,tid[0],tid[1],tid[2]);}  
+  void AddDigit(Int_t c,TVector pad,Int_t q)//for real data digits
+       {TClonesArray &tmp=*((TClonesArray*)fDigitsNew->At(0));new(tmp[fNdigitsNew[0]++])AliRICHdigit(c,pad,q,0,-1,-1,-1);}  
   void AddCluster(AliRICHcluster &cl)                     
-       {Int_t c=cl.C()-1;/*cout<<c<<endl*/;TClonesArray &tmp=*((TClonesArray*)fClusters->At(c));new(tmp[fNclusters[c]++])AliRICHcluster(cl);}
-  void AddReco(Int_t tid,Double_t thetaCherenkov,Int_t nPhotons) 
-       {TClonesArray &tmp=*(TClonesArray*)fRecos;new(tmp[fNrecos++])AliRICHreco(tid,thetaCherenkov,nPhotons);}  
-          
+       {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                  {kCSI=6,kGAP=9};
+  enum                  {kAir=1,kCSI=6,kGAP=9,kAl=10,kCH4=5,kSteel=15,kPerpex=16,kSr90=17};
   AliRICHParam         *fpParam;             //main RICH parametrization     
-  TObjArray            *fChambers;           //list of RICH chambers
-  
                                              //fHits and fDigits belong to AliDetector
-  TClonesArray         *fSdigits;            //! List of sdigits  
-  Int_t                 fNsdigits;           //! Current number of sdigits
+  TClonesArray         *fSdigits;            //! list of sdigits  
+  Int_t                 fNsdigits;           //! current number of sdigits
   
-  TObjArray            *fDigitsNew;          //! Each chamber holds it's one lists of digits
-  Int_t                 fNdigitsNew[kNCH];   //! Array of current numbers of digits
+  TObjArray            *fDigitsNew;          //! each chamber holds it's one lists of digits
+  Int_t                 fNdigitsNew[kNchambers];   //! array of current numbers of digits
   
-  TObjArray            *fClusters;           //! Each chamber holds it's one lists of clusters 
-  Int_t                 fNclusters[kNCH];    //! Array of current numbers of raw clusters
+  TObjArray            *fClusters;           //! each chamber holds it's one lists of clusters 
+  Int_t                 fNclusters[kNchambers];    //! array of current numbers of raw clusters
   
-  TClonesArray         *fRecos;              //! pointer to the list of recos
-  Int_t                 fNrecos;             //! number of recos
-
-  ClassDef(AliRICH,5)                        //Main RICH class 
+  TVector               fCounters;           //Photon history conters, explanation in StepManager() 
+  ClassDef(AliRICH,7)                        //Main RICH class 
 };//class AliRICH  
 
 //__________________________________________________________________________________________________
@@ -261,34 +246,27 @@ void AliRICH::CreateDigits()
 {
   if(fDigitsNew) return;
   if(GetDebug())Info("CreateDigits","creating digits containers.");
-  fDigitsNew = new TObjArray(kNCH);  
-  for(Int_t i=0;i<kNCH;i++) {fDigitsNew->AddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;}
+  fDigitsNew = new TObjArray(kNchambers);  
+  for(Int_t i=0;i<kNchambers;i++) {fDigitsNew->AddAt(new TClonesArray("AliRICHdigit",10000), i); fNdigitsNew[i]=0;}
 }
 //__________________________________________________________________________________________________
 void AliRICH::CreateClusters()
 {
   if(fClusters) return;
   if(GetDebug())Info("CreateClusters","creating clusters containers.");
-  fClusters = new TObjArray(kNCH);  
-  for(Int_t i=0;i<kNCH;i++) {fClusters->AddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;}
+  fClusters = new TObjArray(kNchambers);  
+  for(Int_t i=0;i<kNchambers;i++) {fClusters->AddAt(new TClonesArray("AliRICHcluster",10000), i); fNclusters[i]=0;}
 }
 //__________________________________________________________________________________________________
-void AliRICH::CreateRecos()
-{
-  if(fRecos) return;
-  if(GetDebug())Info("CreateRecos","creating recos containers.");
-  fRecos = new TClonesArray("AliRICHreco",1000);fNrecos=0;  
-}
-//__________________________________________________________________________________________________
-void AliRICH::AddSDigit(Int_t c,Int_t x,Int_t y,Double_t q,Int_t pid,Int_t tid) 
-{   
+void AliRICH::AddSDigit(Int_t c,TVector pad,Double_t q,Int_t pid,Int_t tid) 
+{ 
+  Int_t cfm;  
   switch(pid){
-    case 50000050: pid=1000000;break;//cerenkov
-    case 50000051: pid=1000;   break;//feedback
-    default:       pid=1;      break;//mip
+    case 50000050: cfm=1000000;break;//cerenkov
+    case 50000051: cfm=1000;   break;//feedback
+    default:       cfm=1;      break;//mip
   }   
-  TClonesArray &tmp=*fSdigits;
-  new(tmp[fNsdigits++])AliRICHdigit(c,x,y,q,pid,tid,kBad,kBad);
-}//AddSDigit()   
+  TClonesArray &tmp=*fSdigits; new(tmp[fNsdigits++])AliRICHdigit(c,pad,q,cfm,tid,kBad,kBad);
+}
 //__________________________________________________________________________________________________
 #endif//#ifndef AliRICH_h
index 5fd0198..7f7f86f 100644 (file)
 //  **************************************************************************
 
 #include "AliRICHChamber.h"
-#include "AliRICHParam.h"
 #include <TRotMatrix.h>
 
 ClassImp(AliRICHChamber)       
-//______________________________________________________________________________    
-AliRICHChamber::AliRICHChamber() 
-{
-//default ctor
-  fpParam=0;    
-  fpRotMatrix=0;
-}
 //______________________________________________________________________________
-AliRICHChamber::AliRICHChamber(Int_t iModuleN,AliRICHParam *pParam)
+AliRICHChamber::AliRICHChamber(Int_t iChamber):TNamed()
 {
 //main ctor. Defines all geometry parameters for the given module.
   SetToZenith();//put to up position   
-  switch(iModuleN){
+  switch(iChamber){
     case 1:
-      RotateX(-pParam->AngleYZ());   
-      RotateZ(-pParam->AngleXY());      
+      RotateX(-AliRICHParam::AngleYZ());   
+      RotateZ(-AliRICHParam::AngleXY());      
       fName="RICHc1";fTitle="RICH chamber 1";
       break;      
     case 2:
-      RotateZ(-pParam->AngleXY());      
+      RotateZ(-AliRICHParam::AngleXY());      
       fName="RICHc2";fTitle="RICH chamber 2";
       break;      
     case 3:
-      RotateX(-pParam->AngleYZ());
+      RotateX(-AliRICHParam::AngleYZ());
       fName="RICHc3";fTitle="RICH chamber 3";
       break;      
     case 4:          
       fName="RICHc4";fTitle="RICH chamber 4";  //no turns
       break;      
     case 5:
-      RotateX(pParam->AngleYZ());
+      RotateX(AliRICHParam::AngleYZ());
       fName="RICHc5";fTitle="RICH chamber 5";
       break;      
     case 6:
-      RotateZ(pParam->AngleXY());      
+      RotateZ(AliRICHParam::AngleXY());      
       fName="RICHc6";fTitle="RICH chamber 6";
       break;      
     case 7:
-      RotateX(pParam->AngleYZ());            
-      RotateZ(pParam->AngleXY());      
+      RotateX(AliRICHParam::AngleYZ());            
+      RotateZ(AliRICHParam::AngleXY());      
       fName="RICHc7";fTitle="RICH chamber 7";
       break;      
     default:
-      Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iModuleN);
+      Fatal("named ctor","Wrong chamber number %i, check CreateChamber ctor",iChamber);
   }//switch(iModuleN)
-  RotateZ(pParam->AngleRot());//apply common rotation  
+  RotateZ(AliRICHParam::AngleRot());//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());
-  fpParam=pParam;
 }//main ctor
 //__________________________________________________________________________________________________
-void AliRICHChamber::Print(Option_t *) const
+void AliRICHChamber::Print(Option_t *opt) const
 {
-//debug printout method 
-  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());
+// 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);
 }//Print()
 //__________________________________________________________________________________________________
index c6ff333..08ebdb9 100644 (file)
@@ -9,19 +9,14 @@
 #include <TRotation.h>
 #include <TLorentzVector.h>
 #include "AliRICHParam.h"
-#include "AliSegmentation.h"
-class AliRICHGeometry;
-class AliRICHResponse;
 class TRotMatrix;
 
-typedef enum {kMip, kPhoton} ResponseType;
-class AliRICHParam;
 
 class AliRICHChamber : public TNamed
 {
 public:
-           AliRICHChamber();
-           AliRICHChamber(Int_t iModuleN,AliRICHParam *pParam);
+           AliRICHChamber():TNamed(),fpRotMatrix(0)                      {;}
+           AliRICHChamber(Int_t iChamberN);
            AliRICHChamber(const AliRICHChamber &chamber):TNamed(chamber) {;}
   virtual ~AliRICHChamber()                                              {;}
            AliRICHChamber& operator=(const AliRICHChamber&)              {return *this;}
@@ -29,65 +24,45 @@ public:
   TRotMatrix* RotMatrix()          const{return fpRotMatrix;}
   TString RotMatrixName()          const{return "rot"+fName;}
   TRotation   Rot()                const{return fRot;}
-  Double_t    Rho()                const{return fCenterV3.Mag();} 
-  Double_t    ThetaD()             const{return fCenterV3.Theta()*TMath::RadToDeg();}    
-  Double_t    PhiD()               const{return fCenterV3.Phi()*TMath::RadToDeg();}    
-  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);}
-  void        RotateY(Double_t a)       {fRot.RotateY(a);fCenterV3.RotateY(a);fPcX3.RotateY(a);}
-  void        RotateZ(Double_t a)       {fRot.RotateZ(a);fCenterV3.RotateZ(a);fPcX3.RotateZ(a);}
+  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    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();}
-  TVector3    L2G(TVector3 x3)                       const{x3.Transform(fRot);x3+=fCenterV3;return x3;}
-  TVector3    G2L(TVector3 x3)                       const{x3-=fCenterV3;x3.Transform(fRot.Inverse()); return x3;}
-  inline TVector2  Glob2Loc(TVector3 x3, Bool_t isVector=kFALSE) const;
-  TVector2    Glob2Loc(TLorentzVector x4,Bool_t isVector=kFALSE) const{return Glob2Loc(x4.Vect(),isVector);}
-  TVector3    L2G(Double_t x,Double_t y,Double_t z)  const{return L2G(TVector3(x,y,z));}
-  TVector3    G2L(TLorentzVector x4)                 const{return G2L(x4.Vect());}
-  Float_t     G2Ly(TLorentzVector x4)                const{TVector3 x3=G2L(x4.Vect()); return x3.Z();}
-  TVector3    G2L(Double_t x,Double_t y,Double_t z)  const{return G2L(TVector3(x,y,z));}
-  Float_t     G2Lx(Double_t x,Double_t y,Double_t z) const{TVector3 x3=G2L(x,y,z); return x3.X();}
-  Float_t     G2Ly(Double_t x,Double_t y,Double_t z) const{TVector3 x3=G2L(x,y,z); return x3.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      
    
-  void LocaltoGlobal(Float_t pos[3],Float_t localpos[3]) {
-    //Transformation from local to global coordinates, chamber-dependant
-    TVector3 buf = L2G(localpos[0],localpos[1],localpos[2]);
-    pos[0]=buf.X();pos[1]=buf.Y();pos[2]=buf.Z();
-  }
-  void GlobaltoLocal(Float_t pos[3],Float_t localpos[3]) {
-    //Transformation from Global to local coordinates, chamber-dependant 
-    TVector3 buf = G2L(pos[0],pos[1],pos[2]);
-    localpos[0]=buf.X();localpos[1]=buf.Y();localpos[2]=buf.Z();
-  }
+
   inline void SetToZenith();
   TRotMatrix *GetRotMatrix()       const{return fpRotMatrix;}  
 protected:
-  TVector3      fCenterV3;        //chamber center position in MRS (cm)
+  TVector3      fCenterV3;        //chamber 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 
-  AliRICHParam *fpParam;          //main RICH parameters description  
-  ClassDef(AliRICHChamber,5)      //single RICH chamber description
+  ClassDef(AliRICHChamber,6)      //single RICH chamber description
 };//class AliRICHChamber
 //__________________________________________________________________________________________________
 void AliRICHChamber::SetToZenith()
 {
-  fCenterV3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2,0); 
-  fPcX3.SetXYZ(0,AliRICHParam::Offset()-AliRICHParam::GapThickness()/2+5.276+0.25,0);   
+//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);   
 }
 //__________________________________________________________________________________________________
-TVector2 AliRICHChamber::Glob2Loc(TVector3 x3,Bool_t isVector)const
-{
-  if(!isVector) x3-=fPcX3;
-  x3.Transform(fRot.Inverse()); 
-  return TVector2(x3.X(),x3.Z());//attention Y and Z are misplaced!
-}
-//__________________________________________________________________________________________________  
 #endif //AliRICHChamber_h
index 5b40315..c12ff0e 100644 (file)
@@ -29,219 +29,209 @@ ClassImp(AliRICHClusterFinder)
 //__________________________________________________________________________________________________
 AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)   
 {//main ctor
-//  Info("main ctor","Start.");
+  fRICH = pRICH;//this goes first as GetDebug depends on it.
+  if(GetDebug())Info("main ctor","Start.");
   
-  fRICH = pRICH;
-  
-  fHitMap = 0;
+  fDigitMap = 0;
   fRawCluster.Reset();
   fResolvedCluster.Reset();
-  
+  if(GetDebug())Info("main ctor","Stop.");
 }//main ctor
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::FindLocalMaxima()
-{
-//  Find number of local maxima in the cluster
-//  Info("FindLocalMaxima","Start.");
-  
-  fNlocals = 0;
-  for(Int_t iDig1=0;iDig1<fRawCluster.Size();iDig1++) {
-    Int_t iNotMax = 0;
-    AliRICHdigit *pDig1 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig1);
-    Int_t padX1 = pDig1->X();
-    Int_t padY1 = pDig1->Y();
-    Int_t padQ1 = (Int_t)(pDig1->Q()+0.1);
-    Int_t padC1 = pDig1->ChFbMi();
-    for(Int_t iDig2=0;iDig2<fRawCluster.Size();iDig2++) {
-      AliRICHdigit *pDig2 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig2);
-      Int_t padX2 = pDig2->X();
-      Int_t padY2 = pDig2->Y();
-      Int_t padQ2 = (Int_t)(pDig2->Q()+0.1);
-      if(iDig1==iDig2) continue;
-      Int_t diffx = TMath::Sign(padX1-padX2,1);
-      Int_t diffy = TMath::Sign(padY1-padY2,1);
-      if((diffx+diffy)<=1) {
-         if(padQ2>=padQ1) iNotMax++;
-      }
-    }
-    if(iNotMax==0) {
-      if (fNlocals >= fLocalX.GetSize()) {
-       fLocalX.Set(fNlocals+100);
-       fLocalY.Set(fNlocals+100);
-       fLocalQ.Set(fNlocals+100);
-       fLocalC.Set(fNlocals+100);
-      }
-      TVector2 x2=AliRICHParam::Pad2Loc(padX1,padY1);
-      fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y();
-      fLocalQ[fNlocals] = (Double_t)padQ1;
-      fLocalC[fNlocals] = padC1;
-      fNlocals++;
-    }
-  }
-}//FindLocalMaxima()
-//__________________________________________________________________________________________________
 void AliRICHClusterFinder::Exec()
 {
+//Main method of cluster finder. Loops on  events and chambers, everything else is done in FindClusters()  
   if(GetDebug()) Info("Exec","Start.");
-  
-  
-  Rich()->GetLoader()->LoadDigits(); 
-  
-  Rich()->GetLoader()->GetRunLoader()->LoadHeader();
-  Rich()->GetLoader()->GetRunLoader()->LoadKinematics();
+    
+  R()->GetLoader()                ->LoadDigits();   
+  R()->GetLoader()->GetRunLoader()->LoadHeader();
+  R()->GetLoader()->GetRunLoader()->LoadKinematics();
 
   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
-    if(GetDebug()) Info("Exec","Event %i processed.",iEventN+1);
-//    gAlice->GetRunLoader()->GetEvent(iEventN);
-    Rich()->GetLoader()->GetRunLoader()->GetEvent(iEventN);
+    if(GetDebug()) Info("Exec","Processing event %i...",iEventN);
+    R()->GetLoader()->GetRunLoader()->GetEvent(iEventN);
     
-    Rich()->GetLoader()->MakeTree("R");  Rich()->MakeBranch("R");
-    Rich()->ResetDigits();  Rich()->ResetClusters();
+    R()->GetLoader()->MakeTree("R");  R()->MakeBranch("R");
+    R()->ResetDigits();               R()->ResetClusters();
     
-    Rich()->GetLoader()->TreeD()->GetEntry(0);
+    R()->GetLoader()->TreeD()->GetEntry(0);
     for(Int_t iChamber=1;iChamber<=kNCH;iChamber++){//chambers loop
       FindClusters(iChamber);
     }//chambers loop
-    Rich()->GetLoader()->TreeR()->Fill();
-    Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
+    R()->GetLoader()->TreeR()->Fill();  R()->GetLoader()->WriteRecPoints("OVERWRITE");//write out clusters for current event
   }//events loop  
-  Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();  
-  Rich()->ResetDigits();  Rich()->ResetClusters();
   
-  Rich()->GetLoader()->GetRunLoader()->UnloadHeader();
-  Rich()->GetLoader()->GetRunLoader()->UnloadKinematics();
+  R()->ResetDigits();//reset and unload everything
+  R()->ResetClusters();
+  R()->GetLoader()                ->UnloadDigits(); 
+  R()->GetLoader()                ->UnloadRecPoints();  
+  R()->GetLoader()->GetRunLoader()->UnloadHeader();
+  R()->GetLoader()->GetRunLoader()->UnloadKinematics();
 
   if(GetDebug()) Info("Exec","Stop.");      
 }//Exec()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FindClusters(Int_t iChamber)
 {
-  //finds neighbours and fill the tree with raw clusters
-  Int_t nDigits=Rich()->Digits(iChamber)->GetEntriesFast();
-//  Info("FindClusters","Start for Chamber %i with %i digits.",iChamber,nDigits);  
-  if(nDigits==0)return;
+//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);  
+  
+  if(iNdigits==0)return;//no digits for a given chamber, nothing to do
 
-  fHitMap=new AliRICHMap(Rich()->Digits(iChamber));//create digit map for the given chamber
+  fDigitMap=new AliRICHMap(R()->Digits(iChamber));//create digit map for the given chamber
 
-  for(Int_t iDig=0;iDig<nDigits;iDig++){    
-    AliRICHdigit *dig=(AliRICHdigit*)Rich()->Digits(iChamber)->At(iDig);
+  for(Int_t iDigN=0;iDigN<iNdigits;iDigN++){//digits loop for a given chamber    
+    AliRICHdigit *dig=(AliRICHdigit*)R()->Digits(iChamber)->At(iDigN);
     Int_t i=dig->X();   Int_t j=dig->Y();
-    if(fHitMap->TestHit(i,j)==kUsed) continue;
+    if(fDigitMap->TestHit(i,j)==kUsed) continue;//this digit is already taken, go after next digit
        
-    FormRawCluster(i,j);
-    FindLocalMaxima();
-//    cout << " fNlocals in FindCluster " << fNlocals << endl;
-    fRawCluster.CoG(fNlocals); // first initial approxmation of the CoG...to start minimization.
+    FormRawCluster(i,j);//form raw cluster starting from (i,j) pad 
+    if(GetDebug()) {Info("FindClusters","After FormRawCluster:");fRawCluster.Print();}  
+    FindLocalMaxima();  //find number of local maxima and initial center of gravity
+    if(GetDebug()) {Info("FindClusters","After FindLocalMaxima:");fRawCluster.Print();}  
     
-    if(AliRICHParam::IsResolveClusters()) {
-      ResolveCluster(); // ResolveCluster serialization will happen inside
-    } else {
-      WriteRawCluster(); // simply output of the RawCluster found without deconvolution
+    if(AliRICHParam::IsResolveClusters()&&fRawCluster.Size()>1&&fRawCluster.Size()<6){
+      FitCoG(); //serialization of resolved clusters will happen inside
+    }else{//cluster size=1 or resolving is switched off
+      WriteRawCluster();//simply output the formed raw cluster without deconvolution
     }
-    fRawCluster.Reset();
-    fResolvedCluster.Reset();
-  }//digits loop
+    fRawCluster.Reset(); fResolvedCluster.Reset();
+  }//digits loop for a given chamber
 
-  delete fHitMap;
-//  Info("FindClusters","Stop.");
+  delete fDigitMap;
   
+  if(GetDebug())Info("FindClusters","Stop.");  
 }//FindClusters()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FindClusterContribs(AliRICHcluster *pCluster)
 {
-// finds CombiPid for a given cluster
-// Info("FindClusterContribs","Start");
+//Finds cerenkov-feedback-mip mixture for a given cluster
+  if(GetDebug()) {Info("FindClusterContribs","Start.");pCluster->Print();}
   
   TObjArray *pDigits = pCluster->Digits();
+  if(!pDigits) return; //??????????
   Int_t iNmips=0,iNckovs=0,iNfeeds=0;
   TArrayI contribs(3*pCluster->Size());
   Int_t *pindex = new Int_t[3*pCluster->Size()];
   for(Int_t iDigN=0;iDigN<pCluster->Size();iDigN++) {//loop on digits of a given cluster
-    contribs[3*iDigN]  =((AliRICHdigit*)pDigits->At(iDigN))->Tid(0);
+    contribs[3*iDigN]  =((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(0);
     if (contribs[3*iDigN] >= 10000000) contribs[3*iDigN] = 0;
-    contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(1);
+    contribs[3*iDigN+1]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(1);
     if (contribs[3*iDigN+1] >= 10000000) contribs[3*iDigN+1] = 0;
-    contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->Tid(2);
+    contribs[3*iDigN+2]=((AliRICHdigit*)pDigits->At(iDigN))->GetTrack(2);
     if (contribs[3*iDigN+2] >= 10000000) contribs[3*iDigN+2] = 0;
   }//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 Tid
+  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);
     if(contribs[pindex[iDigN]]!=contribs[pindex[iDigN+1]]) {
-      TParticle* particle = Rich()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]);
+      TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[iDigN]]);
       Int_t code   = particle->GetPdgCode();
       Double_t charge = 0;
-      if (particle->GetPDG()) particle->GetPDG()->Charge();
+      if(particle->GetPDG()) charge=particle->GetPDG()->Charge();
+      if(GetDebug()) Info("FindClusterContribs"," charge of particle %f",charge);
 
       if(code==50000050) iNckovs++;
-      else if(code==50000051) iNfeeds++;
-      else if(charge!=0) iNmips++;
-      if (contribs[pindex[iDigN+1]]==kBad) break;
+      if(code==50000051) iNfeeds++;
+      if(charge!=0) iNmips++;
     }
   }//loop on digits to sort Tid
-  pCluster->SetCombiPid(iNckovs,iNfeeds,iNmips);
-//  pCluster->Print();
+  
+  if (contribs[pindex[3*pCluster->Size()-1]]!=kBad) {
+
+     TParticle* particle = R()->GetLoader()->GetRunLoader()->Stack()->Particle(contribs[pindex[3*pCluster->Size()-1]]);
+     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);
+     if(code==50000050) iNckovs++;
+     if(code==50000051) iNfeeds++;
+     if(charge!=0) iNmips++;
+  }
+    
+  pCluster->CFM(iNckovs,iNfeeds,iNmips);
+//  
   delete [] pindex; 
+  if(GetDebug()) {pCluster->Print();Info("FindClusterContribs","Stop.");}
 }//FindClusterContribs()
 //__________________________________________________________________________________________________
 void  AliRICHClusterFinder::FormRawCluster(Int_t i, Int_t j)
 {
-// Builder of the final Raw Cluster (before deconvolution)  
-  if(GetDebug()) Info("FormRawCluster","Start with digit(%i,%i)",i,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());
   
-  fRawCluster.AddDigit((AliRICHdigit*) fHitMap->GetHit(i,j));
-  fHitMap->FlagHit(i,j);// Flag hit as taken  
+  fRawCluster.AddDigit((AliRICHdigit*) fDigitMap->GetHit(i,j));//take this pad in cluster
+  fDigitMap->FlagHit(i,j);//flag this pad as taken  
 
   Int_t listX[4], listY[4];    //  Now look recursively for all neighbours
-  for (Int_t iNeighbour=0;iNeighbour<Rich()->P()->PadNeighbours(i,j,listX,listY);iNeighbour++)
-    if(fHitMap->TestHit(listX[iNeighbour],listY[iNeighbour])==kUnused) 
-                      FormRawCluster(listX[iNeighbour],listY[iNeighbour]);    
+  for (Int_t iNei=0;iNei<R()->P()->PadNeighbours(i,j,listX,listY);iNei++)
+    if(fDigitMap->TestHit(listX[iNei],listY[iNei])==kUnused) FormRawCluster(listX[iNei],listY[iNei]);    
 }//FormRawCluster()
 //__________________________________________________________________________________________________
-void AliRICHClusterFinder::ResolveCluster()
-{// Decluster algorithm
-  if(GetDebug()) {Info("ResolveCluster","Start."); fRawCluster.Print();}
-  
-  switch (fRawCluster.Size()) {
-  
-  case 1:                     // nothing to decluster: cluster size = 1
-    WriteRawCluster();break; 
-  default:                     // cluster size > 1: if=2 FitCoG; if>2 Resolve and FitCoG
-    FitCoG();break;
-  }     
-}//ResolveCluster()
+void AliRICHClusterFinder::FindLocalMaxima()
+{
+//find number of local maxima in the current raw cluster and then calculates initial center of gravity
+  fNlocals=0;
+  for(Int_t iDig1=0;iDig1<fRawCluster.Size();iDig1++) {
+    Int_t iNotMax = 0;
+    AliRICHdigit *pDig1 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig1);
+    TVector pad1 = pDig1->Pad();
+    Int_t padQ1 = (Int_t)(pDig1->Q()+0.1);
+    Int_t padC1 = pDig1->ChFbMi();
+    for(Int_t iDig2=0;iDig2<fRawCluster.Size();iDig2++) {
+      AliRICHdigit *pDig2 = (AliRICHdigit *)fRawCluster.Digits()->At(iDig2);
+      TVector pad2 = pDig2->Pad();
+      Int_t padQ2 = (Int_t)(pDig2->Q()+0.1);
+      if(iDig1==iDig2) continue;
+      Int_t diffx = TMath::Sign(Int_t(pad1[0]-pad2[0]),1);
+      Int_t diffy = TMath::Sign(Int_t(pad1[1]-pad2[1]),1);
+      if((diffx+diffy)<=1) {
+         if(padQ2>=padQ1) iNotMax++;
+      }
+    }
+    if(iNotMax==0) {
+      TVector2 x2=AliRICHParam::Pad2Loc(pad1);
+      fLocalX[fNlocals]=x2.X();fLocalY[fNlocals]=x2.Y();
+      fLocalQ[fNlocals] = (Double_t)padQ1;
+      fLocalC[fNlocals] = padC1;
+      fNlocals++;
+    }
+  }
+  fRawCluster.CoG(fNlocals); //first initial approximation of the CoG...to start minimization.
+}//FindLocalMaxima()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::WriteRawCluster()
 {
-// out the current raw cluster
+//Add the current raw cluster to the list of clusters
   if(GetDebug()) Info("WriteRawCluster","Start.");
   
-  FindClusterContribs(&fRawCluster);
-  if(GetDebug()) fRawCluster.Dump();
-  Rich()->AddCluster(fRawCluster);
-//  fRawCluster.Print();
+  FindClusterContribs(&fRawCluster);  
+  R()->AddCluster(fRawCluster);
+  
+  if(GetDebug()){fRawCluster.Print(); Info("WriteRawCluster","Stop.");}  
 }//WriteRawCluster()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::WriteResolvedCluster()
 {
-// out the current resolved cluster
+//Add the current resolved cluster to the list of clusters
   if(GetDebug()) Info("WriteResolvedCluster","Start.");
   
-//  FindClusterContribs(&fResolvedCluster);
-  Rich()->AddCluster(fResolvedCluster);
+  FindClusterContribs(&fResolvedCluster);  
+  R()->AddCluster(fResolvedCluster);
   
+  if(GetDebug()){fResolvedCluster.Print(); Info("WriteResolvedCluster","Stop.");}  
 }//WriteResolvedCluster()
 //__________________________________________________________________________________________________
 void AliRICHClusterFinder::FitCoG()
-{// Fit cluster size 2 by single Mathieson
-  if(GetDebug()) Info("FitCoG","Start with size %3i and local maxima %3i",fRawCluster.Size(),fNlocals);
+{
+//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();}
   
   Double_t arglist;
   Int_t ierflag = 0;
 
-//  fRawCluster.Print();
-//  if(fNlocals==0) fRawCluster.Print();
-  if(fNlocals==0||fNlocals>3) {WriteRawCluster();return;}
-  
   TMinuit *pMinuit = new TMinuit(3*fNlocals-1);
   pMinuit->mninit(5,10,7);
   
@@ -265,19 +255,15 @@ void AliRICHClusterFinder::FitCoG()
     lower    = vstart - 2*AliRICHParam::PadSizeX();
     upper    = vstart + 2*AliRICHParam::PadSizeX();
     pMinuit->mnparm(3*i  ,Form("xCoG  %i",i),vstart,stepX,lower,upper,ierflag);
-//    cout << Form("xCoG  %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl;
     vstart   = fLocalY[i];
     lower    = vstart - 2*AliRICHParam::PadSizeY();
     upper    = vstart + 2*AliRICHParam::PadSizeY();
     pMinuit->mnparm(3*i+1,Form("yCoG  %i",i),vstart,stepY,lower,upper,ierflag);
-//    cout << Form("yCoG  %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl;
     if(i==fNlocals-1) break;                    // last parameter is constrained
     vstart = fLocalQ[i]/fRawCluster.Q();
     lower  = 0;
     upper  = 1;
     pMinuit->mnparm(3*i+2,Form("qfrac %i",i),vstart,stepQ,lower,upper,ierflag);
-//    cout << Form("qCoG  %i",i) << "start " << vstart << "lower " << lower << "upper " << upper << endl;
-
   }
   
   arglist = -1;
@@ -303,23 +289,16 @@ void AliRICHClusterFinder::FitCoG()
   qfracCoG[fNlocals-1] = 1 - qfraclast;
   delete pMinuit;
 
-  for(Int_t i=0;i<fNlocals;i++) {
-
-    if(fNlocals==5) {
-    cout << " Start writing  deconvolved cluster n." << i << endl;
-    cout << " fRawCluster " << &fRawCluster << " xCoG " << xCoG[i] << " yCoG " << yCoG[i] << " qfrac " << qfracCoG[i] << endl;
-    cout << " Combipid " << fLocalC[i] << " for maximum n. " << i << endl; 
-  }
+  for(Int_t i=0;i<fNlocals;i++){//resolved positions loop
     fResolvedCluster.Fill(&fRawCluster,xCoG[i],yCoG[i],qfracCoG[i],fLocalC[i]);
-//    fResolvedCluster.Print();
     WriteResolvedCluster();
   }
-if(fNlocals==5)  Info("CoG","Stop.");
+  if(GetDebug()) Info("CoG","Stop.");
 }//FitCoG()
 //__________________________________________________________________________________________________
 void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, Int_t )
 {
-// Mathieson minimization function 
+//Mathieson minimization function 
   
   AliRICHcluster *pRawCluster = ((AliRICHClusterFinder*)gMinuit->GetObjectFit())->GetRawCluster();
 
@@ -338,19 +317,13 @@ void RICHMinMathieson(Int_t &npar, Double_t *, Double_t &chi2, Double_t *par, In
   chi2 = 0;
   Int_t qtot = pRawCluster->Q();
   for(Int_t i=0;i<pRawCluster->Size();i++) {
-    Int_t    padX = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->X();
-    Int_t    padY = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Y();
+    TVector  pad=((AliRICHdigit *)pRawCluster->Digits()->At(i))->Pad();
     Double_t padQ = ((AliRICHdigit *)pRawCluster->Digits()->At(i))->Q();
     Double_t qfracpar=0;
     for(Int_t j=0;j<nFunctions;j++) {
-      qfracpar += q[j]*AliRICHParam::FracQdc(centroid[j],padX,padY);
-//      cout << " function n. " << j+1 << " q " << q[j] << " fracMat " << AliRICHParam::FracQdc(centroid[j],padX,padY) 
-//           << " xpar " << centroid[j].X() << " ypar " << centroid[j].Y() << endl;
+      qfracpar += q[j]*AliRICHParam::FracQdc(centroid[j],pad);
     }
     chi2 += TMath::Power((qtot*qfracpar-padQ),2)/padQ;
-//    cout << " chi2 " << chi2 << " pad n. " << i << " qfracpar " << qfracpar << endl;
   }     
-//  if(iflag==3) {
-//    cout << " chi2 final " << chi2 << endl;
-//  }
 }//RICHMinMathieson()
+//__________________________________________________________________________________________________
index 615e39f..e9a5598 100644 (file)
@@ -5,8 +5,6 @@
  * See cxx source for full Copyright notice                               */
 
 #include "TTask.h"
-#include "TArrayD.h"
-#include "TArrayI.h"
 
 #include "AliRICH.h"
 class AliHitMap;
@@ -17,27 +15,26 @@ public:
            AliRICHClusterFinder(AliRICH *pRICH);
   virtual ~AliRICHClusterFinder()                                          {;}
   
-  AliRICH *Rich() {return fRICH;}                                             //Pointer to RICH  
-  void     Exec();                                                            //Loop on events and chambers  
-  void     FindClusters(Int_t iChamber);                                      //Find all clusters for a given chamber
-  void     FindClusterContribs(AliRICHcluster *pCluster);                     //Find CombiPid for the current cluster
-  void     FormRawCluster(Int_t i, Int_t j);                                  //form a raw cluster
-  void     FindLocalMaxima();                                                 //Find local maxima in a cluster
-  void     ResolveCluster();                                                  //Try to resolve a cluster with maxima > 2
-  void     FitCoG();                                                          //Evaluate the CoG as the best 
-  void     WriteRawCluster();                                                 //write in the list of cluster  
-  void     WriteResolvedCluster();                                            //write in the list of cluster  
-  AliRICHcluster *GetRawCluster() {return &fRawCluster;}                      //Return pointer to the current raw cluster
-  Bool_t   GetDebug()            const {return fRICH->GetDebug();}            //is debug printout needed?
+  AliRICH        *R()                                              {return fRICH;}   //returns pointer to RICH  
+  void            Exec();                                                            //loop on events and chambers  
+  void            FindClusters(Int_t iChamber);                                      //find all clusters for a given chamber
+  void            FindClusterContribs(AliRICHcluster *pCluster);                     //find CFM for the current cluster
+  void            FormRawCluster(Int_t i, Int_t j);                                  //form a raw cluster
+  void            FindLocalMaxima();                                                 //find local maxima in a cluster
+  void            FitCoG();                                                          //evaluate the CoG as the best 
+  void            WriteRawCluster();                                                 //write in the list of cluster  
+  void            WriteResolvedCluster();                                            //write in the list of cluster  
+  AliRICHcluster *GetRawCluster()                           {return &fRawCluster;}   //returns pointer to the current raw cluster
+  Bool_t          GetDebug()                      const{return fRICH->GetDebug();}   //is debug printout needed?
 protected:
-  AliRICH                *fRICH;                         //Pointer to RICH
-  AliHitMap              *fHitMap;                       //Hit Map with digit positions
-  AliRICHcluster         fRawCluster;                    //Current raw cluster before deconvolution
-  AliRICHcluster         fResolvedCluster;               //Current cluster after deconvolution
-  Int_t                  fNlocals;                       // number of local maxima
-  TArrayD                fLocalX,fLocalY;                // list of locals X,Y
-  TArrayD                fLocalQ;                        // list of locals charge Q
-  TArrayI                fLocalC;                        // list of locals CombiPid
-  ClassDef(AliRICHClusterFinder,0) //Finds raw clusters, trasfers them to resolved clusters through declustering.
+  AliRICH                *fRICH;                         //pointer to RICH
+  AliHitMap              *fDigitMap;                     //map of digits positions
+  AliRICHcluster         fRawCluster;                    //current raw cluster before deconvolution
+  AliRICHcluster         fResolvedCluster;               //current cluster after deconvolution
+  Int_t                  fNlocals;                       //number of local maxima
+  Double_t               fLocalX[100],fLocalY[100];      //list of locals X,Y
+  Double_t               fLocalQ[100];                   //list of locals charge Q
+  Int_t                  fLocalC[100];                   //list of locals CombiPid
+  ClassDef(AliRICHClusterFinder,0) //finds raw clusters, trasfers them to resolved clusters through declustering.
 };
 #endif
index 3eaa280..dbca6d8 100644 (file)
@@ -36,16 +36,15 @@ Bool_t AliRICHDigitizer::Init()
 {
 //This methode is called from AliRunDigitizer after the corresponding file is open
   if(GetDebug())Info("Init","Start.");
-  AliRunLoader *pOutAL = 
-    AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
+  AliRunLoader *pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
   if (!pOutAL->GetAliRun()) pOutAL->LoadgAlice();
-  fRich=(AliRICH*)pOutAL->GetAliRun()->GetDetector("RICH");
-  Rich()->P()->GenSigmaThMap();
+  fRich=(AliRICH*)pOutAL->GetAliRun()->GetDetector("RICH");//retrive RICH pointer from OUTPUT stream
   return kTRUE;
 }//Init()
 //__________________________________________________________________________________________________
 void AliRICHDigitizer::Exec(Option_t*)
 {
+//this method invoked   
   if(GetDebug())Info("Exec","Start with %i input(s) for event %i",fManager->GetNinputs(),fManager->GetOutputEventNr());
   
   AliRunLoader *pInAL=0, *pOutAL;//in and out Run loaders
@@ -53,47 +52,52 @@ void AliRICHDigitizer::Exec(Option_t*)
  
   pOutAL = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
   pOutRL = pOutAL->GetLoader("RICHLoader");
-  pOutRL->MakeTree("D");   Rich()->MakeBranch("D"); //create TreeD with RICH branches in output stream
+  pOutRL->MakeTree("D");   R()->MakeBranch("D"); //create TreeD with RICH branches in output stream
   
   TClonesArray tmpCA("AliRICHdigit");//tmp storage for sdigits sum up from all input files
   Int_t total=0;
   for(Int_t inFileN=0;inFileN<fManager->GetNinputs();inFileN++){//files loop
-    pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inFileN)); 
+    pInAL = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(inFileN));//get run loader from current input 
     pInRL = pInAL->GetLoader("RICHLoader"); if(pInRL==0) continue;//no RICH in this input, check the next input
     if (!pInAL->GetAliRun()) pInAL->LoadgAlice();
     AliRICH* rich=(AliRICH*)pInAL->GetAliRun()->GetDetector("RICH");
     pInRL->LoadSDigits(); pInRL->TreeS()->GetEntry(0);
-    Info("Exec","input %i has %i sdigits",inFileN,rich->SDigits()->GetEntries());
-    for(Int_t i=0;i<rich->SDigits()->GetEntries();i++) {
+    if(GetDebug())Info("Exec","input %i has %i sdigits",inFileN,rich->SDigits()->GetEntries());
+    for(Int_t i=0;i<rich->SDigits()->GetEntries();i++){//collect sdigits from current input to tmpCA
       new(tmpCA[total++]) AliRICHdigit(*(AliRICHdigit*)rich->SDigits()->At(i)); 
       ((AliRICHdigit*)tmpCA[total-1])->AddTidOffset(fManager->GetMask(inFileN));
     }
     pInRL->UnloadSDigits();   rich->ResetSDigits();
   }//files loop
   
-  tmpCA.Sort();                     //sort them according to Id() methode
+  tmpCA.Sort();                     //sort them according to Id() method
   
-  Int_t chFbMip=0,chamber=0,x=0,y=0,tid[3],id=0; Double_t q=0;
+  if(GetDebug()) {tmpCA.Print();Info("Exec","Totally %i sdigits in %i inputs",tmpCA.GetEntries(),fManager->GetNinputs());}
+  Int_t cfm=0,chamber=0,id=0; //cfm is cerenkov feedback mip mixture
+  TVector pad(2); pad[0]=0;pad[1]=0;
+  Double_t q=0;
+  Int_t tid[3]={0,0,0};
   Int_t iNdigitsPerPad=0;//how many sdigits for a given pad
   for(Int_t i=0;i<tmpCA.GetEntries();i++){//sdigits loop (sorted)
     AliRICHdigit *pSdig=(AliRICHdigit*)tmpCA.At(i);//get new sdigit
     if(pSdig->Id()==id){//still the same pad
-      iNdigitsPerPad++;         q+=pSdig->Q();       chFbMip+=pSdig->ChFbMi();//sum up charge and cfm
-      if(iNdigitsPerPad<=3)        tid[iNdigitsPerPad-1]=pSdig->Tid(0);
+      iNdigitsPerPad++;         q+=pSdig->Q();       cfm+=pSdig->ChFbMi();//sum up charge and cfm
+      if(pSdig->ChFbMi()==1) tid[0] = pSdig->GetTrack(0); // force the first tid to be mip's tid if it exists in the current pad
+      if(iNdigitsPerPad<=3)        tid[iNdigitsPerPad-1]=pSdig->GetTrack(0);
       else                         if(GetDebug())Warning("Exec","More then 3 sdigits for the given pad");
     }else{//new pad, add the pevious one
-        if(id!=kBad&&Rich()->P()->IsOverTh(chamber,x,y,q)) Rich()->AddDigit(chamber,x,y,(Int_t)q,chFbMip,tid); //add newly created dig
-        chFbMip=pSdig->ChFbMi(); chamber=pSdig->C(); id=pSdig->Id();  x=pSdig->X(); y=pSdig->Y(); q=pSdig->Q();  //init all values by current sdig
-        iNdigitsPerPad=1; tid[0]=pSdig->Tid(0); tid[1]=tid[2]=kBad;
+        if(id!=kBad&&R()->P()->IsOverTh(chamber,pad,q)) R()->AddDigit(chamber,pad,(Int_t)q,cfm,tid); //add newly created dig
+        cfm=pSdig->ChFbMi(); chamber=pSdig->C(); id=pSdig->Id();  pad=pSdig->Pad(); q=pSdig->Q();  //init all values by current sdig
+        iNdigitsPerPad=1; tid[0]=pSdig->GetTrack(0); tid[1]=tid[2]=kBad;
       }
   }//sdigits loop (sorted)
-  if(tmpCA.GetEntries()&&Rich()->P()->IsOverTh(chamber,x,y,q)) Rich()->AddDigit(chamber,x,y,(Int_t)q,chFbMip,tid);//add the last dig
+  if(tmpCA.GetEntries()&&R()->P()->IsOverTh(chamber,pad,q)) R()->AddDigit(chamber,pad,(Int_t)q,cfm,tid);//add the last dig
   
   pOutRL->TreeD()->Fill();              //fill the tree with the list of digits
   pOutRL->WriteDigits("OVERWRITE");     //serialize them to file
   
   tmpCA.Clear();
-  pOutRL->UnloadDigits();   Rich()->ResetDigits();
+  pOutRL->UnloadDigits();   R()->ResetDigits();
             
   if(GetDebug())Info("Exec","Stop.");
 }//Exec()
index c907b65..d6b24ff 100644 (file)
@@ -13,7 +13,7 @@ class AliRICH;
 class AliRICHDigitizer : public AliDigitizer 
 {
 public:
-           AliRICHDigitizer()                                                {;}
+           AliRICHDigitizer():fRich(0)                                       {;}
            AliRICHDigitizer(AliRunDigitizer * manager):AliDigitizer(manager) {if(GetDebug())Info("main ctor","Start.");}
   virtual ~AliRICHDigitizer()                                                {if(GetDebug())Info("dtor","Start.");}
 
@@ -21,7 +21,7 @@ public:
   void     Exec(Option_t* option=0);                //virtual
   Bool_t   Init();                                  //virtual
   Bool_t   GetDebug() const {return ((gAlice) ? gAlice->GetDebug() : kFALSE);}
-  AliRICH* Rich()     const {return fRich;}
+  AliRICH* R()        const {return fRich;}         //returns pointer to RICH
 protected:
   AliRICH* fRich; //pointer to main RICH object
   ClassDef(AliRICHDigitizer,0)
index 610e55f..656872f 100644 (file)
@@ -28,7 +28,7 @@
 ClassImp(AliRICHDisplFast)
 
 //__________________________________________________________________________________________________
-void AliRICHDisplFast::Exec()
+void AliRICHDisplFast::Exec(Option_t *)
 {
   AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH");
   Bool_t isHits    =!pRich->GetLoader()->LoadHits();
@@ -42,14 +42,14 @@ void AliRICHDisplFast::Exec()
   
   TH2F *pHitsH2=0,*pDigitsH2=0,*pClustersH2=0;
   
-  if(isHits)     pHitsH2     = new TH2F("pHitsH2"  ,  "Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2,
-                                                             144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2);
+  if(isHits)     pHitsH2     = new TH2F("pHitsH2"  ,  "Event Display;x,cm;y,cm",165,0,AliRICHParam::PcSizeX(),
+                                                                      144,0,AliRICHParam::PcSizeY());
   if(pHitsH2)    pHitsH2->SetStats(kFALSE);
   
-  if(isDigits)   pDigitsH2   = new TH2F("pDigitsH2"  ,"Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2,
-                                                             144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2);
-  if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,-AliRICHParam::PcSizeX()/2,AliRICHParam::PcSizeX()/2,
-                                                             144,-AliRICHParam::PcSizeY()/2,AliRICHParam::PcSizeY()/2);
+  if(isDigits)   pDigitsH2   = new TH2F("pDigitsH2"  ,"Event Display",165,0,AliRICHParam::PcSizeX(),
+                                                                      144,0,AliRICHParam::PcSizeY());
+  if(isClusters) pClustersH2 = new TH2F("pClustersH2","Event Display",165,0,AliRICHParam::PcSizeX(),
+                                                                      144,0,AliRICHParam::PcSizeY());
 
     
 
@@ -101,7 +101,7 @@ void AliRICHDisplFast::Exec()
         pRich->GetLoader()->TreeD()->GetEntry(0);
         for(Int_t j=0;j<pRich->Digits(iChamber)->GetEntries();j++){//digits loop
           AliRICHdigit *pDig = (AliRICHdigit*)pRich->Digits(iChamber)->At(j);
-         TVector2 x2=AliRICHParam::Pad2Loc(pDig->X(),pDig->Y());
+         TVector2 x2=AliRICHParam::Pad2Loc(pDig->Pad());
          pDigitsH2->Fill(x2.X(),x2.Y(),100);
         }//digits loop
         pDigitsH2->SetMarkerColor(kGreen); pDigitsH2->SetMarkerStyle(29); pDigitsH2->SetMarkerSize(0.4);
@@ -133,70 +133,19 @@ void AliRICHDisplFast::Exec()
 //__________________________________________________________________________________________________
 void AliRICHDisplFast::DrawSectors() 
 { 
-  Double_t x1[5] = {-AliRICHParam::PcSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(),
-                    -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(),
-                    -AliRICHParam::PcSizeX()/2,
-                    -AliRICHParam::PcSizeX()/2};
-  Double_t y1[5] = {AliRICHParam::DeadZone()/2,
-                    AliRICHParam::DeadZone()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::DeadZone()/2};  
-  Double_t x2[5] = {-AliRICHParam::SectorSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2};
-  Double_t y2[5] = {AliRICHParam::DeadZone()/2,
-                    AliRICHParam::DeadZone()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::DeadZone()/2};
-  Double_t x3[5] = { AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(),
-                     AliRICHParam::PcSizeX()/2,
-                     AliRICHParam::PcSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(),
-                     AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone()};
-  Double_t y3[5] = {AliRICHParam::DeadZone()/2,
-                    AliRICHParam::DeadZone()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::PcSizeY()/2,
-                    AliRICHParam::DeadZone()/2};
-  Double_t x4[5] = {-AliRICHParam::PcSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(),
-                    -AliRICHParam::SectorSizeX()/2-AliRICHParam::DeadZone(),
-                    -AliRICHParam::PcSizeX()/2,
-                    -AliRICHParam::PcSizeX()/2};
-  Double_t y4[5] = {-AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::PcSizeY()/2};
-   Double_t x5[5] = {-AliRICHParam::SectorSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2,
-                    -AliRICHParam::SectorSizeX()/2};
-  Double_t y5[5] = {-AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::PcSizeY()/2};
-  Double_t x6[5] = { AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(),
-                     AliRICHParam::PcSizeX()/2,
-                     AliRICHParam::PcSizeX()/2,
-                     AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone(),
-                     AliRICHParam::SectorSizeX()/2+AliRICHParam::DeadZone()};
-  Double_t y6[5] = {-AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::PcSizeY()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::DeadZone()/2,
-                    -AliRICHParam::PcSizeY()/2};
-  TPolyLine *sector1 = new TPolyLine(5,x1,y1);  sector1->SetLineColor(21);  sector1->Draw();
-  TPolyLine *sector2 = new TPolyLine(5,x2,y2);  sector2->SetLineColor(21);  sector2->Draw();
-  TPolyLine *sector3 = new TPolyLine(5,x3,y3);  sector3->SetLineColor(21);  sector3->Draw();
-  TPolyLine *sector4 = new TPolyLine(5,x4,y4);  sector4->SetLineColor(21);  sector4->Draw();
-  TPolyLine *sector5 = new TPolyLine(5,x5,y5);  sector5->SetLineColor(21);  sector5->Draw();
-  TPolyLine *sector6 = new TPolyLine(5,x6,y6);  sector6->SetLineColor(21);  sector6->Draw();
+  AliRICHParam p;
+  Double_t xLeft[5]  = {0,0,p.SectorSizeX(),p.SectorSizeX(),0};
+  Double_t xRight[5] = {p.SectorSizeX()+p.DeadZone(),p.SectorSizeX()+p.DeadZone(),p.PcSizeX(),p.PcSizeX(),p.SectorSizeX()+p.DeadZone()};
+  
+  Double_t yDown[5]   = {0,p.SectorSizeY(),p.SectorSizeY(),0,0};
+  Double_t yCenter[5] = {  p.SectorSizeY()+p.DeadZone(),2*p.SectorSizeY()+p.DeadZone(),2*p.SectorSizeY()+p.DeadZone(),
+                           p.SectorSizeY()+p.DeadZone(),p.SectorSizeY()+p.DeadZone()};  
+  Double_t yUp[5]     = {2*p.SectorSizeY()+2*p.DeadZone(),p.PcSizeY(),p.PcSizeY(),2*p.SectorSizeY()+2*p.DeadZone(),2*p.SectorSizeY()+2*p.DeadZone()};
+  
+  TPolyLine *sec1 = new TPolyLine(5,xLeft ,yDown);    sec1->SetLineColor(21);  sec1->Draw();
+  TPolyLine *sec2 = new TPolyLine(5,xRight,yDown);    sec2->SetLineColor(21);  sec2->Draw();
+  TPolyLine *sec3 = new TPolyLine(5,xLeft ,yCenter);  sec3->SetLineColor(21);  sec3->Draw();
+  TPolyLine *sec4 = new TPolyLine(5,xRight,yCenter);  sec4->SetLineColor(21);  sec4->Draw();
+  TPolyLine *sec5 = new TPolyLine(5,xLeft, yUp);      sec5->SetLineColor(21);  sec5->Draw();
+  TPolyLine *sec6 = new TPolyLine(5,xRight,yUp);      sec6->SetLineColor(21);  sec6->Draw();
 }//DrawSectors()
index 9d4ef00..f46ccd6 100644 (file)
@@ -13,10 +13,10 @@ class AliRICH;
 class AliRICHDisplFast : public TTask 
 {
 public :
-              AliRICHDisplFast() {;}
-  virtual    ~AliRICHDisplFast() {;}      
-  static void DrawSectors();          //Draw sectors in plot 
-  void        Exec();                 //virtual do the main job
+               AliRICHDisplFast() {;}
+  virtual     ~AliRICHDisplFast() {;}      
+  static  void DrawSectors();          //Draw sectors in plot 
+  virtual void Exec(Option_t *opt=0); //virtual do the main job
 protected:  
   ClassDef(AliRICHDisplFast,0)        //Utility class to draw the current event topology
 };
index f8d0a24..ad1679f 100644 (file)
 //  * provided "as is" without express or implied warranty.                  *
 //  **************************************************************************
 #include "AliRICHParam.h"
-#include <Riostream.h>
+#include "AliRICHChamber.h"
 
 ClassImp(AliRICHParam)
 Bool_t   AliRICHParam::fgIsWireSag            =kTRUE;
 Bool_t   AliRICHParam::fgIsResolveClusters    =kTRUE;
+Bool_t   AliRICHParam::fgIsRadioSrc           =kFALSE;
 Double_t AliRICHParam::fgAngleRot             =-60;
-Int_t    AliRICHParam::fgHV[kNsectors]        ={2150,2100,2050,2000,2150,2150};
+Int_t    AliRICHParam::fgHV[kNsectors]        ={2050,2050,2050,2050,2050,2050};
 Int_t    AliRICHParam::fgNsigmaTh             =4;
-Float_t  AliRICHParam::fgSigmaThMean          =1.5;
-Float_t  AliRICHParam::fgSigmaThSpread        =0.5;      
-Float_t  AliRICHParam::fSigmaThMap[kNCH][kNpadsX][kNpadsY];
+Float_t  AliRICHParam::fgSigmaThMean          =1.132; //QDC 
+Float_t  AliRICHParam::fgSigmaThSpread        =0.035; //     
 
-void AliRICHParam::GenSigmaThMap()
-{
-// Generate the map of thresholds sigmas for all pads of all chambers 
-  for(Int_t iChamber=0;iChamber<kNCH;iChamber++)
-    for(Int_t ipadX=0;ipadX<NpadsX();ipadX++)
-      for(Int_t ipadY=0;ipadY<NpadsY();ipadY++) 
-        fSigmaThMap[iChamber][ipadX][ipadY] = SigmaThMean()+(1.-2*gRandom->Rndm())*SigmaThSpread();
-  //  Info("GenSigmaThMap"," Threshold map generated for all RICH chambers");
-}
 //__________________________________________________________________________________________________
-void AliRICHParam::Print()
+void AliRICHParam::Print(Option_t*)
 {
-  cout<<"\nPads in chamber ("<<NpadsX()<<','<<NpadsY()<<") in sector ("<<NpadsXsec()<<','<<NpadsYsec()<<')'<<endl;
-  cout<<"PC size ("<<PcSizeX()<<','<<PcSizeY()<<") sector size ("<<SectorSizeX()<<','<<SectorSizeY()<<") pad size ("<<PadSizeX()<<','
-      <<PadSizeY()<<") Dead zone "<<DeadZone()<<endl;
-  cout<<"Anode wire pitch "<<WirePitch()<<" Anode-Cathode gap "<<AnodeCathodeGap()<<" Protection wires-cathode gap "<<ProximityGap()<<endl<<endl;
+  ::Info("","Pads in chamber (%3i,%3i) in sector (%2i,%2i)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec());
+  fpChambers->Print();
 }
 //__________________________________________________________________________________________________
+void AliRICHParam::CreateChambers()
+{
+//Create all RICH Chambers on each call. Previous chambers deleted.
+  if(fpChambers) delete fpChambers;
+  fpChambers=new TObjArray(kNchambers);
+  fpChambers->SetOwner();
+  for(int iChamberN=0;iChamberN<kNchambers;iChamberN++)  fpChambers->AddAt(new AliRICHChamber(iChamberN+1),iChamberN);  
+}//void AliRICH::CreateChambers()
index 9d0ac82..958e90c 100644 (file)
@@ -4,38 +4,43 @@
 #include <TObject.h>
 #include <TMath.h>
 #include <TVector2.h>
+#include <TVector3.h>
 #include <TRandom.h>
 #include <TError.h>
+#include <TObjArray.h>
 
 
-static const int kNCH=7;           //number of RICH chambers 
-static const int kNpadsX = 144;    //number of pads along X in single chamber
-static const int kNpadsY = 160;    //number of pads along Y in single chamber
+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;      // nb. of sectors per chamber
+static const int kNsectors=6;      //number of sectors per chamber
 
-static const int kadc_satm  = 4096;  //dynamic range (10 bits)
+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
 
+class AliRICHChamber;
 
 class AliRICHParam :public TObject  
 {
 public:
-           AliRICHParam()                    {;}
-  virtual ~AliRICHParam()                    {;}
-  static const Int_t   NpadsX()              {return kNpadsX;}                           //pads along X in chamber
-  static const Int_t   NpadsY()              {return kNpadsY;}                           //pads along Y in chamber
-  static Int_t    NpadsXsec()                {return NpadsX()/3;}                        //pads along X in sector
-  static Int_t    NpadsYsec()                {return NpadsY()/2;}                        //pads along Y in sector
+                  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
+  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
+  static Int_t    NpadsYsec()                {return NpadsY()/3;}                        //pads along Y in sector
   static Double_t DeadZone()                 {return 2.6;}                               //dead zone size in cm  
-  static Double_t PadSizeX()                 {return 0.84;}                              //pad size x in cm 
-  static Double_t PadSizeY()                 {return 0.8;}                               //pad size y in cm   
-  static Double_t SectorSizeX()              {return NpadsX()*PadSizeX()/3;}             //sector size x in cm
-  static Double_t SectorSizeY()              {return NpadsY()*PadSizeY()/2;}             //sector size y in cm 
-  static Double_t PcSizeX()                  {return NpadsX()*PadSizeX()+2*DeadZone();}  //photocathode size x in cm
-  static Double_t PcSizeY()                  {return NpadsY()*PadSizeY()+DeadZone();}    //photocathode size y in cm 
-  static Double_t WirePitch()                {return PadSizeX()/2;}                      //distance between anode wires
+  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;}                             
@@ -45,23 +50,32 @@ public:
   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 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 GapThickness()             {return 8.0;}      
   static Double_t RadiatorToPads()           {return FreonThickness()+QuartzThickness()+GapThickness();}   
-  static Double_t ProximityGap()             {return 0.445;}    
-  static Double_t AnodeCathodeGap()          {return 0.2;}
+  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;}                            
-  static TVector2 MathiesonDelta()           {return TVector2(5*0.18,5*0.18);}    
-  static Int_t    MaxQdc()                   {return 4095;}          
-  static Double_t AlphaFeedback(Int_t sec)   {HV(sec);return 0.036;}
-  
+  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?
+  static Bool_t   IsRadioSrc()                {return fgIsRadioSrc;}         //add radioactive source inside CH4?
   static Int_t    HV(Int_t sector)            {
     if (sector>=1 && sector <=6)
       return fgHV[sector-1];
@@ -70,188 +84,193 @@ public:
       return kBad;
     } 
   }       //high voltage for this sector
-  static void     IsResolveClusters(Bool_t a) {fgIsResolveClusters=a;}  
+  static void     SetDeclustering(Bool_t a)   {fgIsResolveClusters=a;}  
+  static void     SetRadioSrc(Bool_t a)       {fgIsRadioSrc=a;}  
   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;}
 
-  inline static void     Loc2Area(TVector2 x2,Int_t &padxMin,Int_t &padyMin,Int_t &padxMax,Int_t &padyMax); //
-  inline static Int_t    Loc2Pad(TVector2 x2,Int_t &padx,Int_t &pady);                             //return sector and pad
-  inline static TVector2 Pad2Loc(Int_t padx,Int_t pady);                                           //return center of the pad
-         static Int_t    Sector(Int_t padx,Int_t pady)          {return Pad2Sec(padx,pady);}       //sector of this pad
-         static Int_t    Sector(TVector2 x2)                    {int x,y;return Loc2Pad(x2,x,y);}  //sector of this point
+  inline static TVector  Loc2Area(TVector2 x2);                                                    //return area affected by hit x2
+  inline static TVector  Loc2Pad(TVector2 x2);                                                     //return pad containing given position
+  inline static TVector2 Pad2Loc(TVector pad);                                                     //return center of the pad
+         static TVector2 Pad2Loc(Int_t x,Int_t y) {TVector pad(2);pad[0]=x;pad[1]=y;return Pad2Loc(pad);}
   inline static Int_t    PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t aListX[4],Int_t aListY[4]);   //number of neighbours for this pad
-  inline static TVector2 ShiftToWirePos(TVector2 x2);                                              //shift to the nearest wire
   
-  inline static Double_t Mathieson(Double_t lx1,Double_t lx2,Double_t ly1,Double_t ly2);           //Mathienson integral over these limits
-  inline static Double_t GainSag(Double_t y,Int_t sector);                                         //gain variations in %
-  inline static Double_t QdcSlope(Int_t sec);                                                      //weight of electon in QDC channels
-  inline static Double_t Gain(TVector2 x2);                                                        //gain for point in ChRS 
-  inline static Double_t FracQdc(TVector2 x2,Int_t padx,Int_t pady);                               //charge fraction to pad from hit
+  inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2);               //Mathienson integral over these limits
+  inline static Double_t GainSag(Double_t x,Int_t sector);                                         //gain variations in %
+         static Double_t QdcSlope(Int_t sec){switch(sec){case kBad: return 0;  default:   return 33;}} //weight of electon in QDC channels
+         static Double_t Gain(TVector2 x2){if(IsWireSag()) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);else return QdcSlope(Loc2Sec(x2));}//gain for point in chamber RS 
+  inline static Double_t FracQdc(TVector2 x2,TVector pad);                                         //charge fraction to pad from hit
   inline static Int_t    TotQdc(TVector2 x2,Double_t eloss);                                       //total charge for hit eloss=0 for photons
-  inline        Bool_t   IsOverTh(Int_t iChamber, Int_t x, Int_t y, Double_t q);                   //  
-         static Int_t   NsigmaTh()                    {return fgNsigmaTh;}                         //
-         static Float_t SigmaThMean()                 {return fgSigmaThMean;}      //
-         static Float_t SigmaThSpread()               {return fgSigmaThSpread;}    //
-                void    GenSigmaThMap();                                           //generate pedestal map
-         static void    Print();                
-protected:
+  inline        Bool_t   IsOverTh(Int_t c,TVector pad,Double_t q);                                 //is QDC of the pad registered by FEE  
+         static Int_t    NsigmaTh()                    {return fgNsigmaTh;}                        //
+         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(Int_t &padx,Int_t &pady);  //return sector, (padx,pady)->Sector RS
-  static Bool_t  fgIsWireSag;                               //is wire sagitta taken into account
-  static Bool_t  fgIsResolveClusters;                       //performs declustering or not
-  static Int_t   fgHV[6];                                   //HV applied to anod wires
-  static Double_t fgAngleRot;                               //rotation of RICH from up postion (0,0,490)cm
-  static Float_t fSigmaThMap[kNCH][kNpadsX][kNpadsY];       //sigma of the pedestal distributions for all pads
-  static Int_t fgNsigmaTh;                                  //n. of sigmas to cut for zero suppression
-  static Float_t fgSigmaThMean;                             //sigma threshold value
-  static Float_t fgSigmaThSpread;                           //spread of sigma
-  ClassDef(AliRICHParam,4)    //RICH main parameters
+  inline static Int_t    Pad2Sec(TVector pad);              //return sector
+protected:
+         TObjArray *fpChambers;                             //list of chambers    
+  static Bool_t     fgIsWireSag;                            //wire sagitta ON/OFF flag
+  static Bool_t     fgIsResolveClusters;                    //declustering ON/OFF flag
+  static Bool_t     fgIsRadioSrc;                           //radioactive source ON/OFF flag
+  static Int_t      fgHV[6];                                //HV applied to anod wires
+  static Double_t   fgAngleRot;                             //module rotation from up postion (0,0,490)cm
+  static Int_t      fgNsigmaTh;                             //n. of sigmas to cut for zero suppression
+  static Float_t    fgSigmaThMean;                          //sigma threshold value
+  static Float_t    fgSigmaThSpread;                        //spread of sigma
+  ClassDef(AliRICHParam,5)                                  //RICH main parameters class
 };
 //__________________________________________________________________________________________________
 Int_t AliRICHParam::PadNeighbours(Int_t iPadX,Int_t iPadY,Int_t listX[4],Int_t listY[4])
 {
-// Determines all the neighbouring pads for the given one. Returns total amount of these pads.
+// Determines all the neighbouring pads for the given one (iPadX,iPadY). Returns total number of these pads.
 // Dead zones are taken into account.    
+//   1  
+// 2   3
+//   4     
   Int_t nPads=0;
-  if(iPadY!=NpadsY()&&iPadY!=NpadsYsec())                      {listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       
-  if(iPadX!=NpadsXsec()&&iPadX!=2*NpadsXsec()&&iPadX!=NpadsX()){listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       
-  if(iPadY!=1&&iPadY!=NpadsYsec()+1)                           {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}      
-  if(iPadX!=1&&iPadX!=NpadsXsec()+1&&iPadX!=2*NpadsXsec()+1)   {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}
+  if(iPadY!=NpadsY()&&iPadY!=2*NpadsYsec()&&iPadY!=NpadsYsec()){listX[nPads]=iPadX;   listY[nPads]=iPadY+1; nPads++;}       //1
+  if(iPadX!=1&&iPadX!=NpadsXsec()+1)                           {listX[nPads]=iPadX-1; listY[nPads]=iPadY;   nPads++;}       //2
+  if(iPadX!=NpadsXsec()&&iPadX!=NpadsX())                      {listX[nPads]=iPadX+1; listY[nPads]=iPadY;   nPads++;}       //3
+  if(iPadY!=1&&iPadY!=NpadsYsec()+1&&2*NpadsYsec()+1)          {listX[nPads]=iPadX;   listY[nPads]=iPadY-1; nPads++;}       //4
 
   return nPads;
 }//Pad2ClosePads()
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Loc2Sec(TVector2 &x2)
+Int_t AliRICHParam::Loc2Sec(TVector2 &v2)
 {
 // Determines sector containing the given point and trasform this point to the local system of that sector.
-// Returns sector code: 1 2 3
-//                      4 5 6
+// Returns sector code:                       
+//y ^  5 6
+//  |  3 4
+//  |  1 2
+//   -------> x  
+  Double_t x0=0; Double_t x1=SectorSizeX(); Double_t x2=SectorSizeX()+DeadZone(); Double_t x3=PcSizeX();
+  Double_t y0=0; Double_t y1=SectorSizeY(); Double_t y2=SectorSizeY()+DeadZone(); Double_t y3=2*SectorSizeY()+DeadZone(); 
+  Double_t y4=PcSizeY()-SectorSizeY();      Double_t y5=PcSizeY();
+  
   Int_t sector=kBad;  
-  Double_t p1=-0.5*PcSizeX();      Double_t p2=-0.5*SectorSizeX()-DeadZone();  Double_t p3=-0.5*SectorSizeX();
-  Double_t p4= 0.5*SectorSizeX();  Double_t p5= 0.5*SectorSizeX()+DeadZone();  Double_t p6= 0.5*PcSizeX();
-  Double_t x,y;  
-  if     (x2.X()>=p1&&x2.X()<=p2)    {sector=1;x=x2.X()+0.5*PcSizeX();}
-  else if(x2.X()>=p3&&x2.X()<=p4)    {sector=2;x=x2.X()+0.5*SectorSizeX();}
-  else if(x2.X()>=p5&&x2.X()<=p6)    {sector=3;x=x2.X()-0.5*SectorSizeX()-DeadZone();}
-  else                               {return kBad;} //in dead zone or out of chamber
+  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;}
   
-  if     (x2.Y()>=-0.5*PcSizeY() &&x2.Y()<=-0.5*DeadZone())  {y=x2.Y()+0.5*PcSizeY();sector+=3;} //sectors 4,5,6 
-  else if(x2.Y()> -0.5*DeadZone()&&x2.Y()<  0.5*DeadZone())  {return kBad;}                      //in dead zone
-  else if(x2.Y()>= 0.5*DeadZone()&&x2.Y()<= 0.5*PcSizeY())   {y=x2.Y()-0.5*DeadZone();}          //sectors 1,2,3
-  else                                                       {return kBad;}                      //out of chamber    
-  x2.Set(x,y);
+  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;}
+  v2.Set(x,y);
   return sector;
 }//Loc2Sec(Double_t x, Double_t y)
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Loc2Pad(TVector2 x2,Int_t &padx,Int_t &pady)
+TVector AliRICHParam::Loc2Pad(TVector2 x2)
 {
-// Determines pad number (padx,pady) containing the given point x2 defined the chamber RS.
+// Determines pad number TVector(padx,pady) containing the given point x2 defined the chamber RS.
 // Pad count starts in lower left corner from 1,1 to 144,160 in upper right corner of a chamber.
 // Returns sector number of the determined pad.      
+//y ^  5 6
+//  |  3 4
+//  |  1 2
+//   -------> x  
+  TVector pad(2);
   Int_t sector=Loc2Sec(x2);//trasforms x2 to sector reference system
-  if(sector==kBad) {padx=pady=kBad; return sector;}
+  if(sector==kBad) {pad[0]=pad[1]=kBad; return pad;}
   
-  padx=Int_t(x2.X()/PadSizeX())+1; if(padx>NpadsXsec()) padx= NpadsXsec();       
-  if(sector==2||sector==5)   padx+=  NpadsXsec();     // 1 2 3
-  if(sector==3||sector==6)   padx+=2*NpadsXsec();     // 4 5 6
+  pad[0]=Int_t(x2.X()/PadSizeX())+1; if(pad[0]>NpadsXsec()) pad[0]= NpadsXsec();       
+  if(sector==2||sector==4||sector==6)   pad[0]+=  NpadsXsec();     
 
-  pady=Int_t(x2.Y()/PadSizeY())+1; if(pady>NpadsYsec()) pady= NpadsYsec();
-  if(sector<4)               pady+=NpadsYsec();    
-  return sector;
+  pad[1]=Int_t(x2.Y()/PadSizeY())+1; if(pad[1]>NpadsYsec()) pad[1]= NpadsYsec();
+  if(sector==3||sector==4)   pad[1]+=NpadsYsec();    
+  if(sector==5||sector==6)   pad[1]+=2*NpadsYsec();     
+  return pad;
 }
 //__________________________________________________________________________________________________
-Int_t AliRICHParam::Pad2Sec(Int_t &padx, Int_t &pady)
+Int_t AliRICHParam::Pad2Sec(TVector pad)
 {
-// Determines sector containing the given pad (padx,pady) and trasform it to the local RS of that sector.
+// Determines sector containing the given pad.
   Int_t sector=kBad;      
-  if     (padx>=1            &&padx<=NpadsXsec())      {sector=1;}
-  else if(padx> NpadsXsec()  &&padx<=NpadsXsec()*2)    {sector=2;padx-=NpadsXsec();}
-  else if(padx> NpadsXsec()*2&&padx<=NpadsX())         {sector=3;padx-=NpadsXsec()*2;}
-  else                                                 {return 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]);
+    
+  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]);
 
-  if     (pady>=1            &&pady<=NpadsYsec())      {return sector+3;}
-  else if(pady>NpadsYsec()   &&pady<=NpadsY())         {pady-=NpadsYsec();return sector;} 
-  else                                                 {return kBad;}
+  return sector;
 }//Pad2Sec()
 //__________________________________________________________________________________________________
-TVector2 AliRICHParam::Pad2Loc(Int_t padx,Int_t pady)
+TVector2 AliRICHParam::Pad2Loc(TVector pad)
 {
-// Returns position of the center of the given pad (padx,pady) in local RS of the chamber    
-  Int_t sector=Pad2Sec(padx,pady);//shifts to sector RS
-  if(sector==kBad) return TVector2(-101,-101);  
-  Double_t x,y;
-  if(sector<=3)
-    y=0.5*DeadZone()+pady*PadSizeY()-0.5*PadSizeY();   // 1 2 3
-  else{                                                // 4 5 6
-    y=-0.5*PcSizeY()+pady*PadSizeY()-0.5*PadSizeY();
-  }
-  if(sector==1||sector==4)
-    x=-0.5*PcSizeX()+padx*PadSizeX()-0.5*PadSizeX();
-  else if(sector==2||sector==5)
-    x=-0.5*SectorSizeX()+padx*PadSizeX()-0.5*PadSizeX();
+// Returns position of the center of the given pad in local system of the chamber    
+// y ^  5 6
+//   |  3 4        chamber structure
+//   |  1 2
+//    -------> x  
+  Double_t x=kBad,y=kBad;
+  if(pad[0] > 0 && pad[0] <= NpadsXsec())//it's 1 or 3 or 5
+    x=(pad[0]-0.5)*PadSizeX();
+  else if(pad[0] > NpadsXsec() && pad[0] <= NpadsX())//it's 2 or 4 or 6
+    x=(pad[0]-0.5)*PadSizeX()+DeadZone();
   else
-    x= 0.5*SectorSizeX()+DeadZone()+padx*PadSizeX()-0.5*PadSizeX();
+    ::Error("Pad2Loc","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();
+  else if(pad[1] > NpadsYsec() && pad[1] <= 2*NpadsYsec())//it's 3 or 4
+    y=(pad[1]-0.5)*PadSizeY()+DeadZone();
+  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]);
+    
   return TVector2(x,y);
 }
 //__________________________________________________________________________________________________
-Double_t AliRICHParam::GainSag(Double_t y,Int_t sector)
+Double_t AliRICHParam::GainSag(Double_t x,Int_t sector)
 {
 // Returns % of gain variation due to wire sagita.
-// All cureves are parametrized per sector basis, so y must be scaled to the Sector RS.    
-  if(y>0) y-=SectorSizeY()/2; else  y+=SectorSizeY()/2; 
+// All curves are parametrized as per sector basis, so y must be apriory transformed to the Sector RS.    
+  x-=SectorSizeX()/2;
+  if(x>SectorSizeX()) x-=SectorSizeX(); 
   switch(HV(sector)){
-    case 2150: return 9e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0316*TMath::Power(y,2)-3e-4*y+25.367;//%
-    case 2100: return 8e-6*TMath::Power(y,4)+2e-7*TMath::Power(y,3)-0.0283*TMath::Power(y,2)-2e-4*y+23.015;
-    case 2050: return 7e-6*TMath::Power(y,4)+1e-7*TMath::Power(y,3)-0.0254*TMath::Power(y,2)-2e-4*y+20.888;
-    case 2000: return 6e-6*TMath::Power(y,4)+8e-8*TMath::Power(y,3)-0.0227*TMath::Power(y,2)-1e-4*y+18.961;
+    case 2150: return 9e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0316*TMath::Power(x,2)-3e-4*x+25.367;//%
+    case 2100: return 8e-6*TMath::Power(x,4)+2e-7*TMath::Power(x,3)-0.0283*TMath::Power(x,2)-2e-4*x+23.015;
+    case 2050: return 7e-6*TMath::Power(x,4)+1e-7*TMath::Power(x,3)-0.0254*TMath::Power(x,2)-2e-4*x+20.888;
+    case 2000: return 6e-6*TMath::Power(x,4)+8e-8*TMath::Power(x,3)-0.0227*TMath::Power(x,2)-1e-4*x+18.961;
     default:   return 0;
   }
 }
 //__________________________________________________________________________________________________
-Double_t AliRICHParam::QdcSlope(Int_t sec)
-{
-// Returns number of QDC channels per single electron at the unknown yet ???? point for a given sector 
-  switch(sec){
-    case kBad: return 0;
-    default:   return 27;
-  }
-}
-//__________________________________________________________________________________________________
-Double_t AliRICHParam::Gain(TVector2 x2)
-{ 
-//   
-  if(IsWireSag()) 
-    return QdcSlope(Sector(x2))*(1+GainSag(x2.Y(),Sector(x2))/100);
-  else
-    return QdcSlope(Sector(x2));
-}
-//__________________________________________________________________________________________________
 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.
+// 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
-  if(Sector(x2)==kBad) return 0; //hit in the dead zone     
+  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;
   for(Int_t i=1;i<=iNelectrons;i++) qdc+=-Gain(x2)*TMath::Log(gRandom->Rndm());
   return Int_t(qdc);
 }
 //__________________________________________________________________________________________________
-Double_t AliRICHParam::FracQdc(TVector2 x2,Int_t padx,Int_t pady)
+Double_t AliRICHParam::FracQdc(TVector2 x2,TVector pad)
 {
-// Calculates the charge fraction for a given pad (padx,pady) from the given hit point.
-// Mathieson distribution integrated is used.  
-  TVector2 center2=Pad2Loc(padx,pady);//gives center of requested 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(Sector(x2)!=Sector(padx,pady)) return 0;//requested pad does not belong to the sector of given point  
-  else                              return Mathieson(normXmin, normYmin, normXmax, normYmax);
+  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 AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Double_t yMax)
@@ -268,32 +287,40 @@ Double_t AliRICHParam::Mathieson(Double_t xMin,Double_t yMin,Double_t xMax,Doubl
   return 4*kX4*(TMath::ATan(ux2)-TMath::ATan(ux1))*kY4*(TMath::ATan(uy2)-TMath::ATan(uy1));
 }  
 //__________________________________________________________________________________________________
-void AliRICHParam::Loc2Area(TVector2 x2,Int_t &iPadXmin,Int_t &iPadYmin,Int_t &iPadXmax,Int_t &iPadYmax)
+TVector AliRICHParam::Loc2Area(TVector2 x2)
 {
 // Calculates the area of disintegration for a given point. It's assumed here that this points lays on anode wire.
 // Area is a rectangulare set of pads defined by its left-down and right-up coners.
-  Loc2Pad(x2-MathiesonDelta(),iPadXmin,iPadYmin);
-  Loc2Pad(x2+MathiesonDelta(),iPadXmax,iPadYmax);    
+  TVector area(4);
+  TVector pad=Loc2Pad(x2); 
+  area[0]=area[2]=pad[0]; area[1]=area[3]=pad[1];//area is just a pad fired  
+  if(pad[0]!=1           && pad[0]!= NpadsXsec()+1                            ) area[0]--; //left down coner X
+  if(pad[1]!=1           && pad[1]!= NpadsYsec()+1 && pad[1]!= 2*NpadsYsec()+1) area[1]--; //left down coner Y 
+  if(pad[0]!=NpadsXsec() && pad[0]!= NpadsX()                                 ) area[2]++; //right up coner X
+  if(pad[1]!=NpadsYsec() && pad[1]!= 2*NpadsYsec() && pad[1]!= NpadsY()       ) area[3]++; //right up coner Y
+  return area;          
 }
 //__________________________________________________________________________________________________
-Bool_t AliRICHParam::IsOverTh(Int_t c,Int_t x,Int_t y,Double_t q)
+Bool_t AliRICHParam::IsOverTh(Int_t ,TVector ,Double_t q)
 {
-// Calculate the new charge subtracting pedestal and if the current digit is over threshold
-  if (c>0 && x>0 && y>0 && c<kNCH && x<kNpadsX && y<kNpadsY)
-    if(q>NsigmaTh()*fSigmaThMap[c-1][x-1][y-1]) return kTRUE;
-  return kFALSE;
+// 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()));
 }
 //__________________________________________________________________________________________________
-TVector2 AliRICHParam::ShiftToWirePos(TVector2 x2)
+void AliRICHParam::PropogateHelix(TVector3 x0,TVector3 p0,Double_t s,TVector3 *x,TVector3 *p)
 {
-// Calculate the position of the wire nearest to the hit
-  Int_t padx,pady;
-  Loc2Pad(x2,padx,pady);
-  Double_t x;
-  TVector2 center2=Pad2Loc(padx,pady);
-  if(x2.X()>center2.X()) x=center2.X()+0.5*WirePitch();
-  else                   x=center2.X()-0.5*WirePitch();
-  x2.Set(x,x2.Y());
-  return x2;
+// 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 b4b9fd8..9bfebf4 100644 (file)
@@ -22,6 +22,8 @@
 #include <TLorentzVector.h>
 #include <TMath.h>
 
+#include <TGeoManager.h>
+
 ClassImp(AliRICHv0)
 
 void AliRICHv0::StepManager()
@@ -93,11 +95,59 @@ void AliRICHv0::StepManager()
   if(gMC->VolId("CSI ")==gMC->CurrentVolID(copy0)){
     Int_t iChamber;
     gMC->CurrentVolOffID(2,iChamber);
-    TVector3 x3=C(iChamber)->G2L(x4);
-    Info("","loc(%+8.3f,%+8.3f,%8.3f) by G2L",         x3.X(),x3.Y(),x3.Z());  
     TVector2 x2=C(iChamber)->Glob2Loc(x4);
-    Info("","loc(%+8.3f,%+8.3f) by Global2Local",      x2.X(),x2.Y());  
+    Info("","loc(%+8.3f,%+8.3f) by Glob2Loc",      x2.X(),x2.Y());  
   }
   Info("","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 47f0a48..3594ff8 100644 (file)
@@ -9,12 +9,13 @@
 class AliRICHv0 : public AliRICH 
 {    
 public:
-  AliRICHv0():AliRICH()                                              {;}
-  AliRICHv0(const char *name, const char *title):AliRICH(name,title) {;}
-  virtual       ~AliRICHv0()                                         {;}
-  virtual void   Init()                                              {;}
-  virtual Int_t  IsVersion()                                    const{return 0;}
-  virtual void   StepManager();
+  AliRICHv0():AliRICH()                                              {;}                      //default ctor
+  AliRICHv0(const char *name, const char *title):AliRICH(name,title) {;}                      //named ctor
+  virtual       ~AliRICHv0()                                         {;}                      //dtor
+  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 a00b589..8d9efda 100644 (file)
@@ -31,19 +31,34 @@ ClassImp(AliRICHv1)
 //__________________________________________________________________________________________________
 void AliRICHv1::StepManager()
 {
-//Full Step Manager
-
+// Full Step Manager.
+// 3- Ckovs absorbed in Collection electrods
+// 5- Ckovs absorbed in Cathode wires
+// 6- Ckovs absorbed in 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
+  }
+          
 //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(IsLostByFresnel()){ gMC->StopTrack(); return;}        
+    if(gMC->Edep()>0){//CF+CSI+DE
+      if(IsLostByFresnel()){ 
+        if(gMC->TrackPid()==kCerenkov) fCounters(7)++;// 7- Ckovs reflected on CsI
+        gMC->StopTrack();
+        return;
+      }        
       gMC->TrackPosition(cerX4); gMC->CurrentVolOffID(2,iCurrentChamber);
        
       AddHit(iCurrentChamber,gAlice->GetMCApp()->GetCurrentTrackNumber(),cerX4.Vect(),cerX4.Vect());//HIT for PHOTON in conditions CF+CSI+DE
+      fCounters(8)++;//4- Ckovs converted in CsI
       GenerateFeedbacks(iCurrentChamber);
     }//CF+CSI+DE
   }//CF in CSI
@@ -68,6 +83,7 @@ void AliRICHv1::StepManager()
 //__________________________________________________________________________________________________
 Bool_t AliRICHv1::IsLostByFresnel()
 {
+// Calculate probability for the photon to be lost by Fresnel reflection.
   TLorentzVector p4;
   Double_t mom[3],localMom[3];
   gMC->TrackMomentum(p4);   mom[0]=p4(1);   mom[1]=p4(2);   mom[2]=p4(3);
@@ -77,9 +93,84 @@ 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","");
+    if(GetDebug()) Info("IsLostByFresnel","Photon lost");
     return kTRUE;
   }else
     return kFALSE;
 }//IsLostByFresnel()
 //__________________________________________________________________________________________________
+void AliRICHv1::GenerateFeedbacks(Int_t iChamber,Float_t eloss)
+{
+// Generate FeedBack photons for the current particle. To be invoked from StepManager().
+// eloss=0 means photon so only pulse height distribution is to be analysed. This one is done in AliRICHParam::TotQdc()
+  
+  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
+  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 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
+  for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
+    Double_t ranf[2];
+    gMC->GetRandom()->RndmArray(2,ranf);    //Sample direction
+    cthf=ranf[0]*2-1.0;
+    if(cthf<0) continue;
+    sthf = TMath::Sqrt((1 - cthf) * (1 + cthf));
+    phif = ranf[1] * 2 * TMath::Pi();
+    
+    if(Double_t randomNumber=gMC->GetRandom()->Rndm()<=0.57)
+      enfp = 7.5e-9;
+    else if(randomNumber<=0.7)
+      enfp = 6.4e-9;
+    else
+      enfp = 7.9e-9;
+    
+
+    dir[0] = sthf * TMath::Sin(phif);    dir[1] = cthf;    dir[2] = sthf * TMath::Cos(phif);
+    gMC->Gdtom(dir, mom, 2);
+    mom[0]*=enfp;    mom[1]*=enfp;    mom[2]*=enfp;
+    mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
+    
+    // Polarisation
+    e1[0]=      0;    e1[1]=-dir[2];    e1[2]= dir[1];
+    e2[0]=-dir[1];    e2[1]= dir[0];    e2[2]=      0;
+    e3[0]= dir[1];    e3[1]=      0;    e3[2]=-dir[0];
+    
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
+    if (!vmod) for(j=0;j<3;j++) {
+      uswop=e1[j];
+      e1[j]=e3[j];
+      e3[j]=uswop;
+    }
+    vmod=0;
+    for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
+    if (!vmod) for(j=0;j<3;j++) {
+      uswop=e2[j];
+      e2[j]=e3[j];
+      e3[j]=uswop;
+    }
+    
+    vmod=0;  for(j=0;j<3;j++) vmod+=e1[j]*e1[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e1[j]*=vmod;    
+    vmod=0;  for(j=0;j<3;j++) vmod+=e2[j]*e2[j];  vmod=TMath::Sqrt(1/vmod);  for(j=0;j<3;j++) e2[j]*=vmod;
+    
+    phi = gMC->GetRandom()->Rndm()* 2 * TMath::Pi();
+    for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
+    gMC->Gdtom(pol, pol, 2);
+    Int_t outputNtracksStored;    
+    gAlice->GetMCApp()->PushTrack(1,                 //do not transport
+                     gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track 
+                     kFeedback,                      //PID
+                    mom[0],mom[1],mom[2],mom[3],    //track momentum  
+                     x4.X(),x4.Y(),x4.Z(),x4.T(),    //track origin 
+                     pol[0],pol[1],pol[2],           //polarization
+                    kPFeedBackPhoton,
+                     outputNtracksStored,
+                     1.0);    
+  }//feedbacks loop
+  if(GetDebug()) Info("GenerateFeedbacks","Stop.");
+}//GenerateFeedbacks()
index f848de9..41a1ff6 100644 (file)
@@ -15,10 +15,11 @@ public:
   virtual void   Init()                                                              {;}
   virtual Int_t  IsVersion()                                                    const{return 1;}
   
-  virtual void   StepManager();
-          Bool_t IsLostByFresnel(); 
+  virtual void   StepManager();                                         //full slow step manager
+          Bool_t IsLostByFresnel();                                     //checks if the photon lost on Fresnel reflection  
+          void   GenerateFeedbacks(Int_t iChamber,Float_t eloss=0);     //generates feedback photons; eloss=0 for photon
 private:
-  ClassDef(AliRICHv1,1)//RICH full version for simulation
+  ClassDef(AliRICHv1,1)                                                 //RICH full version for simulation
 };
                
 #endif
index 4768f48..87321e3 100644 (file)
@@ -9,7 +9,6 @@
 #pragma link C++ class  AliRICHhit+;
 #pragma link C++ class  AliRICHdigit+;
 #pragma link C++ class  AliRICHcluster+;
-#pragma link C++ class  AliRICHreco+;
 #pragma link C++ class  AliRICHChamber+;
 #pragma link C++ class  AliRICHMap+;
 #pragma link C++ class  AliRICHClusterFinder+;