//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
: AliAnalysisTaskSE(name)
- , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
+ , fGeom(0), fGeomName("EMCAL_COMPLETEV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
, fCalibData(0), fPedestalData(0), fOCDBpath("raw://"), fAccessOCDB(kFALSE)
, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
, fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
, fRun(-1), fRecoUtils(0), fConfigName("")
, fCellLabels(), fCellSecondLabels(), fMaxEvent(1000000000), fDoTrackMatching(kFALSE)
- {
+{
//ctor
for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
- for(Int_t j = 0; j < 24*48*11; j++) {
- fCellLabels[j] = -1;
- fCellSecondLabels[j] = -1;
- }
+ for(Int_t j = 0; j < 24*48*11; j++) {
+ fCellLabels[j] = -1;
+ fCellSecondLabels[j] = -1;
+ }
fDigitsArr = new TClonesArray("AliEMCALDigit",200);
fClusterArr = new TObjArray(100);
- fCaloClusterArr = new TObjArray(100);
+ fCaloClusterArr = new TObjArray(1000);
fRecParam = new AliEMCALRecParam;
fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
fRecoUtils = new AliEMCALRecoUtils();
//________________________________________________________________________
AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
- , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
+ , fGeom(0), fGeomName("EMCAL_COMPLETEV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
, fCalibData(0), fPedestalData(0), fOCDBpath("raw://"), fAccessOCDB(kFALSE)
, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
, fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
fCaloClusterArr->Delete();
delete fCaloClusterArr;
}
-
+
if(fClusterizer) delete fClusterizer;
if(fUnfolder) delete fUnfolder;
if(fRecoUtils) delete fRecoUtils;
-
+
}
//-------------------------------------------------------------------
fDoTrackMatching = clus->fDoTrackMatching;
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) ;
if(fOutputAODBranchName.Length()!=0){
fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
- //printf("GOOD channel\n");
}
else {
aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
- //printf("BAD channel\n");
}
}
aodEMcells.Sort();
AliVEvent* event = InputEvent();
AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
-
+
Double_t pos[3] ;
Double_t covVtx[6];
for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
{
// Main loop
// Called for each event
-
+
AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
Int_t eventN = Entry();
if(aodIH) eventN = aodIH->GetReadEntry();
if (eventN > fMaxEvent) return;
//printf("Clusterizer --- Event %d-- \n",eventN);
-
+
//Remove the contents of output list set in the previous event
fOutputAODBranch->Clear("C");
-// if(fOutputAODBranchName.Length()!=0){
-// //Remove the contents of output list set in the previous event
-// fOutputAODBranch->Clear("C");
-// }
-// else{
-// //Reset and put new clusters in standard branch, also cells
-// AODEvent()->ResetStd(0, 1, 0, 0, 0, AODEvent()->GetNumberOfCaloClusters(), 0, 0);
-// fOutputAODBranch = AODEvent()->GetCaloClusters();// Put new clusters in standard branch
-// }
-
//Magic line to write events to AOD file
AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
LoadBranches();
-
//Init pointers, clusterizer, ocdb
if(fAccessOCDB) AccessOCDB();
InitClusterization();
-
+
//Get the event
AliVEvent * event = 0;
AliESDEvent * esdevent = 0;
-
- //Fill output event with headerΓΈ
-
+
+ //Fill output event with header
+
//Check if input event are embedded events
//If so, take output event
if (aodIH && aodIH->GetMergeEvents()) {
if(!aodIH->GetMergeEMCALCells())
AliFatal("Events merged but not EMCAL cells, check analysis settings!");
-// printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
-// InputEvent()->GetEMCALCells()->GetNumberOfCells());
-// printf("MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
-// aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
-// for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
-// AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
-// if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E %f\n",icl,sigCluster->E());
-// }
-// printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
-// AODEvent()->GetEMCALCells()->GetNumberOfCells());
+ // printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
+ // InputEvent()->GetEMCALCells()->GetNumberOfCells());
+ // printf("MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
+ // aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
+ // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
+ // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
+ // if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E %f\n",icl,sigCluster->E());
+ // }
+ // printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
+ // AODEvent()->GetEMCALCells()->GetNumberOfCells());
}
else {
Error("UserExec","Event not available");
return;
}
-
+
//-------------------------------------------------------------------------------------
//Set the geometry matrix, for the first event, skip the rest
//-------------------------------------------------------------------------------------
//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
- AliVCaloCells *cells= event->GetEMCALCells();
- Int_t nClustersOrg = 0;
+ AliVCaloCells *cells = event->GetEMCALCells();
+ Int_t nClustersOrg = 0;
+ Double_t cellAmplitude = 0;
+ Double_t cellTime = 0;
+ Short_t cellNumber = 0;
+
//-------------------------------------------
//---------Unfolding clusters----------------
//-------------------------------------------
{
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()){
+
+ //Calibrate cluster
//printf("Energy before %f ",clus->E());
fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
//printf("; after %f\n",clus->E());
- }
+
+ //CalibrateCells
+ for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
+ {
+ if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
+ break;
+
+
+ Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
+ fGeom->GetCellIndex(cellNumber,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[cellNumber] =-1; //reset the entry in the array for next event
+ fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
+ continue;
+ }
+
+ cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
+
+ }// cells loop
+ }// recalibrate
+
//Cast to ESD or AOD, needed to create the cluster array
AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
//CLEAN-UP
fUnfolder->Clear();
-
- //printf("nClustersOrg %d\n",nClustersOrg);
- }
+ } // just unfold ESD/AOD cluster
+
//-------------------------------------------
//---------- Recluster cells ----------------
//-------------------------------------------
//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, label1 %d\n",j,fCellLabels[j] ) ;
-// if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
-// }
+ // for(Int_t j = 0; j < 24*48*11; j++) {
+ // if(fCellLabels[j] !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j] ) ;
+ // if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
+ // }
Int_t nClusters = event->GetNumberOfCaloClusters();
if(aodIH && aodIH->GetEventToMerge()) //Embedding
clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
else
clus = event->GetCaloCluster(i);
-
+
if(!clus) {
printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
return;
nClustersOrg++;
}
}
-
+
// Create digits
//.................
Int_t idigit = 0;
Float_t amp = -1;
Float_t time = -1;
- Double_t cellAmplitude = 0;
- Double_t cellTime = 0;
- Short_t cellNumber = 0;
-
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)
//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
- fCellSecondLabels[id]=-1; //reset the entry in the array for next event
+ fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
+ fCellLabels[id] =-1; //reset the entry in the array for next event
+ fCellSecondLabels[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++;
//Do the clusterization
//-------------------------------------------------------------------------------------
TTree *clustersTree = new TTree("clustertree","clustertree");
-
+
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);
+ if(!fCaloClusterArr){
+ printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
+ return;
+ }
+
+ if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
+ printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters: Input entries %d - Output entries %d - %d (not fast)\n",
+ fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
+ }
+
//Reset the array with second labels for this event
memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
fClusterizer->Clear();
fDigitsArr ->Clear("C");
fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
-
+
clustersTree->Delete("all");
digitsTree ->Delete("all");
}
//Recalculate track-matching for the new clusters, only with ESDs
if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
-
+
//-------------------------------------------------------------------------------------
//Put the new clusters in the AOD list
//-------------------------------------------------------------------------------------
- Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntries();
+ Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
//printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
//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]);
-// }
-// }
+ // 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());
-// printf("labels: ");
-// for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
-// printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
-// printf("\n ");
+ // printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
+ // printf("labels: ");
+ // for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
+ // printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
+ // printf("\n ");
newCluster->SetID(i);
new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
delete fUnfolder;
fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
return;
- }
-
+ }
+
//First init the clusterizer
delete fClusterizer;
if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerv2) { //FIX this other way.
- AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
- clusterizer->SetNRowDiff(2);
- clusterizer->SetNColDiff(2);
- fClusterizer = clusterizer;
+ AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
+ clusterizer->SetNRowDiff(2);
+ clusterizer->SetNColDiff(2);
+ fClusterizer = clusterizer;
} else{
AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
}
-
+
//Now set the parameters
fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
{
// Restore clusters from recPoints
// Cluster energy, global position, cells and their amplitude fractions are restored
-
+ Int_t j = 0;
for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
{
AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
const Int_t ncells = recPoint->GetMultiplicity();
Int_t ncells_true = 0;
+ if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
+
// cells and their amplitude fractions
UShort_t absIds[ncells];
Double32_t ratios[ncells];
}
if (ncells_true < 1) {
- Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
+ printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",recPoint->GetEnergy(), ncells);
continue;
}
// create a new cluster
//AliAODCaloCluster *clus = new AliAODCaloCluster();
-
- (*clusArray)[i] = new AliAODCaloCluster() ;
- AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(i) ) ;
+ (*clusArray)[j] = new AliAODCaloCluster() ;
+ AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
+ j++;
clus->SetType(AliVCluster::kEMCALClusterv1);
clus->SetE(recPoint->GetEnergy());
clus->SetPosition(g);
Int_t *parentList = recPoint->GetParents(parentMult);
clus->SetLabel(parentList, parentMult);
-// if(parentMult)printf("RecToESD: Labels %d",parentMult);
-// for (Int_t iii = 0; iii < parentMult; iii++) {
-// printf("\t label %d\n",parentList[iii]);
-// }
+ // if(parentMult)printf("RecToESD: Labels %d",parentMult);
+ // for (Int_t iii = 0; iii < parentMult; iii++) {
+ // printf("\t label %d\n",parentList[iii]);
+ // }
+
//Write the second major contributor to each MC cluster.
Int_t iNewLabel ;
for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
}
-
+
newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
//if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
}//positive cell number
}
- //clusArray->Add(clus);
-
} // recPoints loop
+
}