track number corrected
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALReconstructor.cxx
index a5a2d6d..759a763 100644 (file)
@@ -27,6 +27,7 @@
 // --- ROOT system ---
 #include <TList.h>
 #include <TClonesArray.h>
+#include <TH2.h>
 
 // --- Standard library ---
 
@@ -48,7 +49,6 @@
 #include "AliRawReader.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
-#include "AliEMCALRecParam.h"
 #include "AliEMCALGeometry.h"
 #include "AliEMCAL.h"
 #include "AliEMCALHistoUtilities.h"
 
 ClassImp(AliEMCALReconstructor) 
 
-AliEMCALRecParam* AliEMCALReconstructor::fgkRecParam = 0;  // EMCAL rec. parameters
+const AliEMCALRecParam* AliEMCALReconstructor::fgkRecParam = 0;  // EMCAL rec. parameters
 AliEMCALRawUtils* AliEMCALReconstructor::fgRawUtils = 0;   // EMCAL raw utilities class
+AliEMCALClusterizer* AliEMCALReconstructor::fgClusterizer = 0;   // EMCAL clusterizer class
 TClonesArray*     AliEMCALReconstructor::fgDigitsArr = 0;  // shoud read just once at event
-
 //____________________________________________________________________________
 AliEMCALReconstructor::AliEMCALReconstructor() 
   : fDebug(kFALSE), fList(0), fGeom(0) 
 {
   // ctor
-  InitRecParam();
 
   fgRawUtils = new AliEMCALRawUtils;
+  fgClusterizer = new AliEMCALClusterizerv1;
 
   //To make sure we match with the geometry in a simulation file,
   //let's try to get it first.  If not, take the default geometry
@@ -112,24 +112,6 @@ void AliEMCALReconstructor::Init()
 }
 
 //____________________________________________________________________________
-void AliEMCALReconstructor::InitRecParam() const
-{
-  // Check if the instance of AliEMCALRecParam exists, 
-  // if not, get it from OCDB if available, otherwise create a default one
-
- if (!fgkRecParam  && (AliCDBManager::Instance()->IsDefaultStorageSet())) {
-    AliCDBEntry *entry = (AliCDBEntry*) 
-      AliCDBManager::Instance()->Get("EMCAL/Config/RecParam");
-    if (entry) fgkRecParam =  (AliEMCALRecParam*) entry->GetObject();
-  }
-  
-  if(!fgkRecParam){
-    AliWarning("The Reconstruction parameters for EMCAL nonitialized - Used default one");
-    fgkRecParam = new AliEMCALRecParam;
-  }
-}
-
-//____________________________________________________________________________
 void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
 {
   // method called by AliReconstruction; 
@@ -141,18 +123,22 @@ void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree)
   AliCodeTimerAuto("")
 
   ReadDigitsArrayFromTree(digitsTree);
+  fgClusterizer->InitParameters();
+  fgClusterizer->SetOutput(clustersTree);
 
   if(fgDigitsArr && fgDigitsArr->GetEntries()) {
 
-    AliEMCALClusterizerv1 clu;
-    clu.SetInput(digitsTree);
-    clu.SetOutput(clustersTree);
+    fgClusterizer->SetInput(digitsTree);
+    
     if(Debug())
-      clu.Digits2Clusters("deb all") ;
+      fgClusterizer->Digits2Clusters("deb all") ;
     else
-      clu.Digits2Clusters("") ;
+      fgClusterizer->Digits2Clusters("");
+    
+    fgClusterizer->Clear();
 
   }
+
 }
 
 //____________________________________________________________________________
@@ -165,18 +151,18 @@ void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digits
 
   rawReader->Reset() ; 
 
-  TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",100);
+  TClonesArray *digitsArr = new TClonesArray("AliEMCALDigit",200);
   Int_t bufsize = 32000;
   digitsTree->Branch("EMCAL", &digitsArr, bufsize);
 
   //must be done here because, in constructor, option is not yet known
   fgRawUtils->SetOption(GetOption());
 
-  fgRawUtils->SetRawFormatHighLowGainFactor(fgkRecParam->GetHighLowGainFactor());
-  fgRawUtils->SetRawFormatOrder(fgkRecParam->GetOrderParameter());
-  fgRawUtils->SetRawFormatTau(fgkRecParam->GetTau());
-  fgRawUtils->SetNoiseThreshold(fgkRecParam->GetNoiseThreshold());
-  fgRawUtils->SetNPedSamples(fgkRecParam->GetNPedSamples());
+  fgRawUtils->SetRawFormatHighLowGainFactor(GetRecParam()->GetHighLowGainFactor());
+  fgRawUtils->SetRawFormatOrder(GetRecParam()->GetOrderParameter());
+  fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
+  fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
+  fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
 
   fgRawUtils->Raw2Digits(rawReader,digitsArr);
 
@@ -196,7 +182,6 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   // Works on the current event
   //  printf(" ## AliEMCALReconstructor::FillESD() is started ### \n ");
   //return;
-  const double timeScale = 1.e+11; // transition constant from sec to 0.01 ns 
 
   //######################################################
   //#########Calculate trigger and set trigger info###########
@@ -258,8 +243,24 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   //printf("\n triggerAmplitudes ");
   //for(int i=0; i<4; i++) printf(" %i %f : ", i, triggerAmplitudes[i]);
   //printf("\n");
-  tr.Print("");
-
+  //tr.Print("");
+  //
+  // Trigger jet staff
+  //
+  if(tr.GetNJetThreshold()>0) {
+    // Jet phi/eta
+    Int_t n0 = triggerPosition.GetSize();
+    const TH2F *hpatch = tr.GetJetMatrixE();
+    triggerPosition.Set(n0 + 2);
+    for(Int_t i=0; i<2; i++) triggerPosition[n0+i] = hpatch->GetMean(i+1);   
+    // Add jet ampitudes
+    n0 = triggerAmplitudes.GetSize();
+    triggerAmplitudes.Set(n0 + tr.GetNJetThreshold());
+    Double_t *ampJet = tr.GetL1JetThresholds();
+    for(Int_t i=0; i<tr.GetNJetThreshold(); i++){
+      triggerAmplitudes[n0 + i] = Float_t(ampJet[i]);
+    }
+  }
   esd->AddEMCALTriggerPosition(triggerPosition);
   esd->AddEMCALTriggerAmplitudes(triggerAmplitudes);
   // Fill trigger hists
@@ -347,22 +348,24 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
     Float_t *amplFloat = clust->GetEnergiesList();
     Int_t   *digitInts = clust->GetAbsId();
     TArrayS absIdList(cellMult);
-    //Uncomment when unfolding is done
-    //TArrayD fracList(cellMult);
+    TArrayD fracList(cellMult);
 
     Int_t newCellMult = 0;
     for (Int_t iCell=0; iCell<cellMult; iCell++) {
       if (amplFloat[iCell] > 0) {
       absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
-    //Uncomment when unfolding is done
-    //fracList[newCellMult] = amplFloat[iCell]/emcCells.GetCellAmplitude(digitInts[iCell]);
+      //Uncomment when unfolding is done
+      //if(emcCells.GetCellAmplitude(digitInts[iCell])>0)
+      //fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell])*calibration);//get cell calibration value 
+      //else
+      fracList[newCellMult] = 0; 
       newCellMult++;
       }
     }
-    absIdList.Set(newCellMult);
-    //Uncomment when unfolding is done
-    //fracList.Set(newCellMult);
 
+    absIdList.Set(newCellMult);
+    fracList.Set(newCellMult);
+    
     if(newCellMult > 0) { // accept cluster if it has some digit
       nClustersNew++;
       //Primaries
@@ -376,22 +379,19 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       ec->SetNCells(newCellMult);
       //Change type of list from short to ushort
       UShort_t *newAbsIdList  = new UShort_t[newCellMult];
-      //Uncomment when unfolding is done
-      //Double_t *newFracList  = new Double_t[newCellMult];
+      Double_t *newFracList  = new Double_t[newCellMult];
       for(Int_t i = 0; i < newCellMult ; i++) {
         newAbsIdList[i]=absIdList[i];
-      //Uncomment when unfolding is done
-      //newFracList[i]=fracList[i];
+        newFracList[i]=fracList[i];
       }
       ec->SetCellsAbsId(newAbsIdList);
-      //Uncomment when unfolding is done
-      //ec->SetCellsAmplitudeFraction(newFracList);
+      ec->SetCellsAmplitudeFraction(newFracList);
       ec->SetClusterDisp(clust->GetDispersion());
       ec->SetClusterChi2(-1); //not yet implemented
       ec->SetM02(elipAxis[0]*elipAxis[0]) ;
       ec->SetM20(elipAxis[1]*elipAxis[1]) ;
-      ec->SetM11(-1) ;        //not yet implemented
-
+      ec->SetTOF(clust->GetTime()) ; //time-of-fligh
+      ec->SetNExMax(clust->GetNExMax());          //number of local maxima
       TArrayI arrayTrackMatched(1);// Only one track, temporal solution.
       arrayTrackMatched[0]= matchedTrack[iClust];
       ec->AddTracksMatched(arrayTrackMatched);
@@ -400,10 +400,10 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
       ec->AddLabels(arrayParents);
 
       // add the cluster to the esd object
-     esd->AddCaloCluster(ec);
+      esd->AddCaloCluster(ec);
       delete ec;
-     //delete [] newAbsIdList ;
-     //delete [] newFracList ;
+      delete [] newAbsIdList ;
+      delete [] newFracList ;
    }
  } // cycle on clusters
 
@@ -419,13 +419,14 @@ void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
   pid->SetReconstructor(kTRUE);
   pid->RunPID(esd);
   delete pid;
-
+  
   delete digits;
   delete clusters;
-
+  
   // printf(" ## AliEMCALReconstructor::FillESD() is ended : ncl %i -> %i ### \n ",nClusters, nClustersNew); 
 }
 
+//__________________________________________________________________________
 void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
 {
   // See AliEMCALClusterizer::SetInput(TTree *digitsTree);