Removing static data members (TF1 objects) from AliRICHParam to make it working on...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 6 Apr 2006 12:58:52 +0000 (12:58 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 6 Apr 2006 12:58:52 +0000 (12:58 +0000)
RICH/AliRICHHelix.h
RICH/AliRICHParam.cxx
RICH/AliRICHParam.h
RICH/AliRICHTracker.cxx

index a769f1c..ff8a17c 100644 (file)
@@ -82,23 +82,24 @@ Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam)
 // Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection
 // On exit fPosRad contain position of intersection in radiator LORS (cm)    
 //         fPosPc  contains the same for photocathode  
-  for(Int_t iChamN=1;iChamN<=AliRICHParam::kNch;iChamN++){//chamber loop
-    if(Intersection(pParam->Center(iChamN,AliRICHParam::kRad),pParam->Norm(iChamN))){//there is intersection with radiator plane
-      fPosRad=pParam->Mars2Lors(iChamN,fX,AliRICHParam::kRad);//position on radiator plane
+  for(Int_t iCh=1;iCh<=AliRICHParam::kNch;iCh++){//chambers loop
+    TVector3 norm =pParam->Lors2MarsVec(iCh,TVector3(0,0,1));
+    if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kRad),norm)){//there is intersection with radiator plane
+      fPosRad=pParam->Mars2Lors(iCh,fX,AliRICHParam::kRad);//position on radiator plane
       if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone)
         
-        if(Intersection(pParam->Center(iChamN,AliRICHParam::kPc),pParam->Norm(iChamN))){//there is intersection with photocathode
-          fPosPc=pParam->Mars2Lors(iChamN,fX,AliRICHParam::kPc);//position on radiator plane
+        if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kPc),norm)){//there is intersection with photocathode
+          fPosPc=pParam->Mars2Lors(iCh,fX,AliRICHParam::kPc);//position on radiator plane
           if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone)
             
-            fPloc=pParam->Mars2LorsVec(iChamN,fP);//trasform p to local system
-            return iChamN;
+            fPloc=pParam->Mars2LorsVec(iCh,fP);//trasform p to local system
+            return iCh;
           }//if inside PC
         }//if intersects PC
         
       }//if inside radiator
     }//if for radiator       
-  }//chamber loop
+  }//chambers loop
   return 0;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
index 62efb83..863d913 100644 (file)
 
 ClassImp(AliRICHParam)
 AliRICHParam * AliRICHParam::fgInstance             =0x0;     //singleton pointer               
-Bool_t         AliRICHParam::fgIsWireSag            =kTRUE;   //take ware sagita into account?
-Bool_t         AliRICHParam::fgIsResolveClusters    =kTRUE;   //do cluster resolving?
-Bool_t         AliRICHParam::fgIsFeedback           =kTRUE;   //generate feedback photons?
-Bool_t         AliRICHParam::fgIsTestBeam           =kFALSE;  //special test beam configuration
 
-Int_t    AliRICHParam::fgHV[kNsectors]        ={2050,2050,2050,2050,2050,2050};
-Int_t    AliRICHParam::fgNsigmaTh             =4;
-Float_t  AliRICHParam::fgSigmaThMean          =1.132; //QDC 
-Float_t  AliRICHParam::fgSigmaThSpread        =0.035; //     
-Double_t AliRICHParam::fgErrChrom[4][330];                       //
-Double_t AliRICHParam::fgErrGeom[4][330];                        //
-Double_t AliRICHParam::fgErrLoc[4][330];                         //Chromatic, Geometric and Localization array to parametrize SigmaCerenkov
 Double_t AliRICHParam::fgMass[5]              ={0.00051,0.10566,0.13957,0.49360,0.93828};  
 
 
-Double_t AliRICHParam::fEckovMin=5.5e-9; //GeV
-Double_t AliRICHParam::fEckovMax=8.5e-9; //GeV
-TF1 AliRICHParam::fgAbsC6F14("RabsC4F14","6512.39*(x<=7.75e-9)+(x>7.75e-9)*0.039/(-0.166+0.063e9*x-8.01e7*x^2+3.39e5*x^3)" ,fEckovMin,fEckovMax); 
-TF1 AliRICHParam::fgAbsSiO2 ("RabsSiO2" ,"333"                                                                             ,fEckovMin,fEckovMax); 
-TF1 AliRICHParam::fgAbsCH4  ("RabsCH4"  ,"6512.39*(x<=7.75e-9)+(x>7.75e-9)*0.039/(-0.166+0.063e9*x-8.01e7*x^2+3.39e5*x^3)" ,fEckovMin,fEckovMax);
-TF1 AliRICHParam::fgAbsAir  ("RabsAir"  ,"500"                                                                             ,fEckovMin,fEckovMax);  //len ???
-
-TF1 AliRICHParam::fgIdxAir  ("RidxAir"  ,"1+1e-8*(8342.13 + 2406030/(130-(1.23984e-9/x)^2)+15597/(38.9-(1.23984e-9/x)^2))" ,fEckovMin,fEckovMax);  //???
-TF1 AliRICHParam::fgIdxSiO2 ("RidxSiO2" ,"sqrt(1+46.411/(10.666*10.666-x*x*1e18)+228.71/(18.125*18.125-x*x*1e18))"         ,fEckovMin,fEckovMax);  //TDR p.35
-TF1 AliRICHParam::fgIdxCH4  ("RidxCH4"  ,"1+0.12489e-6/(2.62e-4 - (1239.84e-9/x)^-2)"                                      ,fEckovMin,fEckovMax);  //Olav preprint
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliRICHParam::AliRICHParam():TNamed("RichParam","default version") 
 {
@@ -116,7 +95,6 @@ void AliRICHParam::Print(Option_t* opt) const
 {
 // print some usefull (hopefully) info on some internal guts of RICH parametrisation 
   Printf("Pads in chamber (%3i,%3i) in sector (%2i,%2i) pad size (%4.2f,%4.2f)",NpadsX(),NpadsY(),NpadsXsec(),NpadsYsec(),PadSizeX(),PadSizeY());
-  Printf("Resolve clusters %i sagita %i",IsResolveClusters(),IsWireSag()); 
   
   for(Int_t i=0;i<kNchambers;i++) fMatrix[i]->Print(opt);
 }//Print()
@@ -284,142 +262,6 @@ void AliRICHParam::DrawSectors()
   TPolyLine *sec6 = new TPolyLine(5,xRight,yUp);      sec6->SetLineColor(21);  sec6->Draw();
 }//DrawSectors()
 //__________________________________________________________________________________________________
-void AliRICHParam::ReadErrFiles()
-{
-// Read the three files corresponding to Chrom,Geom and Loc They are parameters of a polynomial of 6th order...  ????????? go to CDB?
-// Arguments: none
-//   Returns: none    
-  
-  static Bool_t count = kFALSE;
-  
-  Float_t c0,c1,c2,c3,c;
-  Float_t g0,g1,g2,g3,g;
-  Float_t l0,l1,l2,l3,l;
-  
-  FILE *pChromErr, *pGeomErr, *pLocErr;  
-
-  if(!count) {
-     AliInfoGeneral("ReadErrFiles","reading RICH error parameters...");
-     pChromErr = fopen(Form("%s/RICH/RICHConfig/SigmaChromErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
-     pGeomErr  = fopen(Form("%s/RICH/RICHConfig/SigmaGeomErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
-     pLocErr   = fopen(Form("%s/RICH/RICHConfig/SigmaLocErr.txt",gSystem->Getenv("ALICE_ROOT")),"r");
-     if(!pChromErr||!pGeomErr||!pLocErr) {AliErrorGeneral("ReadErrFiles"," RICH ERROR READING Parameter FILES: can't open files!!! ");return;}
-     for(Int_t i=0;i<330;i++) {
-       fscanf(pChromErr,"%f%f%f%f%f\n",&c0,&c1,&c2,&c3,&c);
-       fscanf(pGeomErr,"%f%f%f%f%f\n",&g0,&g1,&g2,&g3,&g);
-       fscanf(pLocErr,"%f%f%f%f%f\n",&l0,&l1,&l2,&l3,&l);
-       fgErrChrom[0][i] = c0;
-       fgErrChrom[1][i] = c1;
-       fgErrChrom[2][i] = c2;
-       fgErrChrom[3][i] = c3;  
-       fgErrGeom[0][i] = g0;
-       fgErrGeom[1][i] = g1;
-       fgErrGeom[2][i] = g2;
-       fgErrGeom[3][i] = g3;   
-       fgErrLoc[0][i] = l0;
-       fgErrLoc[1][i] = l1;
-       fgErrLoc[2][i] = l2;
-       fgErrLoc[3][i] = l3;    
-     }
-     AliInfoGeneral("ReadErrFiles","DONE successfully!");
-     fclose(pChromErr);
-     fclose(pGeomErr);
-     fclose(pLocErr);
-  }
-  count = kTRUE;
-}//ReadErrFiles()
-//__________________________________________________________________________________________________
-TVector3 AliRICHParam::SigmaSinglePhoton(Int_t partID, Double_t mom, Double_t theta, Double_t phi)
-
-{
-// Find sigma for single photon. It returns the thrree different errors. If you want
-// to have the error---> TVector3.Mag()
-// partID = 0,1,2,3,4 ---> e,mu,pi,k,p in agreement with AliPID
-  TVector3 v(-999,-999,-999);
-  Double_t pmom;
-
-  ReadErrFiles();
-  Double_t mass = fgMass[partID];
-  Double_t massRef = fgMass[4]; // all the files are calculated for protons...so mass ref is proton mass
-  pmom = mom*massRef/mass; // normalized momentum respect to proton...
-  if(pmom>PmodMax()) pmom = PmodMax();
-  Double_t oneOverRefIndex = 1/IdxC6F14(EckovMean());
-  Double_t pmin = mass*oneOverRefIndex/TMath::Sqrt(1-oneOverRefIndex*oneOverRefIndex);
-  if(pmom<pmin) return v;
-  Double_t Theta = theta*TMath::RadToDeg();
-  if(phi<0) phi+=TMath::TwoPi();
-  Double_t Phi = phi*TMath::RadToDeg();
-  v.SetX(Interpolate(fgErrChrom,pmom,Theta,Phi));
-  v.SetY(Interpolate(fgErrGeom,pmom,Theta,Phi));
-  v.SetZ(Interpolate(fgErrLoc,pmom,Theta,Phi));
-//  v*=1.5; // take into account bigger errors due to multiplicity...to change in future
-
-  return v;
-}//SigmaSinglePhoton
-//__________________________________________________________________________________________________
-TVector3 AliRICHParam::SigmaSinglePhoton(Double_t thetaCer, Double_t theta, Double_t phi)
-
-{
-// Find sigma for single photon. It returns the thrree different errors. If you want
-// to have the error---> TVector3.Mag()
-// partID = 0,1,2,3,4 ---> e,mu,pi,k,p in agreement with AliPID
-  TVector3 v(-999,-999,-999);
-  Double_t pmom;
-
-  ReadErrFiles();
-  Double_t massRef = fgMass[4]; // all the files are calculated for protons...so mass ref is proton mass
-  Double_t beta=1./(IdxC6F14(EckovMean())*TMath::Cos(thetaCer));
-  if(beta>=1) {
-    pmom=6.5; // above physical limi the error is calculated at the saturation...
-  } else {
-    Double_t gamma=1./TMath::Sqrt(1-beta*beta);
-    pmom = beta*gamma*massRef; // normalized momentum respect to proton...
-  }
-  if(pmom>PmodMax()) pmom = PmodMax();
-  Double_t oneOverRefIndex = 1/IdxC6F14(EckovMean());
-  Double_t pmin = massRef*oneOverRefIndex/TMath::Sqrt(1-oneOverRefIndex*oneOverRefIndex);
-  if(pmom<pmin) return v;
-  Double_t Theta = theta*TMath::RadToDeg();
-  if(phi<0) phi+=TMath::TwoPi();
-  Double_t Phi = phi*TMath::RadToDeg();
-  v.SetX(Interpolate(fgErrChrom,pmom,Theta,Phi));
-  v.SetY(Interpolate(fgErrGeom,pmom,Theta,Phi));
-  v.SetZ(Interpolate(fgErrLoc,pmom,Theta,Phi));
-//  v*=1.5; // take into account bigger errors due to multiplicity...to change in future
-
-  return v;
-}//SigmaSinglePhoton
-//__________________________________________________________________________________________________
-Double_t AliRICHParam::Interpolate(Double_t par[4][330], Double_t x, Double_t y, Double_t phi)
-  
-{
-  static Double_t amin = 1.15; static Double_t astep  = 0.2;
-  static Double_t bmin = 0; static Double_t bstep  = 1;
-  
-  Double_t Phi = (phi - 180)/300.;
-  
-  Double_t Sigma[30][11];
-  
-  for(Int_t j=0;j<11;j++) { for(Int_t i=0;i<30;i++) {
-    Sigma[i][j] = par[0][j+11*i] + par[1][j+11*i]*Phi*Phi + par[2][j+11*i]*TMath::Power(Phi,4) + par[3][j+11*i]*TMath::Power(Phi,6);
-    }
-  }
-
-  Int_t i=0;Int_t j=0;
-  
-  i = (Int_t)((x-amin)/astep);
-  j = (Int_t)((y-bmin)/bstep);
-  Double_t ai = amin+i*astep;
-  Double_t ai1 = ai+astep;
-  Double_t bj = bmin+j*bstep;
-  Double_t bj1 = bj+bstep;
-  Double_t t = (x-ai)/(ai1-ai);
-  Double_t gj = (1-t)*Sigma[i][j]+t*Sigma[i+1][j];
-  Double_t gj1 = (1-t)*Sigma[i][j+1]+t*Sigma[i+1][j+1];
-  Double_t u = (y-bj)/(bj1-bj);
-  return (1-u)*gj+u*gj1;
-}//Interpolate
-//__________________________________________________________________________________________________
 TVector3 AliRICHParam::ForwardTracing(TVector3 entranceTrackPoint, TVector3 vectorTrack, Double_t thetaC, Double_t phiC)
 {
 // Trace a single Ckov photon from a given emission point up to photocathode taking into account ref indexes of materials it travereses
index 8db972e..0c3278c 100644 (file)
@@ -39,16 +39,12 @@ public:
   static void     DrawSectors();
 //flags staff         
   static inline AliRICHParam* Instance();                                //pointer to AliRICHParam singleton
-  static                void     SetWireSag(Bool_t a)               {fgIsWireSag=a;}
-  static                Bool_t   IsWireSag()                        {return fgIsWireSag;}
-  static                   void     SetResolveClusters(Bool_t a)    {fgIsResolveClusters=a;}
-  static                   Bool_t   IsResolveClusters()             {return fgIsResolveClusters;}
   static        Int_t      Stack(Int_t evt=-1,Int_t tid=-1);              //Print stack info for event and tid
   static        Int_t      StackCount(Int_t pid,Int_t evt);               //Counts stack particles of given sort in given event  
   static inline Double_t   ErrLoc                  (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
   static inline Double_t   ErrGeom                 (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
   static inline Double_t   ErrCrom                 (Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
-  static inline TVector3   SigmaSinglePhotonFormula(Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
+  static inline Double_t   SigmaSinglePhotonFormula(Double_t thetaC,Double_t phiC,Double_t thetaT,Double_t phiT,Double_t beta);
 //Geometrical properties  
   static        Int_t      NpadsX      ()   {return kNpadsX;}                           //number of pads along X in chamber
   static        Int_t      NpadsY      ()   {return kNpadsY;}                           //number of pads along Y in chamber
@@ -76,11 +72,12 @@ public:
   
 //trasformation methodes
          inline TVector3   Lors2Mars     (Int_t c,Double_t x,Double_t y,Int_t p=kPc); //LORS->MARS transform of point [cm] for chamber c and plane p
+         inline TVector3   Lors2MarsVec  (Int_t c,const TVector3 &p                ); //LORS->MARS transform of vector for chamber c
          inline TVector2   Mars2Lors     (Int_t c,const TVector3 &x    ,Int_t p=kPc); //MARS->LORS transform of point [cm] for chamber c and plane p    
+         inline TVector3   Mars2LorsVec  (Int_t c,const TVector3 &p                ); //MARS->LORS transform of vector for chamber c
   
   static inline TVector3   Lors2MarsOld  (Int_t c,Double_t x,Double_t y,Int_t p); //LORS->MARS transform of position (cm) for chamber c and plane p
   static inline TVector2   Mars2LorsOld  (Int_t c,const TVector3 &x,Int_t p    ); //MARS->LORS transform of position (cm) for chamber c and plane p    
-  static inline TVector3   Mars2LorsVec  (Int_t c,const TVector3 &p            ); //MARS->LORS transform of vector for chamber c
   static inline TVector3   Center        (Int_t c,Int_t p                      ); //Center of plane p of chamber c in MARS (cm)
   static inline TVector3   Norm          (Int_t c                              ); //Norm vector to the chamber c in MARS (cm)
   static inline TGeoMatrix*Matrix        (Int_t iCh, Int_t iPlane              ); //TGeoMatrix for the given chamber plain
@@ -103,16 +100,13 @@ public:
   static        Bool_t     IsAccepted    (const TVector2 &x2             ){return ( x2.X()>=0 && x2.X()<=PcSizeX() && x2.Y()>=0 && x2.Y()<=PcSizeY() ) ? kTRUE:kFALSE;}
 //optical properties methodes  
   static        Float_t    EckovMean     (                               ){return 6.766e-9;}                          //mean Ckov energy according to the total trasmission curve
-  static        Float_t    EckovMin      (                               ){return fEckovMin;}                         //min photon energy [GeV] defined in optical curves
-  static        Float_t    EckovMax      (                               ){return fEckovMax;}                         //min photon energy [GeV] defined in optical curves
-  
-  static        Float_t    AbsCH4        (Float_t gev                    );                                           //AbsLen [cm]=f(Eckov) [GeV] for CH4 used as amp gas
-  static        Float_t    AbsAir        (Float_t gev                    ){return fgAbsAir.Eval(gev);}                //AbsLen [cm]=f(Eckov) [GeV] for air
+  static        Float_t    EckovMin      (                               ){return 5.5e-9;}                            //min photon energy [GeV] defined in optical curves
+  static        Float_t    EckovMax      (                               ){return 8.5e-9;}                            //min photon energy [GeV] defined in optical curves
   
                 Float_t    IdxC6F14      (Float_t gev                    ){return fIdxC6F14->Eval(gev,fIdxC6F14->GetUniqueID());}  //n=f(Eckov) [GeV] for C6H14 used as radiator
-  static        Float_t    IdxSiO2       (Float_t gev                    ){return fgIdxSiO2 .Eval(gev);}              //n=f(Eckov) [GeV] for SiO2 used as window TDR p.35
-  static        Float_t    IdxCH4        (Float_t gev                    ){return fgIdxCH4  .Eval(gev);}              //n=f(Eckov) [GeV] for CF4 
-  static        Float_t    IdxAir        (Float_t gev                    ){return fgIdxAir  .Eval(gev);}              //n=f(Eckov) [GeV] for air
+  static        Float_t    IdxSiO2       (Float_t gev                    ){return TMath::Sqrt(1+46.411/(10.666*10.666-gev*gev*1e18)+228.71/(18.125*18.125-gev*gev*1e18));} //n=f(Eckov) [GeV] for SiO2 used as window TDR p.35
+  static        Float_t    IdxCH4        (Float_t gev                    ){return 1+0.12489e-6/(2.62e-4 - TMath::Power(1239.84e-9/gev,-2));}              //n=f(Eckov) [GeV] for CF4 
+  static        Float_t    AbsCH4        (Float_t gev                    );                                                                          //abs len=f(Eckov) [GeV] for CF4 
   
                 void       CdbRead   (Int_t run,Int_t version           );                                           //read all calibration information for requested run
   
@@ -124,15 +118,14 @@ public:
   static Int_t    QthMIP()                   {return 100;}
   static Double_t DmatchMIP()                {return 1.;}
   static Double_t PmodMax()                  {return 6.5;}
-  static Int_t    HV(Int_t sector)           {if (sector>=1 && sector <=6) return fgHV[sector-1];  else return -1;} //high voltage for this sector
-  static void     SetHV(Int_t sector,Int_t hv){fgHV[sector-1]=hv;}  
+  static Int_t    HV(Int_t sector)           {if (sector>=1 && sector <=6) return 2050;  else return -1;} //high voltage for this sector
 //charge response methodes  
   inline static Double_t Mathieson(Double_t x1,Double_t x2,Double_t y1,Double_t y2);               //Mathienson integral over given limits
   
   inline static Double_t GainSag(Double_t x,Int_t sector);                                         //gain variations in %
-         static Double_t Gain(const TVector2 &x2){//gives chamber gain in terms of QDC channels for given point in local ref system
-                          if(fgIsWireSag) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);
-                          else            return QdcSlope(Loc2Sec(x2));}
+         static Double_t Gain(const TVector2 &x2,Bool_t isSag=kTRUE){//gives chamber gain in terms of QDC channels for given point in local ref system
+                          if(isSag) return QdcSlope(Loc2Sec(x2))*(1+GainSag(x2.X(),Loc2Sec(x2))/100);
+                          else      return QdcSlope(Loc2Sec(x2));}
   inline static Double_t FracQdc(const TVector2 &x2,const TVector &pad);                           //charge fraction to pad from hit
   inline static Int_t    TotQdc(TVector2 x2,Double_t e    );                                       //total charge for Eloss (GeV) 0 for photons
          static Double_t QdcSlope(Int_t sec){switch(sec){case -1: return 0;  default:   return 33;}} //weight of electon in QDC channels
@@ -141,23 +134,19 @@ public:
   static        Double_t IonPot        (                                              ){return 26.0e-9;}                            //for CH4 in GeV taken from ????
   static inline Int_t    QdcTot        (Int_t iPad,Double_t e                         );                                            //total QDC generated by Eloss or Etot [GeV]
   static inline Double_t QdcSag        (Int_t iPad                                    );                                            //mean QDC variation due to sagita [0,1]
-  static        Double_t QdcEle        (Int_t iPad                                    ){return fgIsWireSag?33*(1+QdcSag(iPad)):33;} //mean QDC per electron
+  static        Double_t QdcEle        (Int_t iPad,Bool_t isSag=kTRUE                 ){return isSag?33*(1+QdcSag(iPad)):33;}       //mean QDC per electron
   static inline Int_t    Hit2SDigs     (Int_t iPad,  Double_t e,TClonesArray* pSDigLst);                                            //hit->sdigits, returns Qtot
   static inline Int_t    Hit2SDigs     (TVector2 hit,Double_t e,TClonesArray* pSDigLst);                                            //hit->sdigits, returns Qtot, old style
   static        void     TestHit2SDigs (Double_t x,Double_t y,Double_t e,Bool_t isNew=kFALSE);                                      //test hit->sdigits
   
   inline static 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
+         static Int_t    NsigmaTh()                    {return 4;}                        //
+         static Float_t  SigmaThMean()                 {return 1.132;}                    //QDC electronic noise mean
+         static Float_t  SigmaThSpread()               {return 0.035;}                    //QDC electronic noise width
                 
          static Double_t CogCorr(Double_t x) {return 3.31267e-2*TMath::Sin(2*TMath::Pi()/PadSizeX()*x) //correction of cluster CoG due to sinoidal
                                                     -2.66575e-3*TMath::Sin(4*TMath::Pi()/PadSizeX()*x)
                                                     +2.80553e-3*TMath::Sin(6*TMath::Pi()/PadSizeX()*x)+0.0070;}
-         static void     ReadErrFiles();                                                                  //Read Err file parameters
-         TVector3 SigmaSinglePhoton(Int_t Npart, Double_t mom, Double_t theta, Double_t phi);      //Find Sigma for single photon from momentum and particle id
-         TVector3 SigmaSinglePhoton(Double_t thetaCer, Double_t theta, Double_t phi);              //Fing sigma for single photon from thetacer
-  static Double_t Interpolate(Double_t par[4][330],Double_t x, Double_t y, Double_t phi);          //Find the error value from interpolation
        
          TVector3 ForwardTracing(TVector3 entranceTrackPoint,TVector3 vectorTrack, Double_t thetaC, Double_t phiC); //it traces foward a photon from Emission Point to PC
   static TVector3 PlaneIntersect(const TVector3 &lineDir,const TVector3 &linePoint,const TVector3 &planeNorm,const TVector3 &planePoint); //intersection between line and plane
@@ -171,39 +160,13 @@ public:
   static void     TestTrans();                                                           //test the transform group of methodes
 
   static Double_t fgMass[5];                                // mass array
-  static Bool_t     fgIsTestBeam;                           //test beam geometry instead of normal RICH flag
   enum EPlaneId {kCenter,kPc,kRad,kAnod,kNch=7};            //4 planes in chamber and total number of chambers
 protected:
          AliRICHParam();             //default ctor is protected to enforce it to be singleton
   static AliRICHParam *fgInstance;   //static pointer  to instance of AliRICHParam singleton
-//optical curves
-  static Double_t fEckovMin;         //min Eckov
-  static Double_t fEckovMax;         //max Eckov
-      
-         TF2*  fIdxC6F14;            //n=f(Ephot,T) [GeV] for radiator freon   C6F14
-  static TF1   fgIdxSiO2;            //n=f(Ephot) [GeV] for window   quartz  SiO2 
-  static TF1   fgIdxCH4;             //n=f(Ephot) [GeV] for MWPC     amp gas CF4  
-  static TF1   fgIdxAir;             //n=f(Ephot) [GeV] for          air
-  
-  static TF1   fgAbsC6F14;           //abs len curve for radiator freon   C6F14, cm   versus GeV
-  static TF1   fgAbsSiO2;            //abs len curve for window   quartz  SiO2 , cm   versus GeV  
-  static TF1   fgAbsCH4;             //abs len curve for MWPC     methane CF4  , cm   versus GeV
-  static TF1   fgAbsAir;             //abs len curve for air, cm versus GeV
-  
-  static Bool_t     fgIsWireSag;            //wire sagitta ON/OFF flag
-  static Bool_t     fgIsResolveClusters;    //declustering ON/OFF flag
-  static Bool_t     fgIsFeedback;           //generate feedback photon?
-
-  static Int_t      fgHV[6];                //HV applied to anod wires
-  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
-
-  static Double_t     fgErrChrom[4][330];                       //
-  static Double_t     fgErrGeom[4][330];                        //
-  static Double_t     fgErrLoc[4][330];                         //Chromatic, Geometric and Localization array to parametrize SigmaCerenkov
-         TGeoHMatrix *fMatrix[kNchambers];                      //poiners to matrices defining RICH chambers rotations-translations
-  ClassDef(AliRICHParam,0)                                      //RICH main parameters class
+  TF2         *fIdxC6F14;            //n=f(Ephot,T) [GeV] for radiator freon   C6F14
+  TGeoHMatrix *fMatrix[kNchambers];  //poiners to matrices defining RICH chambers rotations-translations
+  ClassDef(AliRICHParam,0)           //RICH main parameters class
 };
 
 AliRICHParam* AliRICHParam::Instance()
@@ -519,6 +482,13 @@ TVector3  AliRICHParam::Lors2Mars(Int_t iChId,Double_t x,Double_t y,Int_t iPlnId
   return TVector3(mars);
 }    
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+TVector3  AliRICHParam::Lors2MarsVec(Int_t iCh,const TVector3 &p)
+{
+  Double_t mars[3], lors[3]; p.GetXYZ(lors);  
+  fMatrix[iCh-1]->LocalToMasterVect(lors,mars);
+  return TVector3(mars);
+}    
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector2  AliRICHParam::Mars2Lors(Int_t iChId,const TVector3 &x,Int_t iPlnId)
 {
 // Trasform from MARS to LORS
@@ -538,6 +508,13 @@ TVector2  AliRICHParam::Mars2Lors(Int_t iChId,const TVector3 &x,Int_t iPlnId)
   return TVector2(lors[0]+0.5*PcSizeX(),lors[1]+0.5*PcSizeY());
 }    
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+TVector3  AliRICHParam::Mars2LorsVec(Int_t iCh,const TVector3 &p)
+{
+  Double_t mars[3], lors[3]; p.GetXYZ(mars);  
+  fMatrix[iCh-1]->MasterToLocalVect(mars,lors);
+  return TVector3(lors);
+}    
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector3  AliRICHParam::Lors2MarsOld(Int_t iChId,Double_t x,Double_t y,Int_t iPlnId)
 {
 // Trasform from LORS to MARS
@@ -556,13 +533,6 @@ TVector2  AliRICHParam::Mars2LorsOld(Int_t iChamN,const TVector3 &x,Int_t iPlane
   return TVector2(lors[0]+0.5*PcSizeX(),lors[1]+0.5*PcSizeY());
 }    
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-TVector3  AliRICHParam::Mars2LorsVec(Int_t iChamN,const TVector3 &x)
-{
-  TGeoMatrix *pMatrix=Matrix(iChamN,kPc);
-  Double_t mars[3]={x.X(),x.Y(),x.Z()}  , lors[3];  pMatrix->MasterToLocalVect(mars,lors);  delete pMatrix;
-  return TVector3(lors);
-}    
-//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 TVector3  AliRICHParam::Center(Int_t iChamN,Int_t iPlaneN)
 {
   TGeoMatrix *pMatrix=Matrix(iChamN,iPlaneN);
@@ -634,15 +604,15 @@ Int_t AliRICHParam::Hit2SDigs(TVector2 hitX2,Double_t e,TClonesArray *pSDigLst)
   return iQtot;
 }//Hit2SDigs() for TVector2
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-TVector3 AliRICHParam::SigmaSinglePhotonFormula(Double_t thetaCer, Double_t phiCer, Double_t theta, Double_t phi, Double_t beta)
+Double_t AliRICHParam::SigmaSinglePhotonFormula(Double_t thetaCer, Double_t phiCer, Double_t theta, Double_t phi, Double_t beta)
 {
   TVector3 v(-999,-999,-999);
 
-  v.SetX(AliRICHParam::ErrLoc(thetaCer,phiCer,theta,phi,beta));
-  v.SetY(AliRICHParam::ErrGeom(thetaCer,phiCer,theta,phi,beta));
-  v.SetZ(AliRICHParam::ErrCrom(thetaCer,phiCer,theta,phi,beta));
+  v.SetX(ErrLoc(thetaCer,phiCer,theta,phi,beta));
+  v.SetY(ErrGeom(thetaCer,phiCer,theta,phi,beta));
+  v.SetZ(ErrCrom(thetaCer,phiCer,theta,phi,beta));
 
-  return v;
+  return v.Mag2();
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Double_t AliRICHParam::ErrLoc(Double_t thetaC, Double_t phiC, Double_t Ptheta, Double_t Pphi, Double_t beta)
@@ -670,7 +640,7 @@ Double_t RefC6F14m = 1.29337;
 
   return  TMath::Sqrt(0.2*0.2*dtdxc*dtdxc + 0.25*0.25*dtdyc*dtdyc);
 }
-
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Double_t AliRICHParam::ErrCrom(Double_t thetaC, Double_t phiC, Double_t Ptheta, Double_t Pphi, Double_t beta)
 {
   Double_t dphi = phiC - Pphi;
index d4cdd5b..6ed9db4 100644 (file)
@@ -60,6 +60,7 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
   AliDebug(1,Form("Start with %i tracks",iNtracks));
   AliRICH *pRich=((AliRICH*)gAlice->GetDetector("RICH"));  
   AliRICHRecon recon;
+  AliPID pid; // needed to retrive all the PID info
    
   for(Int_t iTrk=0;iTrk<iNtracks;iTrk++){//ESD tracks loop
     AliESDtrack *pTrack = pESD->GetTrack(iTrk);// get next reconstructed track    
@@ -95,31 +96,33 @@ Int_t AliRICHTracker::PropagateBack(AliESD *pESD)
 //here comes PID calculations    
     if(pTrack->GetRICHsignal()>0) {
       AliDebug(1,Form("Start to assign the probabilities"));
-      Double_t sigmaPID[AliPID::kSPECIES];
-      Double_t richPID[AliPID::kSPECIES];
-      for (Int_t iPart=0;iPart<AliPID::kSPECIES;iPart++) {
-        sigmaPID[iPart] = 0;
-        fErrPar[iPart] = 0;
-        for(Int_t iphot=0;iphot<pRich->Clus(iChamber)->GetEntries();iphot++) {
-          recon.SetPhotonIndex(iphot);
+      Double_t sigmaPID[AliPID::kSPECIES];  Double_t richPID[AliPID::kSPECIES];
+      for (Int_t iPid=0;iPid<AliPID::kSPECIES;iPid++){//PID loop
+        sigmaPID[iPid] = 0; fErrPar[iPid] = 0;
+        for(Int_t iCkov=0;iCkov<pRich->Clus(iChamber)->GetEntries();iCkov++){//Ckov candidates loop ????????????? somehting fomr AliRICHRecon must be
+          recon.SetPhotonIndex(iCkov);
           if(recon.GetPhotonFlag() == 2) {
-            Double_t theta_g=recon.GetTrackTheta();
-            Double_t phi_g=(recon.GetPhiPoint()-recon.GetTrackPhi());
-            Double_t sigma2 = AliRICHParam::Instance()->SigmaSinglePhoton(iPart,pTrack->GetP(),theta_g,phi_g).Mag2(); 
-            if(sigma2>0) sigmaPID[iPart] += 1/sigma2;
+            Double_t thetaCkov=0.6; //??????????????????
+            Double_t   phiCkov=0.6; //??????????????????
+            Double_t  thetaTrk=recon.GetTrackTheta();
+            Double_t    phiTrk=(recon.GetPhiPoint()-recon.GetTrackPhi());
+            Double_t     beta =pTrack->GetP()/TMath::Sqrt((pTrack->GetP()*pTrack->GetP()+pid.ParticleMass(iPid)*pid.ParticleMass(iPid)));
+            Double_t   sigma2 = AliRICHParam::SigmaSinglePhotonFormula(thetaCkov,phiCkov,thetaTrk,phiTrk,beta);
+            if(sigma2>0) sigmaPID[iPid] += 1/sigma2;
           }
-        }
-        if (sigmaPID[iPart]>0)
-          sigmaPID[iPart] *= (Double_t)(iMipId-recon.GetPhotBKG())/(Double_t)(iMipId); // n total phots, m are background...the sigma are scaled..
-          if(sigmaPID[iPart]>0) sigmaPID[iPart] = 1/TMath::Sqrt(sigmaPID[iPart])*0.001; // sigma from parametrization are in mrad...
-          else                  sigmaPID[iPart] = 0;
-          fErrPar[iPart]=sigmaPID[iPart];
-        AliDebug(1,Form("sigma for %s is %f rad",AliPID::ParticleName(iPart),sigmaPID[iPart]));
+        }//Ckov candidates loop
+        
+        if (sigmaPID[iPid]>0)
+          sigmaPID[iPid] *= (Double_t)(iMipId-recon.GetPhotBKG())/(Double_t)(iMipId); // n total phots, m are background...the sigma are scaled..
+          if(sigmaPID[iPid]>0) sigmaPID[iPid] = 1/TMath::Sqrt(sigmaPID[iPid])*0.001; // sigma from parametrization are in mrad...
+          else                  sigmaPID[iPid] = 0;
+          fErrPar[iPid]=sigmaPID[iPid];
+        AliDebug(1,Form("sigma for %s is %f rad",AliPID::ParticleName(iPid),sigmaPID[iPid]));
       }
       CalcProb(pTrack->GetRICHsignal(),pTrack->GetP(),sigmaPID,richPID);
       pTrack->SetRICHpid(richPID);         
       AliDebug(1,Form("PROBABILITIES ---> %f - %f - %f - %f - %f",richPID[0],richPID[1],richPID[2],richPID[3],richPID[4]));
-    }    
+    }//if(pTrack->GetRICHsignal())    
   }//ESD tracks loop
   AliDebug(1,"Stop pattern recognition");
   return 0; // error code: 0=no error;