// --- ROOT system ---
-#include "TGeoManager.h"
+#include <TGeoManager.h>
+#include <TH2F.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TPad.h>
// --- ANALYSIS system ---
#include "AliCalorimeterUtils.h"
#include "AliVCluster.h"
#include "AliVCaloCells.h"
#include "AliMixedEvent.h"
+#include "AliAODCaloCluster.h"
// --- Detector ---
#include "AliEMCALGeometry.h"
fRecalculateMatching(kFALSE),
fCutR(20), fCutZ(20),
fCutEta(20), fCutPhi(20),
- fLocMaxCutE(0), fLocMaxCutEDiff(0)
+ fLocMaxCutE(0), fLocMaxCutEDiff(0),
+ fPlotCluster(0)
{
//Ctor
}
+//___________________________________________________________________________
+void AliCalorimeterUtils::SplitEnergy(const Int_t absId1, const Int_t absId2,
+ AliVCluster* cluster,
+ AliVCaloCells* cells,
+ //Float_t & e1, Float_t & e2,
+ AliAODCaloCluster* cluster1,
+ AliAODCaloCluster* cluster2,
+ const Int_t nMax,
+ const Int_t eventNumber)
+{
+
+ // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2
+ // maxima are too close and have common cells, split the energy between the 2
+
+ TH2F* hClusterMap = 0 ;
+ TH2F* hClusterLocMax = 0 ;
+ TH2F* hCluster1 = 0 ;
+ TH2F* hCluster2 = 0 ;
+
+ if(fPlotCluster)
+ {
+ hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
+ hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
+ hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
+ hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
+ hClusterMap ->SetXTitle("column");
+ hClusterMap ->SetYTitle("row");
+ hClusterLocMax ->SetXTitle("column");
+ hClusterLocMax ->SetYTitle("row");
+ hCluster1 ->SetXTitle("column");
+ hCluster1 ->SetYTitle("row");
+ hCluster2 ->SetXTitle("column");
+ hCluster2 ->SetYTitle("row");
+ }
+
+ TString calorimeter = "EMCAL";
+ if(cluster->IsPHOS())
+ {
+ calorimeter="PHOS";
+ printf("AliCalorimeterUtils::SplitEnerg() Not supported for PHOS yet \n");
+ return;
+ }
+
+ const Int_t ncells = cluster->GetNCells();
+ Int_t absIdList[ncells];
+
+ Float_t e1 = 0, e2 = 0 ;
+ Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
+ Float_t eCluster = 0;
+ Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
+ for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
+ {
+ absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
+
+ Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
+ RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
+ eCluster+=ec;
+
+ if(fPlotCluster)
+ {
+ //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
+ sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
+ if(icol > maxCol) maxCol = icol;
+ if(icol < minCol) minCol = icol;
+ if(irow > maxRow) maxRow = irow;
+ if(irow < minRow) minRow = irow;
+ hClusterMap->Fill(icol,irow,ec);
+ }
+
+ }
+
+ // Init counters and variables
+ Int_t ncells1 = 1 ;
+ UShort_t absIdList1[9] ;
+ Double_t fracList1 [9] ;
+ absIdList1[0] = absId1 ;
+ fracList1 [0] = 1. ;
+
+ Float_t ecell1 = cells->GetCellAmplitude(absId1);
+ RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
+ e1 = ecell1;
+
+ Int_t ncells2 = 1 ;
+ UShort_t absIdList2[9] ;
+ Double_t fracList2 [9] ;
+ absIdList2[0] = absId2 ;
+ fracList2 [0] = 1. ;
+
+ Float_t ecell2 = cells->GetCellAmplitude(absId2);
+ RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
+ e2 = ecell2;
+
+ if(fPlotCluster)
+ {
+ Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
+ sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
+ hClusterLocMax->Fill(icol1,irow1,ecell1);
+ sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
+ hClusterLocMax->Fill(icol2,irow2,ecell2);
+ }
+
+ // Very rough way to share the cluster energy
+ Float_t eRemain = (eCluster-ecell1-ecell2)/2;
+ Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
+ Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
+
+ for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
+ {
+ Int_t absId = absIdList[iDigit];
+
+ if(absId==absId1 || absId==absId2 || absId < 0) continue;
+
+ Float_t ecell = cells->GetCellAmplitude(absId);
+ RecalibrateCellAmplitude(ecell, calorimeter, absId);
+
+ if(AreNeighbours(calorimeter, absId1,absId ))
+ {
+ absIdList1[ncells1]= absId;
+
+ if(AreNeighbours(calorimeter, absId2,absId ))
+ {
+ fracList1[ncells1] = shareFraction1;
+ e1 += ecell*shareFraction1;
+ }
+ else
+ {
+ fracList1[ncells1] = 1.;
+ e1 += ecell;
+ }
+
+ ncells1++;
+
+ } // neigbour to cell1
+
+ if(AreNeighbours(calorimeter, absId2,absId ))
+ {
+ absIdList2[ncells2]= absId;
+
+ if(AreNeighbours(calorimeter, absId1,absId ))
+ {
+ fracList2[ncells2] = shareFraction2;
+ e2 += ecell*shareFraction2;
+ }
+ else
+ {
+ fracList2[ncells2] = 1.;
+ e2 += ecell;
+ }
+
+ ncells2++;
+
+ } // neigbour to cell2
+
+ }
+
+ if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2 %f, sum f12 = %f \n",
+ nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
+
+ cluster1->SetE(e1);
+ cluster2->SetE(e2);
+
+ cluster1->SetNCells(ncells1);
+ cluster2->SetNCells(ncells2);
+
+ cluster1->SetCellsAbsId(absIdList1);
+ cluster2->SetCellsAbsId(absIdList2);
+
+ cluster1->SetCellsAmplitudeFraction(fracList1);
+ cluster2->SetCellsAmplitudeFraction(fracList2);
+
+ if(calorimeter=="EMCAL")
+ {
+ GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
+ GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
+ }
+
+ if(fPlotCluster)
+ {
+ //printf("Cells of cluster1: ");
+ for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
+ {
+ //printf(" %d ",absIdList1[iDigit]);
+
+ sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
+
+ if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) )
+ hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
+ else
+ hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
+ }
+
+ //printf(" \n ");
+ //printf("Cells of cluster2: ");
+
+ for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
+ {
+ //printf(" %d ",absIdList2[iDigit]);
+
+ sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
+ if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) )
+ hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
+ else
+ hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
+
+ }
+ //printf(" \n ");
+
+ gStyle->SetPadRightMargin(0.1);
+ gStyle->SetPadLeftMargin(0.1);
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(000000);
+
+ if(maxCol-minCol > maxRow-minRow)
+ {
+ maxRow+= maxCol-minCol;
+ }
+ else
+ {
+ maxCol+= maxRow-minRow;
+ }
+
+ TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
+ c->Divide(2,2);
+ c->cd(1);
+ gPad->SetGridy();
+ gPad->SetGridx();
+ hClusterMap ->SetAxisRange(minCol, maxCol,"X");
+ hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
+ hClusterMap ->Draw("colz");
+ c->cd(2);
+ gPad->SetGridy();
+ gPad->SetGridx();
+ hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
+ hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
+ hClusterLocMax ->Draw("colz");
+ c->cd(3);
+ gPad->SetGridy();
+ gPad->SetGridx();
+ hCluster1 ->SetAxisRange(minCol, maxCol,"X");
+ hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
+ hCluster1 ->Draw("colz");
+ c->cd(4);
+ gPad->SetGridy();
+ gPad->SetGridx();
+ hCluster2 ->SetAxisRange(minCol, maxCol,"X");
+ hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
+ hCluster2 ->Draw("colz");
+
+ if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
+ eventNumber,cluster->E(),nMax,ncells1,ncells2));
+
+ delete c;
+ delete hClusterMap;
+ delete hClusterLocMax;
+ delete hCluster1;
+ delete hCluster2;
+ }
+}
+
//________________________________________________________________________________
void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
TObjArray* clusterArray)
{
//Recalculate track matching
- if (fRecalculateMatching) {
+ if (fRecalculateMatching)
+ {
fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
//AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
//if(esdevent){
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff