#include "AliVEventHandler.h"
#include "AliAODInputHandler.h"
#include "AliOADBContainer.h"
+#include "AliAODMCParticle.h"
// --- EMCAL
#include "AliEMCALAfterBurnerUF.h"
, fFillAODFile(kFALSE), fFillAODHeader(0)
, fFillAODCaloCells(0), fRun(-1)
, fRecoUtils(0), fConfigName("")
+, fOrgClusterCellId()
, fCellLabels(), fCellSecondLabels(), fCellTime()
, fCellMatchdEta(), fCellMatchdPhi()
, fRecalibrateWithClusterTime(0)
, fMaxEvent(0), fDoTrackMatching(kFALSE)
, fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
+, fRejectBelowThreshold(kFALSE)
, fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
-, fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
+, fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
, fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
, fCentralityClass(""), fSelectEMCALEvent(0)
, fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
+, fSetCellMCLabelFromCluster(0)
+, fRemapMCLabelForAODs(0)
+, fInputFromFilter(0)
{
// Constructor
for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
- for(Int_t j = 0; j < 24*48*11; j++)
- {
- fCellLabels[j] = -1;
- fCellSecondLabels[j] = -1;
- fCellTime[j] = 0.;
- fCellMatchdEta[j] = -999;
- fCellMatchdPhi[j] = -999;
- }
+
+ ResetArrays();
fCentralityBin[0] = fCentralityBin[1]=-1;
, fFillAODFile(kFALSE), fFillAODHeader(0)
, fFillAODCaloCells(0), fRun(-1)
, fRecoUtils(0), fConfigName("")
+, fOrgClusterCellId()
, fCellLabels(), fCellSecondLabels(), fCellTime()
, fCellMatchdEta(), fCellMatchdPhi()
, fRecalibrateWithClusterTime(0)
, fMaxEvent(0), fDoTrackMatching(kFALSE)
, fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
+, fRejectBelowThreshold(kFALSE)
, fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
-, fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
+, fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
, fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
, fCentralityClass(""), fSelectEMCALEvent(0)
, fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
+, fSetCellMCLabelFromCluster(0)
+, fRemapMCLabelForAODs(0)
+, fInputFromFilter(0)
{
// Constructor
for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
- for(Int_t j = 0; j < 24*48*11; j++)
- {
- fCellLabels[j] = -1;
- fCellSecondLabels[j] = -1;
- fCellTime[j] = 0.;
- fCellMatchdEta[j] = -999;
- fCellMatchdPhi[j] = -999;
- }
-
+
+ ResetArrays();
+
fCentralityBin[0] = fCentralityBin[1]=-1;
}
if(fEMCALEnergyCut <= 0) return kTRUE; // accept
- Int_t nCluster = InputEvent() -> GetNumberOfCaloClusters();
- AliVCaloCells * caloCell = InputEvent() -> GetEMCALCells();
- Int_t bc = InputEvent() -> GetBunchCrossNumber();
-
+ Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
+ AliVCaloCells * caloCell = fEvent -> GetEMCALCells();
+ Int_t bc = fEvent -> GetBunchCrossNumber();
+
for(Int_t icalo = 0; icalo < nCluster; icalo++)
{
- AliVCluster *clus = (AliVCluster*) (InputEvent()->GetCaloCluster(icalo));
+ AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
if (fDebug > 0)
printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Accept : E %2.2f > %2.2f, nCells %d > %d \n",
- clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
+ clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut);
return kTRUE;
}
return kFALSE;
-}
+}
//_______________________________________________
void AliAnalysisTaskEMCALClusterize::AccessOADB()
TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
+ //If it did not exist for this run, get closes one
+ if (!htd)
+ {
+ AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
+ // let's get the closest runnumber instead then..
+ Int_t lower = 0;
+ Int_t ic = 0;
+ Int_t maxEntry = contRFTD->GetNumberOfEntries();
+
+ while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
+ lower = ic;
+ ic++;
+ }
+
+ Int_t closest = lower;
+ if ( (ic<maxEntry) &&
+ (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
+ closest = ic;
+ }
+
+ AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
+ htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
+ }
+
if(htd)
{
printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
if(trecal)
{
- TString passTmp = pass;
- if(pass!="pass1" && pass!="pass2") passTmp = "pass2"; // TEMPORARY FIX FOR LHC11a analysis
- TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passTmp);
+ TString passM = pass;
+ if(pass=="spc_calo") passM = "pass1";
+ TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
if(trecalpass)
{
// Also check if the quality of the event is good if not reject it
fEvent = 0x0;
-
- //Process events if there is a high energy cluster
- if(!AcceptEventEMCAL()) return ;
-
+
AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
Int_t eventN = Entry();
if(aodIH) eventN = aodIH->GetReadEntry();
AODEvent()->GetEMCALCells()->GetNumberOfCells());
}
}
+ else if(fInputFromFilter)
+ {
+ //printf("Get Input From Filtered AOD\n");
+ fEvent = AODEvent();
+ }
else
{
fEvent = InputEvent();
return ;
}
+ //Process events if there is a high energy cluster
+ if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
+
//-------------------------------------------------------------------------------------
// Reject events if LED was firing, use only for LHC11a data
// Reject event if triggered by exotic cell and remove exotic cells if not triggered
AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
if(aodIH && aodIH->GetEventToMerge()) //Embedding
nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
-
+
+ ResetArrays();
+
for (Int_t i = 0; i < nClusters; i++)
{
AliVCluster *clus = 0;
if(!clus) return;
if(clus->IsEMCAL())
- {
+ {
Int_t label = clus->GetLabel();
Int_t label2 = -1 ;
- //printf("Org cluster E %f, Time %e, Id = ", clus->E(), clus->GetTOF() );
+ //printf("Org cluster E %f, Time %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label );
+ //printf("Original list of labels from old cluster : \n");
+ //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
+
if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
UShort_t * index = clus->GetCellsAbsId() ;
for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
{
+ //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell]));
+ fOrgClusterCellId[index[icell]] = i;
fCellLabels[index[icell]] = label;
fCellSecondLabels[index[icell]] = label2;
- fCellTime[index[icell]] = clus->GetTOF();
+ fCellTime[index[icell]] = clus->GetTOF();
fCellMatchdEta[index[icell]] = clus->GetTrackDz();
fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
- //printf(" %d,", index[icell] );
}
nClustersOrg++;
}
Float_t amp = -1;
Double_t time = -1;
- AliVCaloCells *cells = fEvent->GetEMCALCells();
+ AliVCaloCells *cells = fEvent->GetEMCALCells();
Int_t bc = InputEvent()->GetBunchCrossNumber();
+
for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
{
// Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
if( !accept )
{
- fCellLabels[id] =-1; //reset the entry in the array for next event
- fCellSecondLabels[id]=-1; //reset the entry in the array for next event
- fCellTime[id] = 0.;
- fCellMatchdEta[id] =-999;
- fCellMatchdPhi[id] =-999;
if( DebugLevel() > 2 )
printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
continue;
}
- //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);
+ Int_t mcLabel = cells->GetMCLabel(icell);
+ //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel);
+
+ //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]);
+ if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
+ else if( fSetCellMCLabelFromCluster == 0 && fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel);
+ else mcLabel = -1; // found later
- fCellLabels[id] =-1; //reset the entry in the array for next event
+ //printf("\t new label %d\n",mcLabel);
+ // Create the digit, put a fake primary deposited energy to trick the clusterizer
+ // when checking the most likely primary
+
+ Float_t efrac = cells->GetEFraction(icell);
+
+ //When checking the MC of digits, give weight to cells with embedded signal
+ if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
+
+ //printf("******* Cell %d, id %d, e %f, fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel);
+
+ new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
+ // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that
+ // we give more weight to the MC label of the cell with highest energy in the cluster
+
idigit++;
}
fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
}
}
-
- //Reset the array with second labels for this event
- memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
-
}
//_____________________________________________________
Double_t cellAmplitude = 0;
Double_t cellTime = 0;
Short_t cellNumber = 0;
- Short_t cellMCLabel = 0;
+ Int_t cellMCLabel = 0;
Double_t cellEFrac = 0;
Int_t nClustersOrg = 0;
//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
- fCellTime[cellNumber] = 0.;
- fCellMatchdEta[cellNumber] =-999;
- fCellMatchdPhi[cellNumber] =-999;
+ fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
continue;
- }
cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
//Trigger
header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
- if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
- else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
+
+ header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
header->SetTriggerMask(fEvent->GetTriggerMask());
header->SetTriggerCluster(fEvent->GetTriggerCluster());
header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
- Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
+ Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
Float_t diamcov[3];
fEvent->GetDiamondCovXY(diamcov);
header->SetDiamond(diamxy,diamcov);
new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
- if(DebugLevel() > 1 )
- printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
+ if(DebugLevel() > 1 )
+ printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f, mc label %d \n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel());
} // cluster loop
else if (pass.Contains("ass3")) return TString("pass3");
else if (pass.Contains("ass4")) return TString("pass4");
else if (pass.Contains("ass5")) return TString("pass5");
+ else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
+ else if (pass.Contains("calo") || pass.Contains("high_lumi"))
+ {
+ printf("AliAnalysisTaskEMCALClusterize::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n");
+ return TString("pass1");
+ }
// No condition fullfilled, give a default value
printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
if(fMaxEvent <= 0) fMaxEvent = 1000000000;
if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
-
+ fRejectBelowThreshold = kFALSE;
+
//Centrality
if(fCentralityClass == "") fCentralityClass = "V0M";
fCentralityBin[0] = clus->fCentralityBin[0];
fCentralityBin[1] = clus->fCentralityBin[1];
}
-
- // Init geometry, I do not like much to do it like this ...
- if(fImportGeometryFromFile && !gGeoManager)
- {
- if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
- printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
- TGeoManager::Import(fImportGeometryFilePath) ;
- }
-
+
}
//_______________________________________________________
fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
}//end of loop over parameters
-
+ fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
fClusterizer->InitClusterUnfolding();
}// to unfold
if (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
else fGeomName = "EMCAL_COMPLETE12SMV1";
- printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
+ printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",
+ fGeomName.Data(),runnumber);
}
fGeom = AliEMCALGeometry::GetInstance(fGeomName);
+ // Init geometry, I do not like much to do it like this ...
+ if(fImportGeometryFromFile && !gGeoManager)
+ {
+ if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
+ {
+ // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
+ if (runnumber < 140000 &&
+ runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root";
+ if (runnumber >= 140000 &&
+ runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root";
+ else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013
+ }
+ printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Import %s\n",fImportGeometryFilePath.Data());
+ TGeoManager::Import(fImportGeometryFilePath) ;
+ }
+
if(fDebug > 0)
{
printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
}
- // TString triggerclasses = "";
- // if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses();
- // else triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
- // //
+ // TString triggerclasses = event->GetFiredTriggerClasses();
// printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
// if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
// return;
if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
}
- TString triggerclasses = "";
-
- AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
- if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses();
- else triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
+ TString triggerclasses = fEvent->GetFiredTriggerClasses();
Int_t ncellcut = 21;
if(triggerclasses.Contains("EMC")) ncellcut = 35;
absIds[ncellsTrue] = digit->GetId();
ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
-
+
// In case of unfolding, remove digits with unfolded energy too low
if(fSelectCell)
{
if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
{
- if(DebugLevel() > 1)
+ if(DebugLevel() > 1)
{
printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
recPoint->GetEnergiesList()[c],digit->GetAmplitude());
}// Select cells
//Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
- ncellsTrue++;
clusterE +=recPoint->GetEnergiesList()[c];
// In case of embedding, fill ratio with amount of signal,
}//Embedding
+ ncellsTrue++;
+
}// cluster cell loop
if (ncellsTrue < 1)
}
- //MC
- Int_t parentMult = 0;
- Int_t *parentList = recPoint->GetParents(parentMult);
- clus->SetLabel(parentList, parentMult);
+ // MC
+
+ if ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ;
+ else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ;
+ else
+ {
+ // Normal case, trust what the clusterizer has found
+ Int_t parentMult = 0;
+ Int_t *parentList = recPoint->GetParents(parentMult);
+ clus->SetLabel(parentList, parentMult);
+// printf("Label list : ");
+// for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]);
+// printf("\n");
+ }
+
+ } // recPoints loop
+
+}
+
+//___________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label)
+{
+ // MC label for Cells not remapped after ESD filtering, do it here.
+
+ if(label < 0) return ;
+
+ AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
+ if(!evt) return ;
+
+ TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
+ if(!arr) return ;
+
+ if(label < arr->GetEntriesFast())
+ {
+ AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
+ if(!particle) return ;
+
+ if(label == particle->Label()) return ; // label already OK
+ //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
+ }
+ //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
+
+ // loop on the particles list and check if there is one with the same label
+ for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
+ {
+ AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
+ if(!particle) continue ;
+
+ if(label == particle->Label())
+ {
+ label = ind;
+ //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
+ return;
+ }
+ }
+
+ label = -1;
+
+ //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
+
+}
+
+//________________________________________________
+void AliAnalysisTaskEMCALClusterize::ResetArrays()
+{
+ // Reset arrays containing information for all possible cells
+ for(Int_t j = 0; j < 12672; j++)
+ {
+ fOrgClusterCellId[j] = -1;
+ fCellLabels[j] = -1;
+ fCellSecondLabels[j] = -1;
+ fCellTime[j] = 0.;
+ fCellMatchdEta[j] = -999;
+ fCellMatchdPhi[j] = -999;
+ }
+}
+
+//_____________________________________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint * recPoint,
+ AliAODCaloCluster * clus)
+{
+ // Set the cluster MC label, the digizer was filled with most likely MC label for all cells in original cluster
+ // Now check the second most likely MC label and add it to the new cluster
+
+ Int_t parentMult = 0;
+ Int_t *parentList = recPoint->GetParents(parentMult);
+ clus->SetLabel(parentList, parentMult);
+
+ //Write the second major contributor to each MC cluster.
+ Int_t iNewLabel ;
+ for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
+ {
- //Write the second major contributor to each MC cluster.
- Int_t iNewLabel ;
- for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
+ Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
+ if(idCell>=0)
{
-
- Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
- if(idCell>=0)
+ iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
+ for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
{
- iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
- for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
- {
- if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
+ if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
- }
- if (iNewLabel == 1)
+ }
+ if (iNewLabel == 1)
+ {
+ Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
+ for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
{
- Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
- for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
- {
- newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
- }
-
- newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
- clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
- delete [] newLabelArray;
+ newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
}
- }//positive cell number
+
+ newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
+ clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
+ delete [] newLabelArray;
+ }
+ }//positive cell number
+ }
+}
+
+//___________________________________________________________________________________________________
+void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus)
+{
+ // Get the original clusters that contribute to the new cluster, assign the labels of such clusters
+ // to the new cluster.
+ // Only approximatedly valid when output are V1 clusters, handle with care
+
+ TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
+ clArray.Reset();
+ Int_t nClu = 0;
+ Int_t nLabTotOrg = 0;
+ Float_t emax = -1;
+ Int_t idMax = -1;
+
+ AliVEvent * event = fEvent;
+
+ //In case of embedding the MC signal is in the added event
+ AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
+ if(aodIH && aodIH->GetEventToMerge()) //Embedding
+ event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
+
+
+ //Find the clusters that originally had the cells
+ for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
+ {
+ Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
+
+ if(idCell>=0)
+ {
+ Int_t idCluster = fOrgClusterCellId[idCell];
+
+ Bool_t set = kTRUE;
+ for(Int_t icl =0; icl < nClu; icl++)
+ {
+ if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
+ if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
+ // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
+ // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
+ }
+ if( set && idCluster >= 0)
+ {
+ clArray.SetAt(idCluster,nClu++);
+ //printf("******** idCluster %d \n",idCluster);
+ nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
+
+ //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
+
+ //Search highest E cluster
+ AliVCluster * clOrg = event->GetCaloCluster(idCluster);
+ //printf("\t E %f\n",clOrg->E());
+ if(emax < clOrg->E())
+ {
+ emax = clOrg->E();
+ idMax = idCluster;
+ }
+ }
+ }
+ }// cell loop
+
+ if(nClu==0 || nLabTotOrg == 0)
+ {
+ //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f, n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
+ //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
+ //printf("\n");
+ }
+
+ // Put the first in the list the cluster with highest energy
+ if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
+ {
+ Int_t maxIndex = -1;
+ Int_t firstCluster = ((Int_t)clArray.GetAt(0));
+ for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
+ {
+ if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
}
- } // recPoints loop
+ if(firstCluster >=0 && idMax >=0)
+ {
+ clArray.SetAt(idMax,0);
+ clArray.SetAt(firstCluster,maxIndex);
+ }
+ }
+
+ // Get the labels list in the original clusters, assign all to the new cluster
+ TArrayI clMCArray(nLabTotOrg) ;
+ clMCArray.Reset();
+ Int_t nLabTot = 0;
+ for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
+ {
+ Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
+ //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
+ AliVCluster * clOrg = event->GetCaloCluster(idCluster);
+ Int_t nLab = clOrg->GetNLabels();
+
+ for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
+ {
+ Int_t lab = clOrg->GetLabelAt(iLab) ;
+ if(lab>=0)
+ {
+ Bool_t set = kTRUE;
+ //printf("\t \t Set Label %d \n", lab);
+ for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
+ {
+ if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
+ //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
+ // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
+ }
+ if( set ) clMCArray.SetAt(lab,nLabTot++);
+ }
+ }
+ }// cluster loop
+
+ // Set the final list of labels
+
+ clus->SetLabel(clMCArray.GetArray(), nLabTot);
+
+// printf("Final list of labels for new cluster : \n");
+// for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
+// {
+// printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
+// Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
+// printf(" org %d ",label);
+// RemapMCLabelForAODs(label);
+// printf(" new %d \n",label);
+// }
+// for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff