#include "TClonesArray.h"
#include "TTree.h"
#include "TGeoManager.h"
+#include "TROOT.h"
+#include "TInterpreter.h"
// --- AliRoot Analysis Steering
#include "AliAnalysisTask.h"
#include "AliCDBEntry.h"
#include "AliLog.h"
#include "AliVEventHandler.h"
+#include "AliAODInputHandler.h"
// --- EMCAL
#include "AliEMCALRecParam.h"
#include "AliEMCALDigit.h"
#include "AliCaloCalibPedestal.h"
#include "AliEMCALCalibData.h"
+#include "AliEMCALRecoUtils.h"
#include "AliAnalysisTaskEMCALClusterize.h"
, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
, fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
, fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kTRUE)
- , fRun(-1)
+ , fRun(-1), fRecoUtils(0), fConfigName(""), fCellLabels()
{
//ctor
- for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
+ for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
+ for(Int_t j = 0; j < 12672; j++) fCellLabels[j] = -1;
fDigitsArr = new TClonesArray("AliEMCALDigit",200);
fClusterArr = new TObjArray(100);
fCaloClusterArr = new TObjArray(100);
fRecParam = new AliEMCALRecParam;
- fBranchNames="ESD:AliESDHeader.,EMCALCells.";
+ fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
+ fRecoUtils = new AliEMCALRecoUtils();
}
//________________________________________________________________________
, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
, fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
, fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kFALSE)
- , fRun(-1)
+ , fRun(-1), fRecoUtils(0), fConfigName(""), fCellLabels()
{
// Constructor
- for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
+ for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
+ for(Int_t j = 0; j < 24*48*11; j++) fCellLabels[j] = -1;
fDigitsArr = new TClonesArray("AliEMCALDigit",200);
fClusterArr = new TObjArray(100);
fCaloClusterArr = new TObjArray(100);
fRecParam = new AliEMCALRecParam;
- fBranchNames="ESD:AliESDHeader.,EMCALCells.";
+ fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
+ fRecoUtils = new AliEMCALRecoUtils();
}
+
//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
{
delete fCaloClusterArr;
}
- if(fClusterizer) {delete fClusterizer;}
- if(fUnfolder) {delete fUnfolder; }
+ if(fClusterizer) delete fClusterizer;
+ if(fUnfolder) delete fUnfolder;
+ if(fRecoUtils) delete fRecoUtils;
+
}
+//-------------------------------------------------------------------
+void AliAnalysisTaskEMCALClusterize::Init()
+{
+ //Init analysis with configuration macro if available
+
+ if(gROOT->LoadMacro(fConfigName) >=0){
+ printf("Configure analysis with %s\n",fConfigName.Data());
+ AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
+ fGeomName = clus->fGeomName;
+ fLoadGeomMatrices = clus->fLoadGeomMatrices;
+ fOCDBpath = clus->fOCDBpath;
+ fRecParam = clus->fRecParam;
+ fJustUnfold = clus->fJustUnfold;
+ fFillAODFile = clus->fFillAODFile;
+ fRecoUtils = clus->fRecoUtils;
+ fConfigName = clus->fConfigName;
+ fOutputAODBranchName = clus->fOutputAODBranchName;
+ for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
+
+ }
+
+}
+
//-------------------------------------------------------------------
void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
{
// Init geometry, create list of output clusters
-
+
fGeom = AliEMCALGeometry::GetInstance(fGeomName) ;
fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
fOutputAODBranch->SetName(fOutputAODBranchName);
AddAODBranch("TClonesArray", &fOutputAODBranch);
- Info("UserCreateOutputObjects","Create Branch: %s",fOutputAODBranchName.Data());
+
}
//________________________________________________________________________
// Main loop
// Called for each event
+ //printf("--- Event %d -- \n",(Int_t)Entry());
//Remove the contents of output list set in the previous event
fOutputAODBranch->Clear("C");
- AliVEvent * event = InputEvent();
+ //Get the event
+ AliVEvent * event = 0;
+ AliESDEvent * esdevent = 0;
+
+ //Check if input event are embedded events
+ //If so, take output event
+ AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+ if (aodIH && aodIH->GetMergeEvents()) {
+ //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
+ event = AODEvent();
+// printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
+// InputEvent()->GetEMCALCells()->GetNumberOfCells());
+// printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
+// AODEvent()->GetEMCALCells()->GetNumberOfCells());
+ }
+ else {
+ event = InputEvent();
+ esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
+ }
+
if (!event) {
Error("UserExec","Event not available");
return;
fGeomMatrixSet=kTRUE;
}//ESD
}//Load matrices from Data
+
+ //Recover time dependent corrections, put then in recalibration histograms. Do it once
+ fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
+
}//first event
//Get the list of cells needed for unfolding and reclustering
AliVCluster *clus = event->GetCaloCluster(i);
if(clus->IsEMCAL()){
//printf("Org Cluster %d, E %f\n",i, clus->E());
+
+ //recalibrate/remove bad channels/etc if requested
+ if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
+ //printf("Remove cluster\n");
+ continue;
+ }
+
+ if(fRecoUtils->IsRecalibrationOn()){
+ //printf("Energy before %f ",clus->E());
+ fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
+ //printf("; after %f\n",clus->E());
+ }
+ //Cast to ESD or AOD, needed to create the cluster array
AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
if (esdCluster){
//CLEAN-UP
fUnfolder->Clear();
+
//printf("nClustersOrg %d\n",nClustersOrg);
}
//-------------------------------------------
//-------------------------------------------------------------------------------------
//Transform CaloCells into Digits
//-------------------------------------------------------------------------------------
+
+ //In case of MC, first loop on the clusters and fill MC label to array
+ //.....................................................................
+
+// for(Int_t j = 0; j < 24*48*11; j++) {
+// if(fCellLabels[j]!=-1) printf("Not well initialized cell %d, label %d\n",j,fCellLabels[j]) ;
+// }
+
+ for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
+ {
+ AliVCluster *clus = event->GetCaloCluster(i);
+ if(clus->IsEMCAL()){
+
+ Int_t label = clus->GetLabel();
+ UShort_t * index = clus->GetCellsAbsId() ;
+ for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
+ fCellLabels[index[icell]]=label;
+ //printf("1) absID %d, label %d\n",index[icell], fCellLabels[index[icell]]);
+ }
+ nClustersOrg++;
+ }
+ }
+
+ // Create digits
+ //.................
Int_t idigit = 0;
Int_t id = -1;
Float_t amp = -1;
TTree *digitsTree = new TTree("digitstree","digitstree");
digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
+
for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
{
if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
amp = cellAmplitude;
id = cellNumber;
- AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
- digit->SetId(id);
- digit->SetAmplitude(amp);
- digit->SetTime(time);
- digit->SetTimeR(time);
- digit->SetIndexInList(idigit);
- digit->SetType(AliEMCALDigit::kHG);
- //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
+ Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
+ fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
+ fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
+
+ //Do not include bad channels found in analysis?
+ if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
+ fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
+ fCellLabels[id]=-1; //reset the entry in the array for next event
+ //printf("Remove channel %d\n",id);
+ continue;
+ }
+
+ //Recalibrate?
+ if(fRecoUtils->IsRecalibrationOn()){
+ //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
+ amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
+ }
+
+ //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
+ new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
+ //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
+ //else printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
+ fCellLabels[id]=-1; //reset the entry in the array for next event
+
+ //AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
+ //digit->SetId(id);
+ //digit->SetAmplitude(amp);
+ //digit->SetTime(time);
+ //digit->SetTimeR(time);
+ //digit->SetIndexInList(idigit);
+ //digit->SetType(AliEMCALDigit::kHG);
+
+ //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
idigit++;
}
fClusterizer->SetInput(digitsTree);
fClusterizer->SetOutput(clustersTree);
fClusterizer->Digits2Clusters("");
-
+
//-------------------------------------------------------------------------------------
//Transform the recpoints into AliVClusters
//-------------------------------------------------------------------------------------
TBranch *branch = clustersTree->GetBranch("EMCALECARP");
branch->SetAddress(&fClusterArr);
branch->GetEntry(0);
-
+
RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
//---CLEAN UP-----
digitsTree ->Delete("all");
}
+ //Recalculate track-matching for the new clusters, only with ESDs
+ if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
+
+
//-------------------------------------------------------------------------------------
//Put the new clusters in the AOD list
//-------------------------------------------------------------------------------------
Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntries();
- //Info("UserExec","New clusters %d",kNumberOfCaloClusters);
- //if(nClustersOrg!=kNumberOfCaloClusters) Info("UserExec","Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
+ //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
//if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
+
+ //Add matched track, if any, only with ESDs
+ if(esdevent){
+ Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
+ if(trackIndex >= 0){
+ newCluster->AddTrackMatched(event->GetTrack(trackIndex));
+ if(DebugLevel() > 1)
+ Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
+ }
+ }
+
+ //In case of new bad channels, recalculate distance to bad channels
+ if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
+ //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
+ fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
+ //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
+ }
+
+// if(newCluster->GetNLabels()>0){
+// printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
+// UShort_t * newindex = newCluster->GetCellsAbsId() ;
+// for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
+// printf("\t absID %d\n",newindex[icell]);
+// }
+// }
+//
+// printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
+
+ newCluster->SetID(i);
new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
}
AliVEvent * event = InputEvent();
if (!event)
{
- Warning("AccessODCD","Event not available!!!");
+ Warning("AccessOCDB","Event not available!!!");
return kFALSE;
}
Info("AccessODCD()"," Begin");
fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+
+
AliCDBManager *cdb = AliCDBManager::Instance();
- if (fOCDBpath.Length())
+
+
+ if (fOCDBpath.Length()){
cdb->SetDefaultStorage(fOCDBpath.Data());
+ Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
+ }
+
cdb->SetRun(event->GetRunNumber());
+
//
// EMCAL from RAW OCDB
if (fOCDBpath.Contains("alien:"))
cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
}
+
TString path = cdb->GetDefaultStorage()->GetBaseFolder();
-
+
// init parameters:
+
//Get calibration parameters
if(!fCalibData)
{
AliCDBEntry *entry = (AliCDBEntry*)
AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
+
if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
}
{
AliCDBEntry *entry = (AliCDBEntry*)
AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
+
if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
}
AliFatal("Dead map not found in CDB!");
InitClusterization();
+
return kTRUE;
}
fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
fClusterizer->SetInputCalibrated ( kTRUE );
-
+ fClusterizer->SetJustClusters ( kTRUE );
//In case of unfolding after clusterization is requested, set the corresponding parameters
if(fRecParam->GetUnfold()){
Int_t i=0;
clus->SetM02(elipAxis[0]*elipAxis[0]) ;
clus->SetM20(elipAxis[1]*elipAxis[1]) ;
clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
+
+ //MC
+ Int_t parentMult = 0;
+ Int_t *parentList = recPoint->GetParents(parentMult);
+ clus->SetLabel(parentList, parentMult);
+
clusArray->Add(clus);
} // recPoints loop
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff