new unfolding for clusterization
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 17 Oct 2010 20:00:01 +0000 (20:00 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 17 Oct 2010 20:00:01 +0000 (20:00 +0000)
14 files changed:
EMCAL/AliCaloCalibPedestal.cxx
EMCAL/AliEMCALClusterizer.cxx
EMCAL/AliEMCALClusterizer.h
EMCAL/AliEMCALClusterizerNxN.cxx
EMCAL/AliEMCALClusterizerNxN.h
EMCAL/AliEMCALClusterizerv1.cxx
EMCAL/AliEMCALClusterizerv1.h
EMCAL/AliEMCALRecParam.cxx
EMCAL/AliEMCALRecParam.h
EMCAL/AliEMCALReconstructor.cxx
EMCAL/AliEMCALUnfolding.cxx [new file with mode: 0644]
EMCAL/AliEMCALUnfolding.h [new file with mode: 0644]
EMCAL/EMCALrecLinkDef.h
EMCAL/libEMCALrec.pkg

index 44e9d16..f527aaf 100644 (file)
@@ -254,30 +254,37 @@ void AliCaloCalibPedestal::CompressAndSetOwner()
 //_____________________________________________________________________
 AliCaloCalibPedestal::~AliCaloCalibPedestal()
 {
-  if (fReference) delete fReference;//Delete the reference object, if it has been loaded
+  //dtor
+  printf("Dtor\n");
 
-  // delete also TObjArray's 
-  fPedestalLowGain.Delete();
-  fPedestalHighGain.Delete();
-  fPedestalLEDRefLowGain.Delete();
-  fPedestalLEDRefHighGain.Delete();
-  fPeakMinusPedLowGain.Delete();
-  fPeakMinusPedHighGain.Delete();
-  fPeakMinusPedHighGainHisto.Delete();
-  fPedestalLowGainDiff.Delete();
-  fPedestalHighGainDiff.Delete();
-  fPedestalLEDRefLowGainDiff.Delete();
-  fPedestalLEDRefHighGainDiff.Delete();
-  fPeakMinusPedLowGainDiff.Delete();
-  fPeakMinusPedHighGainDiff.Delete();
-  fPedestalLowGainRatio.Delete();
-  fPedestalHighGainRatio.Delete();
-  fPedestalLEDRefLowGainRatio.Delete();
-  fPedestalLEDRefHighGainRatio.Delete();
-  fPeakMinusPedLowGainRatio.Delete();
-  fPeakMinusPedHighGainRatio.Delete();
-  fDeadMap.Delete();
+  
+  if (fReference){   printf("Ref\n"); delete fReference;}//Delete the reference object, if it has been loaded
+  
 
+  printf("Delete\n");
+
+  // delete also TObjArray's 
+  fPedestalLowGain.Delete();  printf("D 1\n");
+  fPedestalHighGain.Delete();printf("D 2\n");
+  fPedestalLEDRefLowGain.Delete();printf("D 3\n");
+  fPedestalLEDRefHighGain.Delete();printf("D 4\n");
+  fPeakMinusPedLowGain.Delete();printf("D 5\n");
+  fPeakMinusPedHighGain.Delete();printf("D 6\n");
+  fPeakMinusPedHighGainHisto.Delete();printf("D 7\n");
+  fPedestalLowGainDiff.Delete();printf("D 8\n");
+  fPedestalHighGainDiff.Delete();printf("D 9\n");
+  fPedestalLEDRefLowGainDiff.Delete();printf("D 10\n");
+  fPedestalLEDRefHighGainDiff.Delete();printf("D 11\n");
+  fPeakMinusPedLowGainDiff.Delete();printf("D 12\n");
+  fPeakMinusPedHighGainDiff.Delete();printf("D 13\n");
+  fPedestalLowGainRatio.Delete();printf("D 14\n");
+  fPedestalHighGainRatio.Delete();printf("D 15\n");
+  fPedestalLEDRefLowGainRatio.Delete();printf("D 16\n");
+  fPedestalLEDRefHighGainRatio.Delete();printf("D 17\n");
+  fPeakMinusPedLowGainRatio.Delete();printf("D 18\n");
+  fPeakMinusPedHighGainRatio.Delete();printf("D 19\n");
+  fDeadMap.Delete();printf("D 20\n");
+  
 }
 
 // copy ctor
@@ -549,6 +556,8 @@ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
 { 
   // Method to process=analyze one event in the data stream
   if (!in) return kFALSE; //Return right away if there's a null pointer
+  in->Reset(); // just in case the next customer forgets to check if the stream was reset..
+
   fNEvents++; // one more event
 
   if (fNEvents==1) ValidateProfiles(); // 1st event, make sure histos/profiles exist
@@ -652,7 +661,6 @@ Bool_t AliCaloCalibPedestal::ProcessEvent(AliCaloRawStreamV3 *in)
     }// end while over channel   
   }//end while over DDL's, of input stream 
 
-  in->Reset(); // just in case the next customer forgets to check if the stream was reset..
  
   return kTRUE;
 }
index 2963c2d..7675f59 100644 (file)
@@ -56,6 +56,7 @@ class TSystem;
 #include "AliEMCALCalibData.h"
 class AliCDBStorage;
 #include "AliCDBEntry.h"
+#include "AliEMCALUnfolding.h"
 
 ClassImp(AliEMCALClusterizer)
 
@@ -71,7 +72,8 @@ AliEMCALClusterizer::AliEMCALClusterizer():
   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
   fDefaultInit(kFALSE),fToUnfold(kFALSE),
   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
-  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
+  fClusterUnfolding(NULL)
 {
   // ctor
   
@@ -91,7 +93,8 @@ AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
   fDefaultInit(kFALSE),fToUnfold(kFALSE),
   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
-  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
+  fClusterUnfolding(NULL)
 {
   // ctor with the indication of the file where header Tree and digits Tree are stored
   // use this contructor to avoid usage of Init() which uses runloader
@@ -106,7 +109,13 @@ AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
   {
     AliFatal("Geometry not initialized.");
   }
-  
+  Int_t i=0;
+  for (i = 0; i < 8; i++)
+    fSSPars[i] = 0.;
+  for (i = 0; i < 3; i++) {
+    fPar5[i] = 0.;
+    fPar6[i] = 0.;
+  }
 }
 
 //____________________________________________________________________________
@@ -121,13 +130,22 @@ AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCal
   fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
   fDefaultInit(kFALSE),fToUnfold(kFALSE),
   fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
-  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.)
+  fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
+  fClusterUnfolding(NULL)
 {
        // ctor, geometry and calibration are initialized elsewhere.
        
        if (!fGeom)
                AliFatal("Geometry not initialized.");
   
+  Int_t i=0;
+  for (i = 0; i < 8; i++)
+    fSSPars[i] = 0.;
+  for (i = 0; i < 3; i++) {
+    fPar5[i] = 0.;
+    fPar6[i] = 0.;
+  }
+  
 }
 
 
@@ -149,6 +167,8 @@ AliEMCALClusterizer::~AliEMCALClusterizer()
 //     fRecPoints->Delete();
 //     delete fRecPoints;
 //  }
+  
+  if(fClusterUnfolding) delete fClusterUnfolding;
 }
 
 //____________________________________________________________________________
@@ -273,6 +293,14 @@ void AliEMCALClusterizer::Init()
   if(!gMinuit) 
     gMinuit = new TMinuit(100) ;
   
+  Int_t i=0;
+  for (i = 0; i < 8; i++)
+    fSSPars[i] = 0.;
+  for (i = 0; i < 3; i++) {
+    fPar5[i] = 0.;
+    fPar6[i] = 0.;
+  }
+  
 }
 
 //____________________________________________________________________________
@@ -301,6 +329,28 @@ void AliEMCALClusterizer::InitParameters()
                     fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
   }
   
+  if(fToUnfold){
+    Int_t i=0;
+    for (i = 0; i < 8; i++) {
+      fSSPars[i] = recParam->GetSSPars(i);
+    }//end of loop over parameters
+    for (i = 0; i < 3; i++) {
+      fPar5[i] = recParam->GetPar5(i);
+      fPar6[i] = recParam->GetPar6(i);
+    }//end of loop over parameters
+    
+    fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6);
+    
+    for (i = 0; i < 8; i++) {
+      AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
+    }
+    for (i = 0; i < 3; i++) {
+      AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
+      AliDebug(1,Form("unfolding parameter 6: fPar5=%f \n",fPar6[i]));
+    }
+    
+  }
+
 }
 
 
index b320ca0..438cd9d 100644 (file)
@@ -22,6 +22,7 @@ class TTree;
 class AliEMCALGeometry ;
 class AliEMCALCalibData ;
 class AliCaloCalibPedestal ;
+class AliEMCALUnfolding ;
 
 class AliEMCALClusterizer : public TObject {
 
@@ -68,6 +69,8 @@ public:
   virtual void Print(Option_t * option)const ;
   virtual void PrintRecPoints(Option_t * option);
   virtual void PrintRecoInfo();                        //*MENU*
+
+
   
   virtual const char * Version() const {Warning("Version", "Not Defined") ; return 0 ; } 
 
@@ -99,13 +102,18 @@ protected:
   Float_t fECAW0 ;                   // logarithmic weight for the cluster center of gravity calculation
   Float_t fMinECut;                  // Minimum energy for a digit to be a member of a cluster
   
+  AliEMCALUnfolding * fClusterUnfolding ; //! pointer to unfolding object
+  Double_t fSSPars[8];// Shower shape parameters 
+  Double_t fPar5[3];  // Shower shape parameter 5
+  Double_t fPar6[3];  // Shower shape parameter 6
+
 private:
   AliEMCALClusterizer(const AliEMCALClusterizer &); //copy ctor
   AliEMCALClusterizer & operator = (const AliEMCALClusterizer &);
   
   
   
-  ClassDef(AliEMCALClusterizer,2)  // Clusterization algorithm class 
+  ClassDef(AliEMCALClusterizer,3)  // Clusterization algorithm class 
 } ;
 
 #endif // AliEMCALCLUSTERIZER_H
index cca2fc8..75c5d92 100644 (file)
@@ -61,6 +61,7 @@
 #include "AliCaloCalibPedestal.h"
 #include "AliEMCALCalibData.h"
 #include "AliESDCaloCluster.h"
+#include "AliEMCALUnfolding.h"
 
 ClassImp(AliEMCALClusterizerNxN)
 
@@ -121,12 +122,13 @@ void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
   
   MakeClusters() ;  //only the real clusters
   
-  if(fToUnfold)
-    MakeUnfolding() ;
+  if(fToUnfold){
+    fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
+    fClusterUnfolding->MakeUnfolding();
+  }
   
+  //Evaluate position, dispersion and other RecPoint properties for EC section 
   Int_t index ;
-  
-  //Evaluate position, dispersion and other RecPoint properties for EC section                      
   for(index = 0; index < fRecPoints->GetEntries(); index++) 
   { 
     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
@@ -358,85 +360,6 @@ void AliEMCALClusterizerNxN::MakeClusters()
   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
 }
 
-//____________________________________________________________________________
-void AliEMCALClusterizerNxN::MakeUnfolding()
-{
-  // Unfolds clusters using the shape of an ElectroMagnetic shower
-  // Performs unfolding of all clusters
-  
-  if(fNumberOfECAClusters > 0){
-    
-    if (fGeom){
-      
-      Int_t nModulesToUnfold = fGeom->GetNCells();
-      
-      Int_t numberofNotUnfolded = fNumberOfECAClusters ;
-      Int_t index ;
-      for(index = 0 ; index < numberofNotUnfolded ; index++){
-        
-        AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
-        if(recPoint){
-          TVector3 gpos;
-          Int_t absId;
-          recPoint->GetGlobalPosition(gpos);
-          fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
-          if(absId > nModulesToUnfold)
-            break ;
-          
-          Int_t nMultipl = recPoint->GetMultiplicity() ;
-          AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
-          Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-          Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
-          
-          if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
-            //UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
-            fRecPoints->Remove(recPoint);
-            fRecPoints->Compress() ;
-            index-- ;
-            fNumberOfECAClusters-- ;
-            numberofNotUnfolded-- ;
-          }
-          else{
-            recPoint->SetNExMax(1) ; //Only one local maximum
-          }
-          
-          delete[] maxAt ;
-          delete[] maxAtEnergy ;
-        }
-        else AliFatal("Could not get RecPoint!") ;
-      }//loop
-    }
-    else AliFatal("Could not get geometry from EMCALLoader!") ;
-
-  }
-  // End of Unfolding of clusters
-}
-
-//____________________________________________________________________________
-Double_t  AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
-{ 
-  // Shape of the shower
-  // If you change this function, change also the gradient evaluation in ChiSquare()
-  //
-  Double_t r = sqrt(x*x+y*y);
-  Double_t r133  = TMath::Power(r, 1.33) ;
-  Double_t r669  = TMath::Power(r, 6.69) ;
-  Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
-  return shape ;
-}
-
-//____________________________________________________________________________
-void  AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, 
-                                           Int_t /*nMax*/, 
-                                           AliEMCALDigit ** /*maxAt*/, 
-                                           Float_t * /*maxAtEnergy*/)
-{
-  //
-  // Performs the unfolding of a cluster with nMax overlapping showers 
-  //
-  AliWarning("Not implemented. To be.");
-}
-
 //___________________________________________________________________
 void   AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
 {
index 1d4c893..4e5457e 100644 (file)
@@ -33,15 +33,11 @@ public:
   virtual Int_t   AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared)const ; 
                                // Checks if digits are in neighbour cells 
 
-   virtual void   Digits2Clusters(Option_t *option);                // Does the job
+  virtual void   Digits2Clusters(Option_t *option);                // Does the job
 
-
-  static Double_t ShowerShape(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; 
-                                            //class member function (not object member function)
   static void     SetInputCalibrated(Bool_t val);
 
   virtual const char * Version() const { return "clu-NxN" ; }  
   
 protected:
 
@@ -53,14 +49,6 @@ private:
   AliEMCALClusterizerNxN(const AliEMCALClusterizerNxN &); //copy ctor
   AliEMCALClusterizerNxN & operator = (const AliEMCALClusterizerNxN &);
 
-
-  virtual void   MakeUnfolding();
-  void           UnfoldCluster(AliEMCALRecPoint * iniEmc, Int_t Nmax, 
-                               AliEMCALDigit ** maxAt,
-                               Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package
-
-private:
-
    ClassDef(AliEMCALClusterizerNxN,1)   // Clusterizer implementation version 1
 
 };
index 2aede01..8f6c4e0 100644 (file)
@@ -49,6 +49,7 @@
 #include "AliCaloCalibPedestal.h"
 #include "AliEMCALCalibData.h"
 #include "AliESDCaloCluster.h"
+#include "AliEMCALUnfolding.h"
 
 ClassImp(AliEMCALClusterizerv1)
 
@@ -107,12 +108,13 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
   
   MakeClusters() ;  //only the real clusters
   
-  if(fToUnfold)
-    MakeUnfolding() ;
-  
+  if(fToUnfold){
+    fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
+    fClusterUnfolding->MakeUnfolding();
+  }
+    
+  //Evaluate position, dispersion and other RecPoint properties for EC section 
   Int_t index ;
-  
-  //Evaluate position, dispersion and other RecPoint properties for EC section                      
   for(index = 0; index < fRecPoints->GetEntries(); index++) {
     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
     if(rp){
@@ -151,99 +153,6 @@ void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
 }
 
 //____________________________________________________________________________
-Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, 
-                                     const Float_t* maxAtEnergy,
-                                     Int_t nPar, Float_t * fitparameters) const
-{
-  // Calls TMinuit to fit the energy distribution of a cluster with several maxima
-  // The initial values for fitting procedure are set equal to the
-  // positions of local maxima.       
-  // Cluster will be fitted as a superposition of nPar/3
-  // electromagnetic showers
-
-  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
-       
-  if(!gMinuit)
-     gMinuit = new TMinuit(100) ;
-
-  gMinuit->mncler();                     // Reset Minuit's list of paramters
-  gMinuit->SetPrintLevel(-1) ;           // No Printout
-  gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
-  // To set the address of the minimization function
-  TList * toMinuit = new TList();
-  toMinuit->AddAt(recPoint,0) ;
-  toMinuit->AddAt(fDigitsArr,1) ;
-  toMinuit->AddAt(fGeom,2) ;
-
-  gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare
-
-  // filling initial values for fit parameters
-  AliEMCALDigit * digit ;
-
-  Int_t ierflg  = 0;
-  Int_t index   = 0 ;
-  Int_t nDigits = (Int_t) nPar / 3 ;
-
-  Int_t iDigit ;
-
-  for(iDigit = 0; iDigit < nDigits; iDigit++){
-    digit = maxAt[iDigit];
-    Double_t x = 0.;
-    Double_t y = 0.;
-    Double_t z = 0.;
-
-    fGeom->RelPosCellInSModule(digit->GetId(), y, x, z);
-
-    Float_t energy = maxAtEnergy[iDigit] ;
-
-    gMinuit->mnparm(index, "x",  x, 0.1, 0, 0, ierflg) ;
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %f", x ) ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index, "z",  z, 0.1, 0, 0, ierflg) ;
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %f", z) ;
-      return kFALSE;
-    }
-    gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;
-    index++ ;
-    if(ierflg != 0){
-      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;
-      return kFALSE;
-    }
-  }
-
-  Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; 
-                      // The number of function call slightly depends on it.
-  //Double_t p1 = 1.0 ;
-  Double_t p2 = 0.0 ;
-
-  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls
-  //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient
-  gMinuit->SetMaxIterations(5);
-  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
-  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize
-
-  if(ierflg == 4){  // Minimum not found
-    Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;
-    return kFALSE ;
-  }
-  for(index = 0; index < nPar; index++){
-    Double_t err = 0. ;
-    Double_t val = 0. ;
-    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index
-    fitparameters[index] = val ;
-  }
-
-  delete toMinuit ;
-  return kTRUE;
-
-}
-
-//____________________________________________________________________________
 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 {
        // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching 
@@ -386,269 +295,3 @@ void AliEMCALClusterizerv1::MakeClusters()
   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); 
 }
 
-//____________________________________________________________________________
-void AliEMCALClusterizerv1::MakeUnfolding()
-{
-  // Unfolds clusters using the shape of an ElectroMagnetic shower
-  // Performs unfolding of all clusters
-  
-  if(fNumberOfECAClusters > 0){
-    if (fGeom==0)
-      AliFatal("Did not get geometry from EMCALLoader") ;
-    Int_t nModulesToUnfold = fGeom->GetNCells();
-    
-    Int_t numberofNotUnfolded = fNumberOfECAClusters ;
-    Int_t index ;
-    for(index = 0 ; index < numberofNotUnfolded ; index++){
-      
-      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
-      if(recPoint){
-        TVector3 gpos;
-        Int_t absId = -1;
-        recPoint->GetGlobalPosition(gpos);
-        fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
-        if(absId > nModulesToUnfold)
-          break ;
-        
-        Int_t nMultipl = recPoint->GetMultiplicity() ;
-        AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
-        Float_t * maxAtEnergy = new Float_t[nMultipl] ;
-        Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
-        
-        if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
-          UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
-          fRecPoints->Remove(recPoint);
-          fRecPoints->Compress() ;
-          index-- ;
-          fNumberOfECAClusters-- ;
-          numberofNotUnfolded-- ;
-        }
-        else{
-          recPoint->SetNExMax(1) ; //Only one local maximum
-        }
-        
-        delete[] maxAt ;
-        delete[] maxAtEnergy ;
-      }
-      else AliFatal("Null recpoint in Array!");
-    }
-  }
-  // End of Unfolding of clusters
-}
-
-//____________________________________________________________________________
-Double_t  AliEMCALClusterizerv1::ShowerShape(Double_t x, Double_t y)
-{ 
-  // Shape of the shower
-  // If you change this function, change also the gradient evaluation in ChiSquare()
-
-  Double_t r = sqrt(x*x+y*y);
-  Double_t r133  = TMath::Power(r, 1.33) ;
-  Double_t r669  = TMath::Power(r, 6.69) ;
-  Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
-  return shape ;
-}
-
-//____________________________________________________________________________
-void  AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, 
-                                          Int_t nMax, 
-                                          AliEMCALDigit ** maxAt, 
-                                          Float_t * maxAtEnergy)
-{
-  // Performs the unfolding of a cluster with nMax overlapping showers 
-  Int_t nPar = 3 * nMax ;
-  Float_t * fitparameters = new Float_t[nPar] ;
-  
-  if (fGeom==0)
-    AliFatal("Did not get geometry from EMCALLoader") ;
-  
-  Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;
-  if( !rv ) {
-    // Fit failed, return and remove cluster
-    iniTower->SetNExMax(-1) ;
-    delete[] fitparameters ;
-    return ;
-  }
-  
-  // create unfolded rec points and fill them with new energy lists
-  // First calculate energy deposited in each sell in accordance with
-  // fit (without fluctuations): efit[]
-  // and later correct this number in acordance with actual energy
-  // deposition
-  
-  Int_t nDigits = iniTower->GetMultiplicity() ;
-  Float_t * efit = new Float_t[nDigits] ;
-  Double_t xDigit=0.,yDigit=0.,zDigit=0. ;
-  Float_t xpar=0.,zpar=0.,epar=0.  ;
-  
-  AliEMCALDigit * digit = 0 ;
-  Int_t * digitsList = iniTower->GetDigitsList() ;
-  
-  Int_t iparam = 0 ;
-  Int_t iDigit ;
-  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;
-    if(digit){
-      fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
-      efit[iDigit] = 0;
-      
-      while(iparam < nPar ){
-        xpar = fitparameters[iparam] ;
-        zpar = fitparameters[iparam+1] ;
-        epar = fitparameters[iparam+2] ;
-        iparam += 3 ;
-        efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
-      }
-    }
-    else AliFatal("Null digit in array!");
-  }
-  
-  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
-  // so that energy deposited in each cell is distributed between new clusters proportionally
-  // to its contribution to efit
-  
-  Float_t * energiesList = iniTower->GetEnergiesList() ;
-  Float_t ratio = 0 ;
-  
-  iparam = 0 ;
-  while(iparam < nPar ){
-    xpar = fitparameters[iparam] ;
-    zpar = fitparameters[iparam+1] ;
-    epar = fitparameters[iparam+2] ;
-    iparam += 3 ;
-    
-    AliEMCALRecPoint * recPoint = 0 ;
-    
-    if(fNumberOfECAClusters >= fRecPoints->GetSize())
-      fRecPoints->Expand(2*fNumberOfECAClusters) ;
-    
-    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
-    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
-    if(recPoint){
-      fNumberOfECAClusters++ ;
-      recPoint->SetNExMax((Int_t)nPar/3) ;
-      
-      Float_t eDigit = 0. ;
-      for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
-        digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;
-        if(digit){
-          fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
-          
-          ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
-          eDigit = energiesList[iDigit] * ratio ;
-          recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case
-        }
-        else AliFatal("Null digit in array!");
-      }
-    }
-    else AliFatal("Null recpoint in array!");
-  }
-  
-  delete[] fitparameters ;
-  delete[] efit ;
-  
-}
-
-//_____________________________________________________________________________
-void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad,
-                                              Double_t & fret,
-                                              Double_t * x, Int_t iflag)
-{
-  // Calculates the Chi square for the cluster unfolding minimization
-  // Number of parameters, Gradient, Chi squared, parameters, what to do
-  
-  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
-  if(!toMinuit){
-    printf("Unfolding not possible!\n");
-    return;
-  }
-  
-  AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;
-  TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;
-  // A bit buggy way to get an access to the geometry
-  // To be revised!
-  AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));
-  
-  if(!recPoint || !digits || !geom){
-    printf("Unfolding not possible!\n");
-    return;
-  }
-  
-  Int_t * digitsList     = recPoint->GetDigitsList() ;
-  
-  Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;
-  
-  Float_t * energiesList = recPoint->GetEnergiesList() ;
-  
-  fret = 0. ;
-  Int_t iparam ;
-  
-  if(iflag == 2)
-    for(iparam = 0 ; iparam < nPar ; iparam++)
-      Grad[iparam] = 0 ; // Will evaluate gradient
-  
-  Double_t efit = 0. ;
-  
-  AliEMCALDigit * digit ;
-  Int_t iDigit ;
-  
-  for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
-    
-    digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );
-    if (digit) {
-      
-      Double_t xDigit=0 ;
-      Double_t zDigit=0 ;
-      Double_t yDigit=0 ;//not used yet, assumed to be 0
-      
-      geom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit);
-      
-      if(iflag == 2){  // calculate gradient
-        Int_t iParam = 0 ;
-        efit = 0. ;
-        while(iParam < nPar ){
-          Double_t dx = (xDigit - x[iParam]) ;
-          iParam++ ;
-          Double_t dz = (zDigit - x[iParam]) ;
-          iParam++ ;
-          efit += x[iParam] * ShowerShape(dx,dz) ;
-          iParam++ ;
-        }
-        Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)
-        iParam = 0 ;
-        while(iParam < nPar ){
-          Double_t xpar = x[iParam] ;
-          Double_t zpar = x[iParam+1] ;
-          Double_t epar = x[iParam+2] ;
-          Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
-          Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
-          Double_t r133 =  TMath::Power(dr, 1.33);
-          Double_t r669 = TMath::Power(dr,6.69);
-          Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )
-                                                              - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;
-          
-          Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ;  // Derivative over x
-          iParam++ ;
-          Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ;  // Derivative over z
-          iParam++ ;
-          Grad[iParam] += shape ;                                  // Derivative over energy
-          iParam++ ;
-        }
-      }
-      efit = 0;
-      iparam = 0 ;
-      
-      
-      while(iparam < nPar ){
-        Double_t xpar = x[iparam] ;
-        Double_t zpar = x[iparam+1] ;
-        Double_t epar = x[iparam+2] ;
-        iparam += 3 ;
-        efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
-      }
-      
-      fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;
-      // Here we assume, that sigma = sqrt(E) 
-    }
-  }
-}
index f7b94eb..fe7151e 100644 (file)
@@ -38,10 +38,6 @@ public:
                                // Checks if digits are in neighbour cells
   virtual void    Digits2Clusters(Option_t *option);                // Does the job
 
-  static Double_t ShowerShape(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; 
-                                            //class member function (not object member function)
-  static void UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;
-                                            // Chi^2 of the fit. Should be static to be passes to MINUIT
   virtual const char * Version() const { return "clu-v1" ; }  
 
 protected:
@@ -52,18 +48,6 @@ private:
   AliEMCALClusterizerv1(const AliEMCALClusterizerv1 &); //copy ctor
   AliEMCALClusterizerv1 & operator = (const AliEMCALClusterizerv1 &);
 
-
-  Bool_t  FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** MaxAt, const Float_t * maxAtEnergy, 
-                 Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit
-
-  virtual void   MakeUnfolding();
-  void           UnfoldCluster(AliEMCALRecPoint * iniEmc, Int_t Nmax, 
-                               AliEMCALDigit ** maxAt,
-                               Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package
-
-private:
-
-
   ClassDef(AliEMCALClusterizerv1,10)   // Clusterizer implementation version 1
 
 };
index ca73559..675911c 100644 (file)
@@ -221,7 +221,22 @@ AliEMCALRecParam::AliEMCALRecParam() :
   fHadronEnergyProb[3]=  2.030230e+00;
   fHadronEnergyProb[4]= -6.402242e-02;
   
-  
+  //unfolding  
+  fSSPars[0] = 0.9262;
+  fSSPars[1] = 3.365;
+  fSSPars[2] = 1.548;
+  fSSPars[3] = 0.1625;
+  fSSPars[4] = -0.4195;
+  fSSPars[5] = 0.;
+  fSSPars[6] = 0.;
+  fSSPars[7] = 2.332;
+  fPar5[0] = 12.31;
+  fPar5[1] = -0.007381;
+  fPar5[2] = -0.06936;
+  fPar6[0] = 0.05452; 
+  fPar6[1] = 0.0001228; 
+  fPar6[2] = 0.001361; 
+
 }
 
 //-----------------------------------------------------------------------------
@@ -273,6 +288,15 @@ AliEMCALRecParam::AliEMCALRecParam(const AliEMCALRecParam& rp) :
     
   }
   
+  //unfolding  
+  for (i = 0; i < 8; i++) {
+    fSSPars[i] = rp.fSSPars[i];
+  }
+  for (i = 0; i < 3; i++) {
+    fPar5[i] = rp.fPar5[i];
+    fPar6[i] = rp.fPar6[i];
+  }
+
 }
 
 //-----------------------------------------------------------------------------
@@ -323,7 +347,15 @@ AliEMCALRecParam& AliEMCALRecParam::operator = (const AliEMCALRecParam& rp)
       fPiZeroEnergyProb[i] = rp.fPiZeroEnergyProb[i];
       fHadronEnergyProb[i] = rp.fHadronEnergyProb[i];
     }
-    
+    //unfolding  
+    for (i = 0; i < 8; i++) {
+      fSSPars[i] = rp.fSSPars[i];
+    }
+    for (i = 0; i < 3; i++) {
+      fPar5[i] = rp.fPar5[i];
+      fPar6[i] = rp.fPar6[i];
+    }
+
   }    
   
   return *this;
@@ -563,6 +595,20 @@ void AliEMCALRecParam::Print(Option_t * opt) const
     
     AliInfo(Form("Track-matching cuts :\n x %f, y %f, z %f, R %f \n alphaMin %f, alphaMax %f, Angle %f, NITS %f, NTPC %f\n", 
                                 fTrkCutX, fTrkCutY, fTrkCutZ, fTrkCutR,fTrkCutAlphaMin,fTrkCutAlphaMax, fTrkCutAngle,fTrkCutNITS,fTrkCutNTPC));
+
+    AliInfo(Form("Unfolding parameters, Shower shape function :\n")); 
+    for(Int_t i = 0; i < 8; i++){
+       printf(" %f, ", fSSPars[i]);
+    }
+    printf("\n Parameter 5 : ");
+    for(Int_t i = 0; i < 3; i++){
+      printf(" %f, ", fPar5[i]);
+    }
+    printf("\n Parameter 6 : ");
+    for(Int_t i = 0; i < 3; i++){
+      printf(" %f, ", fPar6[i]);
+    }
+    printf("\n");
   }
   
   if(!strcmp("",opt) || !strcmp("pid",opt)){
index 3e09c71..3e2be8b 100644 (file)
@@ -119,6 +119,15 @@ class AliEMCALRecParam : public AliDetectorRecoParam
   Bool_t   UseFALTRO()            const {return fUseFALTRO; }
   Bool_t   FitLEDEvents()         const {return fFitLEDEvents; }
 
+  //Unfolding (Adam)
+  Double_t GetSSPars(Int_t i) const   {return fSSPars[i];}
+  Double_t GetPar5(Int_t i) const     {return fPar5[i];}
+  Double_t GetPar6(Int_t i) const     {return fPar6[i];}
+  void SetSSPars(Int_t i, Double_t param )     {fSSPars[i]=param;}
+  void SetPar5(Int_t i, Double_t param )       {fPar5[i]=param;}
+  void SetPar6(Int_t i, Double_t param )       {fPar6[i]=param;}
+
+
   virtual void Print(Option_t * option="") const ;
   
   static AliEMCALRecParam* GetDefaultParameters();
@@ -177,9 +186,14 @@ class AliEMCALRecParam : public AliDetectorRecoParam
   Bool_t   fUseFALTRO;             // get FALTRO (trigger) and put it on trigger digits.
   Bool_t   fFitLEDEvents;          // fit LED events or not
        
+  //Shower shape parameters (Adam)
+  Double_t fSSPars[8]; // Unfolding shower shape parameters
+  Double_t fPar5[3];   // UF SSPar nr 5
+  Double_t fPar6[3];   // UF SSPar nr 6
+
   static TObjArray* fgkMaps;       // ALTRO mappings for RCU0..RCUX
   
-  ClassDef(AliEMCALRecParam,13)     // Reconstruction parameters
+  ClassDef(AliEMCALRecParam,14)     // Reconstruction parameters
     
     } ;
 
index 566cf6f..8895b1f 100644 (file)
@@ -124,14 +124,14 @@ AliEMCALReconstructor::AliEMCALReconstructor()
        
   if(!fGeom) AliFatal(Form("Could not get geometry!"));
 
-  AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
+  //AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
 
-  const AliEMCALTriggerDCSConfig* dcsConfig = dcsConfigDB->GetTriggerDCSConfig();
+  //const AliEMCALTriggerDCSConfig* dcsConfig = dcsConfigDB->GetTriggerDCSConfig();
 
-  if (!dcsConfig) AliFatal("No Trigger DCS Configuration from OCDB!");
-  fgTriggerProcessor = new AliEMCALTriggerElectronics( dcsConfig );
+  //if (!dcsConfig) AliFatal("No Trigger DCS Configuration from OCDB!");
+  //fgTriggerProcessor = new AliEMCALTriggerElectronics( dcsConfig );
        
-  fTriggerData = new AliEMCALTriggerData();
+  //fTriggerData = new AliEMCALTriggerData();
 
  //Init temporary list of digits
   fgDigitsArr   = new TClonesArray("AliEMCALDigit",1000);
@@ -160,7 +160,7 @@ AliEMCALReconstructor::~AliEMCALReconstructor()
   
   if(fgRawUtils)         delete fgRawUtils;
   if(fgClusterizer)      delete fgClusterizer;
-  if(fgTriggerProcessor) delete fgTriggerProcessor;
+  //if(fgTriggerProcessor) delete fgTriggerProcessor;
   
   AliCodeTimer::Instance()->Print();
 } 
diff --git a/EMCAL/AliEMCALUnfolding.cxx b/EMCAL/AliEMCALUnfolding.cxx
new file mode 100644 (file)
index 0000000..f9cd63b
--- /dev/null
@@ -0,0 +1,617 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//_________________________________________________________________________\r
+//  Base class for the cluster unfolding algorithm \r
+//*-- Author: Adam Matyja (SUBATECH)\r
+//  Based on unfolding in clusterizerv1 done by Cynthia Hadjidakis\r
+//-- Unfolding for eta~0: Cynthia Hadjidakis - still in AliEMCALCLusterizerv1\r
+//-- Unfolding extension for whole EMCAL: Adam Matyja (SUBATECH & INP PAN)\r
+//\r
+//  unfolds the clusters having several local maxima. \r
+//////////////////////////////////////////////////////////////////////////////\r
+\r
+// --- ROOT system ---\r
+#include "TClonesArray.h"\r
+//#include "TTree.h"\r
+//#include <TFile.h> \r
+//class TFolder;\r
+#include <TMath.h> \r
+#include <TMinuit.h>\r
+//#include <TTree.h> \r
+//class TSystem; \r
+//#include <TBenchmark.h>\r
+//#include <TBrowser.h>\r
+//#include <TROOT.h>\r
+\r
+// --- Standard library ---\r
+#include <cassert>\r
+\r
+// --- AliRoot header files ---\r
+#include "AliEMCALUnfolding.h"\r
+#include "AliEMCALGeometry.h"\r
+#include "AliRunLoader.h"\r
+#include "AliRun.h"\r
+#include "AliEMCAL.h"\r
+#include "AliEMCALRecParam.h"\r
+#include "AliEMCALRecPoint.h"\r
+#include "AliEMCALDigit.h"\r
+#include "AliEMCALReconstructor.h"\r
+//#include "AliEMCALClusterizer.h"\r
+\r
+\r
+\r
+#include "AliLog.h"\r
+\r
+#include "AliCDBManager.h"\r
+//#include "AliCaloCalibPedestal.h"\r
+//#include "AliEMCALCalibData.h"\r
+class AliCDBStorage;\r
+#include "AliCDBEntry.h"\r
+\r
+Double_t AliEMCALUnfolding::fSSPars[8]={0.9262,3.365,1.548,0.1625,-0.4195,0.,0.,2.332};\r
+Double_t AliEMCALUnfolding::fPar5[3]={12.31,-0.007381,-0.06936};\r
+Double_t AliEMCALUnfolding::fPar6[3]={0.05452,0.0001228,0.001361};\r
+\r
+ClassImp(AliEMCALUnfolding)\r
+  \r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding():\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(0),\r
+  fGeom(NULL),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
+  \r
+  Init() ;\r
+}\r
+\r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry):\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(0),\r
+  fGeom(geometry),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
+  // use this contructor to avoid usage of Init() which uses runloader\r
+  // change needed by HLT - MP\r
+  if (!fGeom)\r
+  {\r
+    AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
+  }\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+AliEMCALUnfolding::AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6):\r
+  fNumberOfECAClusters(0),\r
+  fECALocMaxCut(ECALocMaxCut),\r
+  fGeom(geometry),\r
+  fRecPoints(NULL),\r
+  fDigitsArr(NULL)\r
+{\r
+  // ctor with the indication of the file where header Tree and digits Tree are stored\r
+  // use this contructor to avoid usage of Init() which uses runloader\r
+  // change needed by HLT - MP\r
+  if (!fGeom)\r
+  {\r
+    AliFatal("AliEMCALUnfolding: Geometry not initialized.");\r
+  }\r
+  Int_t i=0;\r
+  for (i = 0; i < 8; i++) fSSPars[i] = SSPars[i];\r
+  for (i = 0; i < 3; i++) {\r
+    fPar5[i] = Par5[i];\r
+    fPar6[i] = Par6[i];\r
+  }\r
+\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::Init()\r
+{\r
+  // Make all memory allocations which can not be done in default constructor.\r
+  // Attach the Clusterizer task to the list of EMCAL tasks\r
+  \r
+  AliRunLoader *rl = AliRunLoader::Instance();\r
+  if (rl->GetAliRun()){\r
+    AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));\r
+    if(emcal)fGeom = emcal->GetGeometry();\r
+  }\r
+  \r
+  if(!fGeom)\r
+    fGeom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());\r
+  \r
+  AliDebug(1,Form("geom %p",fGeom));\r
+  \r
+  if(!gMinuit) \r
+    gMinuit = new TMinuit(100) ;\r
+  \r
+}\r
+\r
+//____________________________________________________________________________\r
+  AliEMCALUnfolding::~AliEMCALUnfolding()\r
+{\r
+  // dtor\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr)\r
+{\r
+  //\r
+  //Set input for unfolding purposes\r
+  SetNumberOfECAClusters(numberOfECAClusters);\r
+  SetRecPoints(recPoints);\r
+  SetDigitsArr(digitsArr);\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::MakeUnfolding()\r
+{\r
+  // Unfolds clusters using the shape of an ElectroMagnetic shower\r
+  // Performs unfolding of all clusters\r
+\r
+  if(fNumberOfECAClusters > 0){\r
+    if (fGeom==0)\r
+      AliFatal("Did not get geometry from EMCALLoader") ;\r
+    //Int_t nModulesToUnfold = fGeom->GetNCells();\r
+\r
+    Int_t numberofNotUnfolded = fNumberOfECAClusters ;\r
+    Int_t index ;\r
+    for(index = 0 ; index < numberofNotUnfolded ; index++){\r
+      AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;\r
+\r
+      //do we really need it?\r
+//      TVector3 gpos;\r
+//      Int_t absId = -1;\r
+//      recPoint->GetGlobalPosition(gpos);\r
+//      fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);\r
+//      if(absId > nModulesToUnfold)\r
+//        break ;\r
+\r
+      Int_t nMultipl = recPoint->GetMultiplicity() ;\r
+      AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;\r
+      Float_t * maxAtEnergy = new Float_t[nMultipl] ;\r
+      Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;\r
+\r
+      if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0\r
+        if(UnfoldClusterV2(recPoint, nMax, maxAt, maxAtEnergy) ){\r
+         fRecPoints->Remove(recPoint);\r
+         fRecPoints->Compress() ;//is it really needed\r
+         index-- ;\r
+         fNumberOfECAClusters-- ;\r
+         numberofNotUnfolded-- ;\r
+       }\r
+      }\r
+      else{\r
+        recPoint->SetNExMax(1) ; //Only one local maximum\r
+      }\r
+\r
+      delete[] maxAt ;\r
+      delete[] maxAtEnergy ;\r
+    }\r
+  }\r
+  // End of Unfolding of clusters\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t AliEMCALUnfolding::UnfoldClusterV2(AliEMCALRecPoint * iniTower, \r
+                                         Int_t nMax, \r
+                                         AliEMCALDigit ** maxAt, \r
+                                         Float_t * maxAtEnergy)\r
+{\r
+  // Extended to whole EMCAL \r
+  // Performs the unfolding of a cluster with nMax overlapping showers \r
+  \r
+  Int_t nPar = 3 * nMax ;\r
+  Float_t * fitparameters = new Float_t[nPar] ;\r
+\r
+  if (fGeom==0)\r
+    AliFatal("Did not get geometry from EMCALLoader") ;\r
+\r
+  Bool_t rv = FindFitV2(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ;\r
+  if( !rv ) {\r
+    // Fit failed, return (and remove cluster? - why? I leave the cluster)\r
+    iniTower->SetNExMax(-1) ;\r
+    delete[] fitparameters ;\r
+    return kFALSE;\r
+  }\r
+\r
+  // create unfolded rec points and fill them with new energy lists\r
+  // First calculate energy deposited in each sell in accordance with\r
+  // fit (without fluctuations): efit[]\r
+  // and later correct this number in acordance with actual energy\r
+  // deposition\r
+\r
+  Int_t nDigits = iniTower->GetMultiplicity() ;\r
+  Float_t * efit = new Float_t[nDigits] ;//new fitted energy in cells\r
+  Float_t xpar=0.,zpar=0.,epar=0.  ;//center of gravity in cell units\r
+\r
+  AliEMCALDigit * digit = 0 ;\r
+  Int_t * digitsList = iniTower->GetDigitsList() ;\r
+  \r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+\r
+  Int_t iparam = 0 ;\r
+  Int_t iDigit = 0 ;\r
+\r
+  for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
+    digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At(digitsList[iDigit] ) ) ;\r
+    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                      iIphi, iIeta,iphi,ieta);\r
+    EvalParsPhiDependence(digit->GetId(),fGeom);\r
+\r
+    efit[iDigit] = 0.;\r
+    iparam = 0;\r
+    while(iparam < nPar ){\r
+      xpar = fitparameters[iparam] ;\r
+      zpar = fitparameters[iparam+1] ;\r
+      epar = fitparameters[iparam+2] ;\r
+      iparam += 3 ;\r
+\r
+      efit[iDigit] += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+    }\r
+\r
+  }\r
+\r
+  // Now create new RecPoints and fill energy lists with efit corrected to fluctuations\r
+  // so that energy deposited in each cell is distributed between new clusters proportionally\r
+  // to its contribution to efit\r
+\r
+  Float_t * energiesList = iniTower->GetEnergiesList() ;\r
+  Float_t ratio = 0 ;\r
+\r
+  iparam = 0 ;\r
+  while(iparam < nPar ){\r
+    xpar = fitparameters[iparam] ;\r
+    zpar = fitparameters[iparam+1] ;\r
+    epar = fitparameters[iparam+2] ;\r
+    iparam += 3 ;\r
+\r
+    AliEMCALRecPoint * recPoint = 0 ;\r
+\r
+    if(fNumberOfECAClusters >= fRecPoints->GetSize())\r
+      fRecPoints->Expand(2*fNumberOfECAClusters) ;\r
+\r
+    //add recpoint\r
+    (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;\r
+    recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;\r
+    fNumberOfECAClusters++ ;\r
+    recPoint->SetNExMax((Int_t)nPar/3) ;\r
+\r
+    Float_t eDigit = 0. ;\r
+    for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){\r
+      digit = dynamic_cast<AliEMCALDigit*>( fDigitsArr->At( digitsList[iDigit] ) ) ;\r
+      fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+      fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                        iIphi, iIeta,iphi,ieta);\r
+      EvalParsPhiDependence(digit->GetId(),fGeom);\r
+      if(efit[iDigit]==0) continue;//just for sure\r
+      ratio = epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) / efit[iDigit] ;\r
+      eDigit = energiesList[iDigit] * ratio ;\r
+      recPoint->AddDigit( *digit, eDigit, kFALSE ) ; //FIXME, need to study the shared case\r
+    }\r
+\r
+  }\r
+\r
+  delete[] fitparameters ;\r
+  delete[] efit ;\r
+\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t AliEMCALUnfolding::FindFitV2(AliEMCALRecPoint * recPoint, AliEMCALDigit ** maxAt, \r
+                                       const Float_t* maxAtEnergy,\r
+                                       Int_t nPar, Float_t * fitparameters) const\r
+{\r
+  // Calls TMinuit to fit the energy distribution of a cluster with several maxima\r
+  // The initial values for fitting procedure are set equal to the\r
+  // positions of local maxima.       \r
+  // Cluster will be fitted as a superposition of nPar/3\r
+  // electromagnetic showers\r
+\r
+  cout<<"inside FindFitV2"<<endl;\r
+\r
+  if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");\r
+       \r
+  if(!gMinuit)\r
+    gMinuit = new TMinuit(100) ;//max 100 parameters\r
+\r
+  gMinuit->mncler();                     // Reset Minuit's list of paramters\r
+  gMinuit->SetPrintLevel(-1) ;           // No Printout\r
+  gMinuit->SetFCN(AliEMCALUnfolding::UnfoldingChiSquareV2) ;\r
+  // To set the address of the minimization function\r
+  TList * toMinuit = new TList();\r
+  toMinuit->AddAt(recPoint,0) ;\r
+  toMinuit->AddAt(fDigitsArr,1) ;\r
+  toMinuit->AddAt(fGeom,2) ;\r
+\r
+  gMinuit->SetObjectFit(toMinuit) ;         // To tranfer pointer to UnfoldingChiSquare\r
+\r
+  // filling initial values for fit parameters\r
+  AliEMCALDigit * digit ;\r
+\r
+  Int_t ierflg  = 0;\r
+  Int_t index   = 0 ;\r
+  Int_t nDigits = (Int_t) nPar / 3 ;\r
+\r
+  Int_t iDigit ;\r
+\r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+\r
+  for(iDigit = 0; iDigit < nDigits; iDigit++){\r
+    digit = maxAt[iDigit];\r
+    fGeom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+    fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                      iIphi, iIeta,iphi,ieta);\r
+\r
+    Float_t energy = maxAtEnergy[iDigit] ;\r
+\r
+    //gMinuit->mnparm(index, "x",  iphi, 0.1, 0, 0, ierflg) ;//original\r
+    gMinuit->mnparm(index, "x",  iphi, 0.05, 0, 0, ierflg) ;\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : x = %d", iphi ) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+    //gMinuit->mnparm(index, "z",  ieta, 0.1, 0, 0, ierflg) ;//original\r
+    gMinuit->mnparm(index, "z",  ieta, 0.05, 0, 0, ierflg) ;\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : z = %d", ieta) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+    //gMinuit->mnparm(index, "Energy",  energy , 0.05*energy, 0., 4.*energy, ierflg) ;//original\r
+    gMinuit->mnparm(index, "Energy",  energy , 0.001*energy, 0., 5.*energy, ierflg) ;//was 0.05\r
+    index++ ;\r
+    if(ierflg != 0){\r
+      Error("FindFit", "EMCAL Unfolding  Unable to set initial value for fit procedure : energy = %f", energy) ;\r
+      toMinuit->Clear();\r
+      delete toMinuit ;\r
+      return kFALSE;\r
+    }\r
+  }\r
+\r
+  Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; \r
+                      // The number of function call slightly depends on it.\r
+  //  Double_t p1 = 1.0 ;// par to gradient \r
+  Double_t p2 = 0.0 ;\r
+  //  Double_t p3 = 3.0 ;\r
+  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls\r
+  //  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient\r
+  gMinuit->SetMaxIterations(5);//was 5\r
+  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings\r
+  //gMinuit->mnexcm("SET PRI", &p3 , 3, ierflg) ;  // printouts\r
+\r
+  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize\r
+  //gMinuit->mnexcm("MINI", &p0, 0, ierflg) ;    // minimize\r
+  if(ierflg == 4){  // Minimum not found\r
+    Error("FindFit", "EMCAL Unfolding  Fit not converged, cluster abandoned " ) ;\r
+    toMinuit->Clear();\r
+    delete toMinuit ;\r
+    return kFALSE ;\r
+  }\r
+  for(index = 0; index < nPar; index++){\r
+    Double_t err = 0. ;\r
+    Double_t val = 0. ;\r
+    gMinuit->GetParameter(index, val, err) ;    // Returns value and error of parameter index\r
+    fitparameters[index] = val ;\r
+  }\r
+\r
+  toMinuit->Clear();\r
+  delete toMinuit ;\r
+  return kTRUE;\r
+\r
+}\r
+\r
+//____________________________________________________________________________\r
+Double_t  AliEMCALUnfolding::ShowerShapeV2(Double_t x, Double_t y)\r
+{ \r
+  // extended to whole EMCAL \r
+  // Shape of the shower\r
+  // If you change this function, change also the gradient evaluation in ChiSquare()\r
+\r
+  Double_t r = fSSPars[7]*TMath::Sqrt(x*x+y*y);\r
+  Double_t rp1  = TMath::Power(r, fSSPars[1]) ;\r
+  Double_t rp5  = TMath::Power(r, fSSPars[5]) ;\r
+  Double_t shape = fSSPars[0]*TMath::Exp( -rp1 * (1. / (fSSPars[2] + fSSPars[3] * rp1) + fSSPars[4] / (1 + fSSPars[6] * rp5) ) ) ;\r
+  return shape ;\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad,\r
+                                                Double_t & fret,\r
+                                                Double_t * x, Int_t iflag)\r
+{\r
+  // Calculates the Chi square for the cluster unfolding minimization\r
+  // Number of parameters, Gradient, Chi squared, parameters, what to do\r
+\r
+  TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;\r
+\r
+  AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint*>( toMinuit->At(0) )  ;\r
+  TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) )  ;\r
+  // A bit buggy way to get an access to the geometry\r
+  // To be revised!\r
+  AliEMCALGeometry *geom = dynamic_cast<AliEMCALGeometry *>(toMinuit->At(2));\r
+\r
+  Int_t * digitsList     = recPoint->GetDigitsList() ;\r
+\r
+  Int_t nOdigits = recPoint->GetDigitsMultiplicity() ;\r
+\r
+  Float_t * energiesList = recPoint->GetEnergiesList() ;\r
+\r
+  fret = 0. ;\r
+  Int_t iparam = 0 ;\r
+\r
+  if(iflag == 2)\r
+    for(iparam = 0 ; iparam < nPar ; iparam++)\r
+      Grad[iparam] = 0 ; // Will evaluate gradient\r
+\r
+  Double_t efit = 0. ;\r
+\r
+  AliEMCALDigit * digit ;\r
+  Int_t iDigit ;\r
+\r
+  Int_t iSupMod =  0 ;\r
+  Int_t iTower  =  0 ;\r
+  Int_t iIphi   =  0 ;\r
+  Int_t iIeta   =  0 ;\r
+  Int_t iphi    =  0 ;//x direction\r
+  Int_t ieta    =  0 ;//z direstion\r
+\r
+\r
+  for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {\r
+    if(energiesList[iDigit]==0) continue;\r
+\r
+    digit = dynamic_cast<AliEMCALDigit*>( digits->At( digitsList[iDigit] ) );\r
+\r
+    geom->GetCellIndex(digit->GetId(),iSupMod,iTower,iIphi,iIeta); \r
+    geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,\r
+                                      iIphi, iIeta,iphi,ieta);\r
+    EvalParsPhiDependence(digit->GetId(),geom);\r
+\r
+    if(iflag == 2){  // calculate gradient\r
+      Int_t iParam = 0 ;\r
+      efit = 0. ;\r
+      while(iParam < nPar ){\r
+        Double_t dx = ((Float_t)iphi - x[iParam]) ;\r
+        iParam++ ;\r
+        Double_t dz = ((Float_t)ieta - x[iParam]) ;\r
+        iParam++ ;\r
+        efit += x[iParam] * ShowerShapeV2(dx,dz) ;\r
+        iParam++ ;\r
+      }\r
+\r
+      Double_t sum = 2. * (efit - energiesList[iDigit]) / energiesList[iDigit] ; // Here we assume, that sigma = sqrt(E)\r
+      iParam = 0 ;\r
+      while(iParam < nPar ){\r
+        Double_t xpar = x[iParam] ;\r
+        Double_t zpar = x[iParam+1] ;\r
+        Double_t epar = x[iParam+2] ;\r
+\r
+        Double_t dr = fSSPars[7]*TMath::Sqrt( ((Float_t)iphi - xpar) * ((Float_t)iphi - xpar) + ((Float_t)ieta - zpar) * ((Float_t)ieta - zpar) );\r
+        Double_t shape = sum * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+       Double_t rp1  = TMath::Power(dr, fSSPars[1]) ;\r
+       Double_t rp5  = TMath::Power(dr, fSSPars[5]) ;\r
+\r
+       Double_t deriv = -2 * TMath::Power(dr,fSSPars[1]-2.) * fSSPars[7] * fSSPars[7] * \r
+         (fSSPars[1] * ( 1/(fSSPars[2]+fSSPars[3]*rp1) + fSSPars[4]/(1+fSSPars[6]*rp5) ) - \r
+          (fSSPars[1]*fSSPars[3]*rp1/( (fSSPars[2]+fSSPars[3]*rp1)*(fSSPars[2]+fSSPars[3]*rp1) ) + \r
+           fSSPars[4]*fSSPars[5]*fSSPars[6]*rp5/( (1+fSSPars[6]*rp5)*(1+fSSPars[6]*rp5) ) ) );\r
+\r
+        //Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) )\r
+       //                                                   - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ;\r
+\r
+        Grad[iParam] += epar * shape * deriv * ((Float_t)iphi - xpar) ;  // Derivative over x\r
+        iParam++ ;\r
+        Grad[iParam] += epar * shape * deriv * ((Float_t)ieta - zpar) ;  // Derivative over z\r
+        iParam++ ;\r
+        Grad[iParam] += shape ;                                  // Derivative over energy\r
+       iParam++ ;\r
+      }\r
+    }\r
+    efit = 0;\r
+    iparam = 0 ;\r
+\r
+    while(iparam < nPar ){\r
+      Double_t xpar = x[iparam] ;\r
+      Double_t zpar = x[iparam+1] ;\r
+      Double_t epar = x[iparam+2] ;\r
+      iparam += 3 ;\r
+      efit += epar * ShowerShapeV2((Float_t)iphi - xpar,(Float_t)ieta - zpar) ;\r
+    }\r
+\r
+    fret += (efit-energiesList[iDigit])*(efit-energiesList[iDigit])/energiesList[iDigit] ;\r
+    // Here we assume, that sigma = sqrt(E) \r
+  }\r
+\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetShowerShapeParams(Double_t *pars){\r
+  for(UInt_t i=0;i<7;++i)\r
+    fSSPars[i]=pars[i];\r
+  if(pars[2]==0. && pars[3]==0.) fSSPars[2]=1.;//to avoid dividing by 0\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetPar5(Double_t *pars){\r
+  for(UInt_t i=0;i<3;++i)\r
+    fPar5[i]=pars[i];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::SetPar6(Double_t *pars){\r
+  for(UInt_t i=0;i<3;++i)\r
+    fPar6[i]=pars[i];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalPar5(Double_t phi){\r
+  //\r
+  //Evaluate the 5th parameter of the shower shape function\r
+  //phi in degrees range (-10,10)\r
+  //\r
+  //fSSPars[5] = 12.31 - phi*0.007381 - phi*phi*0.06936;\r
+  fSSPars[5] = fPar5[0] + phi * fPar5[1] + phi*phi * fPar5[2];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalPar6(Double_t phi){\r
+  //\r
+  //Evaluate the 6th parameter of the shower shape function\r
+  //phi in degrees range (-10,10)\r
+  //\r
+  //fSSPars[6] = 0.05452 + phi*0.0001228 + phi*phi*0.001361;\r
+  fSSPars[6] = fPar6[0] + phi * fPar6[1] + phi*phi * fPar6[2];\r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliEMCALUnfolding::EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom){\r
+  //\r
+  // calculate params p5 and p6 depending on the phi angle in global coordinate\r
+  // for the cell with given absId index\r
+  //\r
+  Double_t etaGlob = 0.;//eta in global c.s. - unused\r
+  Double_t phiGlob = 0.;//phi in global c.s. in radians\r
+  geom->EtaPhiFromIndex(absId, etaGlob, phiGlob);\r
+  phiGlob*=180./TMath::Pi();\r
+  phiGlob-=90.;\r
+  phiGlob-= (Double_t)((Int_t)geom->GetSuperModuleNumber(absId)/2 * 20);\r
+\r
+  EvalPar5(phiGlob);\r
+  EvalPar6(phiGlob);\r
+}\r
+\r
diff --git a/EMCAL/AliEMCALUnfolding.h b/EMCAL/AliEMCALUnfolding.h
new file mode 100644 (file)
index 0000000..aada310
--- /dev/null
@@ -0,0 +1,91 @@
+#ifndef ALIEMCALUNFOLDING_H\r
+#define ALIEMCALUNFOLDING_H\r
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice                               */\r
+     \r
+//_________________________________________________________________________\r
+//  Base class for the cluster unfolding algorithm \r
+//*-- Author: Adam Matyja (SUBATECH)\r
+\r
+// --- ROOT system ---\r
+#include "AliLog.h"\r
+#include "TObject.h" \r
+//class TTree;\r
+\r
+// --- Standard library ---\r
+\r
+// --- AliRoot header files ---\r
+class AliEMCALGeometry ;\r
+//class AliEMCALCalibData ;\r
+//class AliCaloCalibPedestal ;\r
+class AliEMCALRecPoint ; \r
+class AliEMCALDigit ;\r
+\r
+\r
+class AliEMCALUnfolding : public TObject {\r
+\r
+public:\r
+\r
+  AliEMCALUnfolding() ;        // default ctor\r
+  virtual ~AliEMCALUnfolding() ; // dtorEM\r
+  AliEMCALUnfolding(AliEMCALGeometry* geometry);// constructor\r
+  AliEMCALUnfolding(AliEMCALGeometry* geometry,Float_t ECALocMaxCut,Double_t *SSPars,Double_t *Par5,Double_t *Par6);// constructor\r
+\r
+  virtual void Init() ;\r
+  virtual void SetInput(Int_t numberOfECAClusters,TObjArray *recPoints,TClonesArray *digitsArr);\r
+\r
+  //setters and getters\r
+  virtual void SetNumberOfECAClusters(Int_t n) { fNumberOfECAClusters = n; }\r
+  virtual Int_t GetNumberOfECAClusters() const { return fNumberOfECAClusters; }\r
+  virtual void SetRecPoints(TObjArray *rec) { fRecPoints = rec; }\r
+  virtual TObjArray * GetRecPoints() const { return fRecPoints; }\r
+  virtual void SetDigitsArr(TClonesArray *digit) { fDigitsArr = digit; }\r
+  virtual TClonesArray * GetDigitsArr() const { return fDigitsArr; }\r
+  virtual void SetECALocalMaxCut(Float_t cut) { fECALocMaxCut = cut ; }\r
+  virtual Float_t GetECALocalMaxCut() const { return fECALocMaxCut; }\r
+\r
+  //unfolding main methods\r
+  virtual void   MakeUnfolding();\r
+  static Double_t ShowerShapeV2(Double_t x, Double_t y) ; // Shape of EM shower used in unfolding; \r
+                                              //class member function (not object member function)\r
+  static void UnfoldingChiSquareV2(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)  ;\r
+                                            // Chi^2 of the fit. Should be static to be passes to MINUIT\r
+  virtual void SetShowerShapeParams(Double_t *pars) ;\r
+  virtual Double_t* GetShowerShapeParams() const { return fSSPars ; }\r
+  virtual void SetPar5(Double_t *pars) ;\r
+  virtual Double_t* GetPar5() const { return fPar5 ; }\r
+  virtual void SetPar6(Double_t *pars) ;\r
+  virtual Double_t* GetPar6() const { return fPar6 ; }\r
+\r
+protected:\r
+  Int_t   fNumberOfECAClusters ;     // number of clusters found in EC section\r
+  Float_t fECALocMaxCut ;            // minimum energy difference to distinguish local maxima in a cluster\r
+  AliEMCALGeometry     * fGeom;       //! pointer to geometry for utilities\r
+  TObjArray    *fRecPoints; // Array with EMCAL clusters\r
+  TClonesArray *fDigitsArr; // Array with EMCAL digits\r
+\r
+private:\r
+  AliEMCALUnfolding(const AliEMCALUnfolding &); //copy ctor\r
+  AliEMCALUnfolding & operator = (const AliEMCALUnfolding &);\r
+  \r
+  Bool_t         UnfoldClusterV2(AliEMCALRecPoint * iniEmc, Int_t Nmax, \r
+                                AliEMCALDigit ** maxAt,\r
+                                Float_t * maxAtEnergy ); //Unfolds cluster using TMinuit package\r
+  Bool_t  FindFitV2(AliEMCALRecPoint * emcRP, AliEMCALDigit ** MaxAt, const Float_t * maxAtEnergy, \r
+                   Int_t NPar, Float_t * FitParametres) const; //Used in UnfoldClusters, calls TMinuit\r
+\r
+  static Double_t fSSPars[8];//! Unfolding shower shape parameters\r
+  // function:\r
+  // f(r)=exp(-(p0*r)^p1 * (1/(p2+p3*(p0*r)^p1)+p4/(1+p6*(p0*r)^p5) ) )\r
+  // p0,p1,p2,p3,p4 are fixed\r
+  // params p5 and p6 are phi-dependent and set in ShowerShapeV2\r
+  static Double_t fPar5[3];//! UF SSPar nr 5 = p0 + phi*p1 + phi^2 *p2\r
+  static Double_t fPar6[3];//! UF SSPar nr 6 = p0 + phi*p1 + phi^2 *p2\r
+  static void EvalPar5(Double_t phi);\r
+  static void EvalPar6(Double_t phi);\r
+  static void EvalParsPhiDependence(Int_t absId, AliEMCALGeometry *geom);\r
+\r
+  ClassDef(AliEMCALUnfolding,1)  // Unfolding algorithm class \r
+} ;\r
+\r
+#endif // AliEMCALUNFOLDING_H\r
index 2053121..73c9d6e 100644 (file)
@@ -7,6 +7,7 @@
 #pragma link C++ class AliEMCALClusterizer+;
 #pragma link C++ class AliEMCALClusterizerv1+;
 #pragma link C++ class AliEMCALClusterizerNxN+;
+#pragma link C++ class AliEMCALUnfolding+;
 #pragma link C++ class AliEMCALTrack+;
 #pragma link C++ class AliEMCALTracker+;
 #pragma link C++ class AliEMCALPID+;
index ae76b12..c1b2ec8 100644 (file)
@@ -6,6 +6,7 @@ AliEMCALReconstructor.cxx \
 AliEMCALClusterizer.cxx \
 AliEMCALClusterizerv1.cxx \
 AliEMCALClusterizerNxN.cxx \
+AliEMCALUnfolding.cxx \
 AliEMCALTrack.cxx \
 AliEMCALTracker.cxx \
 AliEMCALPID.cxx \