Protection instead fatal in case of absent galice.root or Kinematics.root

[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx

diff --git a/PHOS/AliPHOSClusterizerv1.cxx b/PHOS/AliPHOSClusterizerv1.cxx
index 5e4fdcd725f20dfbbf0960d24e41d77d153ec2c7..5f5d221c391dbda69a0743a75f729512dd7d4622 100644 (file)
--- a/PHOS/AliPHOSClusterizerv1.cxx
+++ b/PHOS/AliPHOSClusterizerv1.cxx
@@ -198,11 +198,10 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), 
   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
-  fW0CPV(0),                  fEmcTimeGate(0)
+  fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 {
   // default ctor (to be used mainly by Streamer)
   
-  InitParameters() ; 
   fDefaultInit = kTRUE ; 
 }
 
@@ -214,11 +213,10 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
   fNumberOfEmcClusters(0),    fNumberOfCpvClusters(0),
   fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), 
   fEmcLocMaxCut(0),           fW0(0),                   fCpvLocMaxCut(0),
-  fW0CPV(0),                  fEmcTimeGate(0)
+  fW0CPV(0),                  fEmcTimeGate(0),          fEcoreRadius(0)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   
-  InitParameters() ;
   Init() ;
   fDefaultInit = kFALSE ; 
 }
@@ -264,6 +262,8 @@ void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
     AliInfo(Form("took %f seconds for Clusterizing\n",
                 gBenchmark->GetCpuTime("PHOSClusterizer"))); 
   }
+  fEMCRecPoints->Delete();
+  fCPVRecPoints->Delete();
 }
 
 //____________________________________________________________________________
@@ -381,24 +381,24 @@ void AliPHOSClusterizerv1::InitParameters()
   fNumberOfCpvClusters     = 0 ; 
   fNumberOfEmcClusters     = 0 ; 
 
-  const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
-  if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
+  const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
+  if(!recoParam) AliFatal("Reconstruction parameters are not set!");
 
-  const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv(); 
-  if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
+  recoParam->Print();
 
-  fCpvClusteringThreshold  = parCpv->GetClusteringThreshold();
-  fEmcClusteringThreshold  = parEmc->GetClusteringThreshold();
+  fEmcClusteringThreshold  = recoParam->GetEMCClusteringThreshold();
+  fCpvClusteringThreshold  = recoParam->GetCPVClusteringThreshold();
   
-  fEmcLocMaxCut            = parEmc->GetLocalMaxCut();
-  fCpvLocMaxCut            = parCpv->GetLocalMaxCut();
+  fEmcLocMaxCut            = recoParam->GetEMCLocalMaxCut();
+  fCpvLocMaxCut            = recoParam->GetCPVLocalMaxCut();
 
-  fW0                      = parEmc->GetLogWeight();
-  fW0CPV                   = parCpv->GetLogWeight();
+  fW0                      = recoParam->GetEMCLogWeight();
+  fW0CPV                   = recoParam->GetCPVLogWeight();
 
   fEmcTimeGate             = 1.e-6 ; 
+  fEcoreRadius             = recoParam->GetEMCEcoreRadius();
   
-  fToUnfold                = parEmc->ToUnfold() ;
+  fToUnfold                = recoParam->EMCToUnfold() ;
     
   fWrite                   = kTRUE ;
 }
@@ -486,7 +486,8 @@ void AliPHOSClusterizerv1::WriteRecPoints()
   Int_t index ;
   //Evaluate position, dispersion and other RecPoint properties..
   Int_t nEmc = fEMCRecPoints->GetEntriesFast();
-  Float_t emcMinE= AliPHOSReconstructor::GetRecoParamEmc()->GetMinE(); //Minimal digit energy
+  Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
+  TVector3 fakeVtx(0.,0.,0.) ;
   for(index = 0; index < nEmc; index++){
     AliPHOSEmcRecPoint * rp =
       dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
@@ -498,8 +499,8 @@ void AliPHOSClusterizerv1::WriteRecPoints()
     }
 
     // No vertex is available now, calculate corrections in PID
-    rp->EvalAll(fW0,fDigitsArr) ;
-    TVector3 fakeVtx(0.,0.,0.) ;
+    rp->EvalAll(fDigitsArr) ;
+    rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
     rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
     rp->EvalLocal2TrackingCSTransform();
   }
@@ -516,7 +517,8 @@ void AliPHOSClusterizerv1::WriteRecPoints()
   //Now the same for CPV
   for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
     AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
-    rp->EvalAll(fW0CPV,fDigitsArr) ;
+    rp->EvalAll(fDigitsArr) ;
+    rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
     rp->EvalLocal2TrackingCSTransform();
   }
   fCPVRecPoints->Sort() ;
@@ -541,14 +543,20 @@ void AliPHOSClusterizerv1::MakeClusters()
   fNumberOfEmcClusters     = 0 ;
 
   //Mark all digits as unused yet
+  const Int_t maxNDigits = 1500;
   Int_t nDigits=fDigitsArr->GetEntriesFast() ;
+  if (nDigits > maxNDigits) {
+    AliWarning(Form("Skip event with too high digit occupancy: nDigits=%d",nDigits));
+    return;
+  }
+
   for(Int_t i=0; i<nDigits; i++){
     fDigitsUsed[i]=0 ;
   }
   Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
                      //e.g. first digit in this module, first CPV digit etc.
   AliPHOSDigit * digit ; 
-  TArrayI clusterdigitslist(1500) ;   
+  TArrayI clusterdigitslist(maxNDigits) ;   
   AliPHOSRecPoint * clu = 0 ; 
   for(Int_t i=0; i<nDigits; i++){
     if(fDigitsUsed[i])
@@ -1049,15 +1057,22 @@ void AliPHOSClusterizerv1::SetDistancesToBadChannels()
   TVector3 dR;
 
   Float_t dist,minDist;
-  Int_t relid[4] ;
-
+  Int_t relid[4]={0,0,0,0} ;
+  TVector3 lpos ;
   for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
     rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
-    minDist = 1.e+07;
+    //evaluate distance to border
+    relid[0]=rp->GetPHOSMod() ;
+    relid[2]=1 ;
+    relid[3]=1 ;
+    Float_t xcorner,zcorner;
+    fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
+    rp->GetLocalPosition(lpos) ;
+    minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
     for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
       fGeom->AbsToRelNumbering(badIds[iBad],relid)  ;
-      if(relid[0]!=rp->GetPHOSMod())
-        continue ;
+      if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since 
+        continue ;                   //bad channels can be in the module which does not exist in simulations.
       rp->GetGlobalPosition(gposRecPoint,gmat);
       fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
       AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",