Code revisited after debugging. Optimization for functions in V3. Calculations of...
authordibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Aug 2008 18:55:23 +0000 (18:55 +0000)
committerdibari <dibari@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 19 Aug 2008 18:55:23 +0000 (18:55 +0000)
HMPID/AliHMPIDParam.cxx
HMPID/AliHMPIDParam.h
HMPID/AliHMPIDPreprocessor.cxx
HMPID/AliHMPIDPreprocessor.h
HMPID/AliHMPIDReconstructor.cxx
HMPID/AliHMPIDTracker.cxx

index 0e03013..60f0468 100644 (file)
@@ -17,6 +17,8 @@
 #include "AliLog.h"         //general
 #include <AliRunLoader.h>   //Stack()
 #include <AliStack.h>       //Stack()
+#include "AliCDBManager.h"  //ctor
+#include "AliCDBEntry.h"    //ctor
 #include <TLatex.h>         //TestTrans()  
 #include <TView.h>          //TestTrans()
 #include <TPolyMarker3D.h>  //TestTrans()
@@ -24,6 +26,8 @@
 #include <TParticle.h>      //Stack()    
 #include <TGeoPhysicalNode.h> //ctor
 #include <TGeoBBox.h>
+#include <TF1.h>                 //ctor
+
 ClassImp(AliHMPIDParam)
 
 
@@ -65,14 +69,33 @@ Int_t AliHMPIDParam::fgSigmas=4;
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 AliHMPIDParam::AliHMPIDParam(Bool_t noGeo=kFALSE):
   TNamed("HmpidParam","default version"),
-  fX(0), fY(0), fRadNmean(0)
+  fX(0), fY(0), fRefIdx(1.28947),fPhotEMean(6.675),fTemp(25)                          //just set a refractive index for C6F14 at ephot=6.675 eV @ T=25 C
 {
 // Here all the intitializition is taken place when AliHMPIDParam::Instance() is invoked for the first time.
 // In particular, matrices to be used for LORS<->MARS trasnformations are initialized from TGeo structure.    
 // Note that TGeoManager should be already initialized from geometry.root file  
 
+  AliCDBManager *pCDB = AliCDBManager::Instance();
+  if(!pCDB) {
+     AliWarning("No Nmean C6F14 from OCDB. Default is taken from ctor.");
+  } else {
+    AliCDBEntry *pNmeanEnt =pCDB->Get("HMPID/Calib/Nmean"); //contains TObjArray of 42 TF1 + 1 EPhotMean
+    if(!pNmeanEnt) {
+      AliWarning("No Nmean C6F14 from OCDB. Default is taken from ctor.");
+    } else {
+      TObjArray *pNmean = (TObjArray*)pNmeanEnt->GetObject();
+      if(pNmean->GetEntries()==43) {                                               //for backward compatibility
+        Double_t tmin,tmax;
+        ((TF1*)pNmean->At(42))->GetRange(tmin,tmax);
+        fPhotEMean = ((TF1*)pNmean->At(42))->Eval(tmin);                          //photon eMean from OCDB
+        AliInfo(Form("EPhotMean = %f eV successfully loaded from OCDB",fPhotEMean));
+      } else {
+        AliWarning("For backward compatibility EPhotMean is taken from ctor.");
+      }
+    }
+  }
 
-  fRadNmean = MeanIdxRad(); //initialization of the running ref. index of freon
+  fRefIdx = MeanIdxRad(); //initialization of the running ref. index of freon
   
   Float_t dead=2.6;// cm of the dead zones between PCs-> See 2CRC2099P1
 
index 533ab69..262c6b7 100644 (file)
@@ -43,8 +43,7 @@ public:
   static Float_t SizeAllX    (                               )     {return fgAllX;                                   }  //all PCs size x, [cm]        
   static Float_t SizeAllY    (                               )     {return fgAllY;                                   }  //all PCs size y, [cm]    
 
-  static Float_t LorsX       (Int_t pc,Int_t padx             )     {return (padx    +0.5)*SizePadX()+fgkMinPcX[pc]; }  //center of the pad x, [cm]
-
+  static Float_t LorsX       (Int_t pc,Int_t padx             )    {return (padx    +0.5)*SizePadX()+fgkMinPcX[pc]; }   //center of the pad x, [cm]
   static Float_t LorsY       (Int_t pc,Int_t pady            )     {return (pady    +0.5)*SizePadY()+fgkMinPcY[pc];  }  //center of the pad y, [cm]
 
   inline static void   Lors2Pad(Float_t x,Float_t y,Int_t &pc,Int_t &px,Int_t &py);                                     //(x,y)->(pc,px,py) 
@@ -58,24 +57,43 @@ public:
 
   static Bool_t  IsOverTh    (Float_t q                      )     {return q >= fgSigmas;                            }  //is digit over threshold?
   
-  Double_t GetRefIdx         (                               )const{return fRadNmean;                                }  //refractive index of freon
   Bool_t  GetInstType        (                               )const{return fgInstanceType;                            }  //return if the instance is from geom or ideal                        
   
   inline static Bool_t IsInDead(Float_t x,Float_t y        );                                                           //is the point in dead area?
   inline static Int_t  InHVSector(           Float_t y     );                                                           //find HV sector
+  static Int_t  Radiator(          Float_t y               )       {if (InHVSector(y)<0) return -1; return InHVSector(y)/2;}
   static Bool_t  IsInside    (Float_t x,Float_t y,Float_t d=0)     {return  x>-d&&y>-d&&x<fgkMaxPcX[kMaxPc]+d&&y<fgkMaxPcY[kMaxPc]+d; } //is point inside chamber boundaries?
 
-            Double_t   MeanIdxRad              ()const {return 1.29204;}   //<--TEMPORAR--> to be removed in future. Mean ref index C6F14
-            Double_t   MeanIdxWin              ()const {return 1.57819;}   //<--TEMPORAR--> to be removed in future. Mean ref index quartz
-            Float_t    DistCut                 ()const {return 1.0;}       //<--TEMPORAR--> to be removed in future. Cut for MIP-TRACK residual 
-            Float_t    QCut                    ()const {return 100;}       //<--TEMPORAR--> to be removed in future. Separation PHOTON-MIP charge 
-            Float_t    MultCut                 ()const {return 200;}       //<--TEMPORAR--> to be removed in future. Multiplicity cut to activate WEIGHT procedure 
+  //For optical properties
+  static Double_t   EPhotMin()                       {return 5.5;}           //
+  static Double_t   EPhotMax()                       {return 8.5;}           //Photon energy range,[eV]
+  static Double_t NIdxRad(Double_t eV,Double_t temp) {return TMath::Sqrt(1+0.554*(1239.84/eV)*(1239.84/eV)/((1239.84/eV)*(1239.84/eV)-5769)-0.0005*(temp-20));}
+  static Double_t NIdxWin(Double_t eV)               {return TMath::Sqrt(1+46.411/(10.666*10.666-eV*eV)+228.71/(18.125*18.125-eV*eV));}  
+  static Double_t NMgF2Idx(Double_t eV)              {return 1.7744 - 2.866e-3*(1239.842609/eV) + 5.5564e-6*(1239.842609/eV)*(1239.842609/eV);}          // MgF2 idx of trasparency system
+  static Double_t NIdxGap(Double_t eV)               {return 1+0.12489e-6/(2.62e-4 - eV*eV/1239.84/1239.84);}
+  static Double_t LAbsRad(Double_t eV)               {return (eV<7.8)*(GausPar(eV,3.20491e16,-0.00917890,0.742402)+GausPar(eV,3035.37,4.81171,0.626309))+(eV>=7.8)*0.0001;}
+  static Double_t LAbsWin(Double_t eV)               {return (eV<8.2)*(818.8638-301.0436*eV+36.89642*eV*eV-1.507555*eV*eV*eV)+(eV>=8.2)*0.0001;}//fit from DiMauro data 28.10.03
+  static Double_t LAbsGap(Double_t eV)               {return (eV<7.75)*6512.399+(eV>=7.75)*3.90743e-2/(-1.655279e-1+6.307392e-2*eV-8.011441e-3*eV*eV+3.392126e-4*eV*eV*eV);}
+  static Double_t QEffCSI(Double_t eV)               {return (eV>6.07267)*0.344811*(1-exp(-1.29730*(eV-6.07267)));}//fit from DiMauro data 28.10.03
+  static Double_t GausPar(Double_t x,Double_t a1,Double_t a2,Double_t a3) {return a1*TMath::Exp(-0.5*((x-a2)/a3)*((x-a2)/a3));}
+  inline static Double_t FindTemp(Double_t tLow,Double_t tUp,Double_t y);    //find the temperature of the C6F14 in a given point with coord. y (in x is uniform)
+  
+  
+  Double_t   GetEPhotMean            ()const {return fPhotEMean;} 
+  Double_t   GetRefIdx               ()const {return fRefIdx;}                       //running refractive index
+  
+  Double_t   MeanIdxRad              ()const {return NIdxRad(fPhotEMean,fTemp);}
+  Double_t   MeanIdxWin              ()const {return NIdxWin(fPhotEMean);}
+  //
+  Float_t    DistCut                 ()const {return 1.0;}       //<--TEMPORAR--> to be removed in future. Cut for MIP-TRACK residual 
+  Float_t    QCut                    ()const {return 100;}       //<--TEMPORAR--> to be removed in future. Separation PHOTON-MIP charge 
+  Float_t    MultCut                 ()const {return 200;}       //<--TEMPORAR--> to be removed in future. Multiplicity cut to activate WEIGHT procedure 
 
-            Double_t   RadThick                ()const {return 1.5;}       //<--TEMPORAR--> to be removed in future. Radiator thickness
-            Double_t   WinThick                ()const {return 0.5;}       //<--TEMPORAR--> to be removed in future. Window thickness
-            Double_t   GapThick                ()const {return 8.0;}       //<--TEMPORAR--> to be removed in future. Proximity gap thickness
-            Double_t   WinIdx                  ()const {return 1.5787;}    //<--TEMPORAR--> to be removed in future. Mean refractive index of WIN material (SiO2) 
-            Double_t   GapIdx                  ()const {return 1.0005;}    //<--TEMPORAR--> to be removed in future. Mean refractive index of GAP material (CH4)
+  Double_t   RadThick                ()const {return 1.5;}       //<--TEMPORAR--> to be removed in future. Radiator thickness
+  Double_t   WinThick                ()const {return 0.5;}       //<--TEMPORAR--> to be removed in future. Window thickness
+  Double_t   GapThick                ()const {return 8.0;}       //<--TEMPORAR--> to be removed in future. Proximity gap thickness
+  Double_t   WinIdx                  ()const {return 1.5787;}    //<--TEMPORAR--> to be removed in future. Mean refractive index of WIN material (SiO2) 
+  Double_t   GapIdx                  ()const {return 1.0005;}    //<--TEMPORAR--> to be removed in future. Mean refractive index of GAP material (CH4)
 
   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  
@@ -90,9 +108,13 @@ public:
                                                                                            ph=TMath::ATan2(l[1],l[0]);}    
   TVector3 Norm        (Int_t c                                             )const{Double_t n[3]; Norm(c,n); return TVector3(n);               }//norm 
   void     Norm        (Int_t c,Double_t *n                                 )const{Double_t l[3]={0,0,1};fM[c]->LocalToMasterVect(l,n);        }//norm
-  void     Point       (Int_t c,Double_t *p,Int_t plane                     )const{Lors2Mars(c,0,0,p,plane);}      //point of given chamber plane
+  void     Point       (Int_t c,Double_t *p,Int_t plane                     )const{Lors2Mars(c,0,0,p,plane);}         //point of given chamber plane
 
-  void     SetRefIdx      (Double_t refRadIdx                                  ) {fRadNmean = refRadIdx;}             //set refractive index of freon
+  void     SetTemp        (Double_t temp                                       ) {fTemp = temp;}                      //set actual temperature of the C6F14
+  void     SetEPhotMean   (Double_t ePhotMean                                  ) {fPhotEMean = ePhotMean;}            //set mean photon energy
+  
+  void     SetRefIdx      (Double_t refRadIdx                                  ) {fRefIdx = refRadIdx;}               //set running refractive index
+  
   void     SetSigmas      (Int_t sigmas                                        ) {fgSigmas = sigmas;}                 //set sigma cut    
   void     SetInstanceType(Bool_t inst                                         ) {fgInstanceType = inst;}             //kTRUE if from geomatry kFALSE if from ideal geometry
   //For PID
@@ -143,14 +165,15 @@ protected:
 
   TGeoHMatrix *fM[7];                 //pointers to matrices defining HMPID chambers rotations-translations
   Float_t fX;                         //x shift of LORS with respect to rotated MARS 
-  Float_t fY;                         //y shift of LORS with respect to rotated MARS   
-  Double_t fRadNmean;                 //C6F14 mean index as a running parameter
-  
+  Float_t fY;                         //y shift of LORS with respect to rotated MARS
+  Double_t fRefIdx;                   //running refractive index of C6F14
+  Double_t fPhotEMean;                //mean energy of photon
+  Double_t fTemp;                     //actual temparature of C6F14  
 private:
   AliHMPIDParam(const AliHMPIDParam& r);              //dummy copy constructor
   AliHMPIDParam &operator=(const AliHMPIDParam& r);   //dummy assignment operator
       
-  ClassDef(AliHMPIDParam,0)           //HMPID main parameters class
+  ClassDef(AliHMPIDParam,1)           //HMPID main parameters class
 };
 
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -214,4 +237,15 @@ Int_t AliHMPIDParam::InHVSector(Float_t y)
    return hvsec;
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDParam::FindTemp(Double_t tLow,Double_t tHigh,Double_t y)
+{
+//  Model for gradient in temperature
+  
+//  Double_t gradT = (t2-t1)/SizePcY();  // linear gradient
+//  return gradT*y+t1;
+  Double_t halfPadSize = 0.5*SizePadY();
+  Double_t gradT = (TMath::Log(SizePcY()) - TMath::Log(halfPadSize))/(TMath::Log(tHigh)-TMath::Log(tLow));
+  return tLow*TMath::Power(y/halfPadSize,1./gradT);  
+}
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 #endif
index 4596192..22bc559 100644 (file)
@@ -56,7 +56,7 @@ UInt_t AliHMPIDPreprocessor::Process(TMap* pMap)
        Log("HMPID - Pedestal processing successful!!");                    return kFALSE;  // ok for pedestals
     }
   } else if ( runType=="STANDALONE" || runType=="PHYSICS"){
-    if (!ProcDcs(pMap)){
+  if (!ProcDcs(pMap)){
        Log("HMPID - ERROR - DCS processing failed!!");                     return kTRUE;   // error in DCS processing
     } else {
        Log("HMPID - DCS processing successful!!");                         return kFALSE;  // ok for DCS
@@ -68,7 +68,9 @@ UInt_t AliHMPIDPreprocessor::Process(TMap* pMap)
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
 {
-// Process: 1. inlet and outlet C6F14 temperature, stores TObjArray of 21 TF1, where TF1 is Nmean=f(t), one per radiator
+// Process: 1 (old). inlet and outlet C6F14 temperature, stores TObjArray of 21 TF1, where TF1 is Nmean=f(t), one per radiator
+// Process: 1. inlet and outlet C6F14 temperature, stores TObjArray of 42 TF1, where TF1 are Tin and Tout per radiator
+//             + one function for the mean energy photon (in total 43).
 //          2. CH4 pressure and HV                 stores TObjArray of 7 TF1 where TF1 is thr=f(t), one per chamber
 // Arguments: pDcsMap - map of structure "alias name" - TObjArray of AliDCSValue
 // Assume that: HV is the same during the run for a given chamber, different chambers might have different HV
@@ -77,22 +79,22 @@ Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
 
   Bool_t stDcsStore=kFALSE;
 
-  TF2 idx("RidxC4F14","sqrt(1+0.554*(1239.84/x)^2/((1239.84/x)^2-5769)-0.0005*(y-20))",5.5 ,8.5 ,0  ,50);  //N=f(Ephot,T) [eV,grad C] DiMauro mail
-  
 // Qthr=f(HV,P) [V,mBar]  logA0=k*HV+b is taken from p. 64 TDR plot 2.59 for PC32 
 //                           A0=f(P) is taken from DiMauro mail
 // Qthr is estimated as 3*A0
+
   TF2 thr("RthrCH4"  ,"3*10^(3.01e-3*x-4.72)+170745848*exp(-y*0.0162012)"             ,2000,3000,900,1200); 
   
-  TObjArray arTmean(21);       arTmean.SetOwner(kTRUE);     //21 Tmean=f(time) one per radiator
-  TObjArray arPress(7);        arPress.SetOwner(kTRUE);     //7  Press=f(time) one per chamber
-  TObjArray arNmean(21);       arNmean.SetOwner(kTRUE);     //21 Nmean=f(time) one per radiator
+  TObjArray arNmean(43);       arNmean.SetOwner(kTRUE);     //42 Tin and Tout one per radiator + 1 for ePhotMean
   TObjArray arQthre(42);       arQthre.SetOwner(kTRUE);     //42 Qthre=f(time) one per sector
   
   AliDCSValue *pVal; Int_t cnt=0;
   
   Double_t xP,yP;
 
+  TF1 **pTin  = new TF1*[21];
+  TF1 **pTout = new TF1*[21];
+
 // evaluate Environment Pressure
   
   TObjArray *pPenv=(TObjArray*)pMap->GetValue("HMP_DET/HMP_ENV/HMP_ENV_PENV.actual.value");
@@ -151,11 +153,13 @@ Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
       arQthre.AddAt(new TF1(Form("HMP_QthreC%iS%i",iCh,iSec),
           Form("3*10^(3.01e-3*HV%i_%i - 4.72)+170745848*exp(-(P%i+Penv)*0.0162012)",iCh,iSec,iCh),fStartTime,fEndTime),6*iCh+iSec);
     }
-    
 // evaluate Temperatures: in and out of the radiators    
-    
     // T in
     for(Int_t iRad=0;iRad<3;iRad++){
+      
+      pTin[3*iCh+iRad]  = new TF1(Form("Tin%i%i" ,iCh,iRad),"[0]+[1]*x",fStartTime,fEndTime);
+      pTout[3*iCh+iRad] = new TF1(Form("Tout%i%i",iCh,iRad),"[0]+[1]*x",fStartTime,fEndTime);          
+      
       TObjArray *pT1=(TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_MP%i/HMP_MP%i_LIQ_LOOP.actual.sensors.Rad%iIn_Temp",iCh,iCh,iRad));  
       Log(Form(" Temperatures for module %i inside data  ---> %3i entries",iCh,pT1->GetEntries()));
       if(pT1->GetEntries()) {
@@ -163,9 +167,10 @@ Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
         TGraph *pGrT1=new TGraph; cnt=0;  while((pVal=(AliDCSValue*)nextT1())) pGrT1->SetPoint(cnt++,pVal->GetTimeStamp(),pVal->GetFloat()); //T inlet
         if(cnt==1) { 
           pGrT1->GetPoint(0,xP,yP);
-          new TF1(Form("Tin%i%i",iCh,iRad),Form("%f",yP),fStartTime,fEndTime);
+          pTin[3*iCh+iRad]->SetParameter(0,yP);
+          pTin[3*iCh+iRad]->SetParameter(1,0);
         } else {
-          pGrT1->Fit(new TF1(Form("Tin%i%i",iCh,iRad),"[0]+[1]*x+[2]*sin([3]*x)",fStartTime,fEndTime),"Q");
+          pGrT1->Fit(pTin[3*iCh+iRad],"Q");
         }
         delete pGrT1;
       } else {AliWarning(" No Data Points from HMP_MP0-6_LIQ_LOOP.actual.sensors.Rad0-2In_Temp!");return kFALSE;}
@@ -177,20 +182,25 @@ Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
         TGraph *pGrT2=new TGraph; cnt=0;  while((pVal=(AliDCSValue*)nextT2())) pGrT2->SetPoint(cnt++,pVal->GetTimeStamp(),pVal->GetFloat()); //T outlet 
         if(cnt==1) { 
           pGrT2->GetPoint(0,xP,yP);
-          new TF1(Form("Tou%i%i",iCh,iRad),Form("%f",yP),fStartTime,fEndTime);
+          pTout[3*iCh+iRad]->SetParameter(0,yP);
+          pTout[3*iCh+iRad]->SetParameter(1,0);
         } else {
-          pGrT2->Fit(new TF1(Form("Tou%i%i",iCh,iRad),"[0]+[1]*x+[2]*sin([3]*x)",fStartTime,fEndTime),"Q");
+          pGrT2->Fit(pTout[3*iCh+iRad],"Q");
         }
         delete pGrT2;
       } else {AliWarning(" No Data Points from HMP_MP0-6_LIQ_LOOP.actual.sensors.Rad0-2Out_Temp!");return kFALSE;}
        
 // evaluate Mean Refractive Index
       
-      arNmean.AddAt(new TF1(Form("HMP_Nmean%i-%i",iCh,iRad),"1.292",fStartTime,fEndTime),3*iCh+iRad); //Nmean=f(t)
+      arNmean.AddAt(pTin[3*iCh+iRad] ,6*iCh+2*iRad  ); //Tin =f(t)
+      arNmean.AddAt(pTout[3*iCh+iRad],6*iCh+2*iRad+1); //Tout=f(t)
       
     }//radiators loop
   }//chambers loop
   
+  Double_t eMean = ProcTrans(pMap);
+  arNmean.AddAt(new TF1("HMP_PhotEmean",Form("%f",eMean),fStartTime,fEndTime),42); //Photon energy mean
+    
   AliCDBMetaData metaData; 
   metaData.SetBeamPeriod(0); 
   metaData.SetResponsible("AliHMPIDPreprocessor"); 
@@ -201,6 +211,13 @@ Bool_t AliHMPIDPreprocessor::ProcDcs(TMap* pMap)
   if(!stDcsStore) {
     Log("HMPID - failure to store DCS data results in OCDB");    
   }
+  
+//  arNmean.Delete();
+//  arQthre.Delete();
+
+//  for(Int_t i=0;i<21;i++) delete pTin[i]; delete []pTin;
+//  for(Int_t i=0;i<21;i++) delete pTout[i]; delete []pTout;
+  
   return stDcsStore;
 }//Process()
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -270,10 +287,107 @@ Bool_t AliHMPIDPreprocessor::ProcPed()
   
 }//ProcPed()  
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t ProcTrans()
+Double_t AliHMPIDPreprocessor::ProcTrans(TMap* pMap)
 {
-// Process transparency monitoring data and calculates Emean  
-  Double_t eMean=6.67786;     //mean energy of photon defined  by transperancy window
+  //  Process transparency monitoring data and calculates Emean  
+
+  
+  Double_t sEnergProb=0, sProb=0;
+
+  Double_t tRefCR5 = 19. ;                                      // mean temperature of CR5 where the system is in place
+
+  Double_t eMean = 0;
+      
+  AliDCSValue *pVal;
+
+  for(Int_t i=0; i<30; i++){
+
+    // evaluate wavelenght 
+    TObjArray *pWaveLenght  = (TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_INFR/HMP_INFR_TRANPLANT/HMP_INFR_TRANPLANT_MEASURE.measure[%i].waveLenght",i));
+    TIter NextWl(pWaveLenght); pVal=(AliDCSValue*)NextWl();
+    Double_t lambda = pVal->GetFloat();
+
+    Double_t photEn = 1239.842609/lambda;     // 1239.842609 from nm to eV
+    
+    if(photEn<AliHMPIDParam::EPhotMin() || photEn>AliHMPIDParam::EPhotMax()) continue;
+    
+    // evaluate phototube current for argon reference
+    TObjArray *pArgonRef  = (TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_INFR/HMP_INFR_TRANPLANT/HMP_INFR_TRANPLANT_MEASURE.measure[%i].argonReference",i));
+    TIter NextArRef(pArgonRef); pVal=(AliDCSValue*)NextArRef();
+    Double_t aRefArgon = pVal->GetFloat();
+
+    // evaluate phototube current for argon cell
+    TObjArray *pArgonCell = (TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_INFR/HMP_INFR_TRANPLANT/HMP_INFR_TRANPLANT_MEASURE.measure[%i].argonCell",i));
+    TIter NextArCell(pArgonCell); pVal=(AliDCSValue*)NextArCell();
+    Double_t aCellArgon = pVal->GetFloat();
+
+    //evaluate phototube current for freon reference
+    TObjArray *pFreonRef  = (TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_INFR/HMP_INFR_TRANPLANT/HMP_INFR_TRANPLANT_MEASURE.measure[%i].c6f14Reference",i));
+    TIter NextFrRef(pFreonRef); pVal=(AliDCSValue*)NextFrRef();
+    Double_t aRefFreon = pVal->GetFloat();
+
+    //evaluate phototube current for freon cell
+    TObjArray *pFreonCell = (TObjArray*)pMap->GetValue(Form("HMP_DET/HMP_INFR/HMP_INFR_TRANPLANT/HMP_INFR_TRANPLANT_MEASURE.measure[%i].c6f14Cell",i));
+    TIter NextFrCell(pFreonCell); pVal=(AliDCSValue*)NextFrCell();
+    Double_t aCellFreon = pVal->GetFloat();
+   //evaluate correction factor to calculate trasparency (Ref. NIMA 486 (2002) 590-609)
+    
+    Double_t aN1 = AliHMPIDParam::NIdxRad(photEn,tRefCR5);
+    Double_t aN2 = AliHMPIDParam::NMgF2Idx(photEn);
+    Double_t aN3 = 1;                              // Argon Idx
+
+    Double_t aR1               = ((aN1 - aN2)*(aN1 - aN2))/((aN1 + aN2)*(aN1 + aN2));
+    Double_t aR2               = ((aN2 - aN3)*(aN2 - aN3))/((aN2 + aN3)*(aN2 + aN3));
+    Double_t aT1               = (1 - aR1);
+    Double_t aT2               = (1 - aR2);
+    Double_t aCorrFactor       = (aT1*aT1)/(aT2*aT2);
+
+    // evaluate 15 mm of thickness C6F14 Trans
+    Double_t aTransRad;
+    
+    if(aRefFreon*aRefArgon>0) {
+      aTransRad  = TMath::Power((aCellFreon/aRefFreon)/(aCellArgon/aRefArgon)*aCorrFactor,1.5);
+    } else {
+      return DefaultEMean();
+    }
+
+    // evaluate 0.5 mm of thickness SiO2 Trans
+    Double_t aTransSiO2 = TMath::Exp(-0.5/AliHMPIDParam::LAbsWin(photEn));
+
+    // evaluate 80 cm of thickness Gap (low density CH4) transparency
+    Double_t aTransGap  = TMath::Exp(-80./AliHMPIDParam::LAbsGap(photEn));
+
+    // evaluate CsI quantum efficiency
+    Double_t aCsIQE            = AliHMPIDParam::QEffCSI(photEn);
+
+    // evaluate total convolution of all material optical properties
+    Double_t aTotConvolution   = aTransRad*aTransSiO2*aTransGap*aCsIQE;
+
+    sEnergProb+=aTotConvolution*photEn;  
+  
+    sProb+=aTotConvolution;  
+}
+
+  if(sProb>0) {
+    eMean = sEnergProb/sProb;
+  } else {
+    return DefaultEMean();
+  }
+
+  Log(Form(" Mean energy photon calculated ---> %f eV ",eMean));
+
+  if(eMean<AliHMPIDParam::EPhotMin() || eMean>AliHMPIDParam::EPhotMax()) return DefaultEMean();
+  
   return eMean;
+
 }   
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+Double_t AliHMPIDPreprocessor::DefaultEMean()
+{
+    Double_t eMean = 6.675;      //just set a refractive index for C6F14 at ephot=6.675 eV @ T=25 C
+    AliWarning(Form("Mean energy for photons out of range [%f,%f] in Preprocessor. Default value Eph=%f eV taken.",AliHMPIDParam::EPhotMin(),
+                                                                                                                   AliHMPIDParam::EPhotMax(),
+                                                                                                                   eMean));
+    return eMean;
+}
index af6a7af..63be3d3 100644 (file)
@@ -21,10 +21,12 @@ public:
            }
   virtual ~AliHMPIDPreprocessor(                             )                                 {}
 protected:
-  virtual void   Initialize(Int_t run, UInt_t startTime, UInt_t endTime); //
-  virtual UInt_t Process   (TMap* pDcsMap                              ); //process everthing
-          Bool_t ProcDcs   (TMap* pDcsMap                              ); //process DCS data points
-          Bool_t ProcPed   (                                           ); //process pedestal files
+  virtual void     Initialize(Int_t run, UInt_t startTime, UInt_t endTime); //
+  virtual UInt_t   Process   (TMap* pDcsMap                              ); //process everthing
+          Bool_t   ProcDcs   (TMap* pDcsMap                              ); //process DCS data points
+          Bool_t   ProcPed   (                                           ); //process pedestal files
+          Double_t ProcTrans (TMap *pDcsMap                              );
+          Double_t DefaultEMean();                                          //set a default value in ePhotMean                             
   ClassDef(AliHMPIDPreprocessor, 0);
 };
 
index e1fe219..5d50600 100644 (file)
@@ -63,28 +63,6 @@ AliHMPIDReconstructor::AliHMPIDReconstructor():AliReconstructor(),fUserCut(0),fD
     }
   } 
 
-  AliCDBEntry *pQthreEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Qthre"); //contains TObjArray of 7 TF1
-  if(!pQthreEnt) AliFatal("No Qthre available");
-  TObjArray *pQthre = (TObjArray*)pQthreEnt->GetObject();
-  for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++) {
-    if(pQthre->GetEntriesFast()==AliHMPIDParam::kMaxCh+1) {        // for backward compatibility...
-      TF1 *pfQthre = (TF1*)pQthre->At(iCh); 
-      Double_t tMin,tMax;
-      pfQthre->GetRange(tMin,tMax);
-      Double_t qthre=pfQthre->Eval(tMin);
-      AliDebug(1,Form("Qthre successfully loaded for chamber %i -> %f ",iCh,qthre));
-    } else {
-      for(Int_t isec=0;isec<=5;isec++) {
-       TF1 *pfQthre = (TF1*)pQthre->At(6*iCh+isec); 
-       Double_t tMin,tMax;
-       pfQthre->GetRange(tMin,tMax);
-       Double_t qthre=pfQthre->Eval(tMin);
-       AliDebug(1,Form("Qthre successfully loaded for chamber %i  sector %i -> %f ",iCh,isec,qthre));
-      }
-    }
-  }
-
-     
   AliCDBEntry *pDaqSigEnt =AliCDBManager::Instance()->Get("HMPID/Calib/DaqSig");  //contains TObjArray of TObjArray 14 TMatrixF sigmas values for pads 
   if(!pDaqSigEnt) AliFatal("No pedestals from DAQ!");
   fDaqSig = (TObjArray*)pDaqSigEnt->GetObject();
index 8b7af80..02a3009 100644 (file)
@@ -80,7 +80,7 @@ Int_t AliHMPIDTracker::PropagateBack(AliESDEvent *pEsd)
 // Interface pure virtual in AliTracker. Invoked from AliReconstruction::RunTracking() after invocation of AliTracker::LoadClusters() once per event
 // Agruments: pEsd - pointer to ESD
 //   Returns: error code    
-  AliCDBEntry *pNmeanEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Nmean"); //contains TObjArray of 21 TF1
+  AliCDBEntry *pNmeanEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Nmean"); //contains TObjArray of 42 TF1 + 1 EPhotMean
   AliCDBEntry *pQthreEnt =AliCDBManager::Instance()->Get("HMPID/Calib/Qthre"); //contains TObjArray of 42 (7ch * 6sec) TF1
   if(!pNmeanEnt) AliFatal("No Nmean C6F14 ");
   if(!pQthreEnt) AliFatal("No Qthre");
@@ -105,14 +105,23 @@ Int_t AliHMPIDTracker::Recon(AliESDEvent *pEsd,TObjArray *pClus,TObjArray *pNmea
       continue;                                                                         
     }
     pTrk->SetHMPIDtrk(xRa,yRa,theta,phi);                                                        //store initial infos
-    Double_t nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp());                       //C6F14 Nmean for this chamber
+    Double_t nmean;
+    if(pNmean->GetEntries()==21) {                                                               //for backward compatibility
+       nmean=((TF1*)pNmean->At(3*cham))->Eval(pEsd->GetTimeStamp());                             //C6F14 Nmean for this chamber
+     } else {
+       Int_t iRad     = AliHMPIDParam::Radiator(yRa);                                            //evaluate the radiator involved
+       Double_t tLow  = ((TF1*)pNmean->At(6*cham+2*iRad  ))->Eval(pEsd->GetTimeStamp());         //C6F14 low  temp for this chamber
+       Double_t tHigh = ((TF1*)pNmean->At(6*cham+2*iRad+1))->Eval(pEsd->GetTimeStamp());         //C6F14 high temp for this chamber
+       Double_t tExp  = AliHMPIDParam::FindTemp(tLow,tHigh,yRa);                                 //estimated temp for that chamber at that y
+       nmean = AliHMPIDParam::NIdxRad(AliHMPIDParam::Instance()->GetEPhotMean(),tExp);           //mean ref idx @ a given temp
+     }
     Double_t qthre = 0; 
     if(pQthre->GetEntriesFast()==AliHMPIDParam::kMaxCh+1)                                        // just for backward compatibility
       qthre=((TF1*)pQthre->At(cham))->Eval(pEsd->GetTimeStamp());                                //
     else {                                                                                       // in the past just 1 qthre
       Int_t hvsec = AliHMPIDParam::InHVSector(yPc);                                              //  per chamber
       if (hvsec>=0)
-       qthre=((TF1*)pQthre->At(6*cham+hvsec))->Eval(pEsd->GetTimeStamp());                        //
+       qthre=((TF1*)pQthre->At(6*cham+hvsec))->Eval(pEsd->GetTimeStamp());                      //
     }                                                                                            //
     recon.SetImpPC(xPc,yPc);                                                                     //store track impact to PC
     recon.CkovAngle(pTrk,(TClonesArray *)pClus->At(cham),nmean,qthre);                           //search for Cerenkov angle of this track