#include "AliEMCALGeometry.h"
#include "AliVCluster.h"
#include "AliVCaloCells.h"
+#include "AliEMCALRecoUtils.h"
+
//#include "AliEMCALAodCluster.h"
//#include "AliEMCALCalibData.h"
//__________________________________________________
AliAnalysisTaskEMCALPi0CalibSelection::AliAnalysisTaskEMCALPi0CalibSelection(const char* name) :
AliAnalysisTaskSE(name),fEMCALGeo(0x0),//fCalibData(0x0),
- fEmin(0.5), fEmax(15.), fMinNCells(2), fGroupNCells(0),
- fLogWeight(4.5), fCopyAOD(kFALSE), fEMCALGeoName("EMCAL_FIRSTYEAR"),
+ fEmin(0.5), fEmax(15.), fAsyCut(1.),fMinNCells(2), fGroupNCells(0),
+ fLogWeight(4.5), fCopyAOD(kFALSE), fSameSM(kFALSE), fOldAOD(kFALSE),
+ fEMCALGeoName("EMCAL_FIRSTYEAR"), fNCellsFromEMCALBorder(0),
fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
- fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
- fNbins(300), fMinBin(0.), fMaxBin(300.),
- fOutputContainer(0x0),fHmgg(0x0), fhNEvents(0x0),fCuts(0x0)
+ fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
+ fRecoUtils(new AliEMCALRecoUtils),
+ fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
+ fHmgg(0x0), fHmggDifferentSM(0x0),
+ fHOpeningAngle(0x0), fHOpeningAngleDifferentSM(0x0),
+ fHIncidentAngle(0x0), fHIncidentAngleDifferentSM(0x0),
+ fHAsymmetry(0x0), fHAsymmetryDifferentSM(0x0),
+ fhNEvents(0x0),fCuts(0x0)
{
//Named constructor which should be used.
for(Int_t iMod=0; iMod < 12; iMod++) {
for(Int_t iX=0; iX<24; iX++) {
for(Int_t iZ=0; iZ<48; iZ++) {
- fHmpi0[iMod][iX][iZ]=0;
+ fHmpi0[iMod][iZ][iX]=0;
}
}
}
+ for(Int_t iSM=0; iSM<4; iSM++) {
+ fHmggSM[iSM] =0;
+ fHmggPairSM[iSM] =0;
+ fHOpeningAngleSM[iSM] =0;
+ fHOpeningAnglePairSM[iSM] =0;
+ fHAsymmetrySM[iSM] =0;
+ fHAsymmetryPairSM[iSM] =0;
+ fHIncidentAngleSM[iSM] =0;
+ fHIncidentAnglePairSM[iSM]=0;
+ fhTowerDecayPhotonHit[iSM] =0;
+ fhTowerDecayPhotonEnergy[iSM]=0;
+ fhTowerDecayPhotonAsymmetry[iSM]=0;
+ }
+
DefineOutput(1, TList::Class());
DefineOutput(2, TList::Class()); // will contain cuts or local params
fEMCALRecalibrationFactors->Clear();
delete fEMCALRecalibrationFactors;
}
+
+ if(fRecoUtils) delete fRecoUtils ;
+
}
//_____________________________________________________
// Local Initialization
// Create cuts/param objects and publish to slot
-
- char onePar[255] ;
+ const Int_t buffersize = 255;
+ char onePar[buffersize] ;
fCuts = new TList();
- sprintf(onePar,"Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d;", fEmin,fEmax, fMinNCells) ;
+ snprintf(onePar,buffersize, "Custer cuts: %2.2f < E < %2.2f GeV; min number of cells %d; Assymetry cut %1.2f", fEmin,fEmax, fMinNCells, fAsyCut) ;
+ fCuts->Add(new TObjString(onePar));
+ snprintf(onePar,buffersize, "Group %d cells;", fGroupNCells) ;
fCuts->Add(new TObjString(onePar));
- sprintf(onePar,"Group %d cells;", fGroupNCells) ;
+ snprintf(onePar,buffersize, "Cluster maximal cell away from border at least %d cells;", fNCellsFromEMCALBorder) ;
fCuts->Add(new TObjString(onePar));
- sprintf(onePar,"Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
+ snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
fCuts->Add(new TObjString(onePar));
- sprintf(onePar,"Switchs: Remove Bad Channels? %d; Copy AODs? %d; Recalibrate %d?",fRemoveBadChannels,fCopyAOD,fRecalibration) ;
+ snprintf(onePar,buffersize, "Switchs: Remove Bad Channels? %d; Copy AODs? %d; Recalibrate %d?, Analyze Old AODs? %d, Mass per channel same SM clusters? %d ",
+ fRemoveBadChannels,fCopyAOD,fRecalibration, fOldAOD, fSameSM) ;
fCuts->Add(new TObjString(onePar));
- sprintf(onePar,"EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
+ snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
fCuts->Add(new TObjString(onePar));
// Post Data
// Copy AOD header, vertex, CaloClusters and CaloCells to output AOD
AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ if(!aod) {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
+ return;
+ }
+
// set arrays and pointers
Float_t posF[3];
Double_t pos[3];
// Access to the AOD container of clusters
TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
Int_t jClusters=0;
- printf("nCaloClus %d\n",nCaloClus);
+ //printf("nCaloClus %d\n",nCaloClus);
for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
AliAODCaloCluster * cluster = aod->GetCaloCluster(iClust);
//Check if it is a EMCAL cluster
- if(!cluster->IsEMCAL()) continue ;
- printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
- if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
- printf("copy\n");
+ if(!IsEMCALCluster(cluster)) continue ;
+ //printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
+ //if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
+ //printf("copy\n");
Int_t id = cluster->GetID();
Float_t energy = cluster->E();
cluster->GetPosition(posF);
// Copy Header, Vertex, CaloClusters and CaloCells from ESDs to AODs
AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
+ if(!esd) {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::CreateAODFromAOD() - This event does not contain AODs?");
+ return;
+ }
+
// set arrays and pointers
Float_t posF[3];
Double_t pos[3];
// Access to the AOD container of clusters
TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
Int_t jClusters=0;
- printf("nCaloClus %d\n",nCaloClus);
+ //printf("nCaloClus %d\n",nCaloClus);
for (Int_t iClust=0; iClust<nCaloClus; ++iClust) {
AliESDCaloCluster * cluster = esd->GetCaloCluster(iClust);
//Check which calorimeter information we want to keep.
- if(!cluster->IsEMCAL()) continue ;
- printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
+ if(!IsEMCALCluster(cluster)) continue ;
+ //printf("EMCAL cluster %d, ncells %d\n",iClust, cluster->GetNCells());
if(ClusterContainsBadChannel(cluster->GetCellsAbsId(), cluster->GetNCells())) continue;
- printf("copy\n");
+ //printf("copy\n");
Int_t id = cluster->GetID();
Float_t energy = cluster->E();
}
+//_________________________________________________________________
+Int_t AliAnalysisTaskEMCALPi0CalibSelection::GetEMCALClusters(AliVEvent * event, TRefArray *clusters) const
+{
+ // fills the provided TRefArray with all found emcal clusters
+
+ clusters->Clear();
+ AliVCluster *cl = 0;
+ Bool_t first = kTRUE;
+ for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) {
+ if ( (cl = event->GetCaloCluster(i)) ) {
+ if (IsEMCALCluster(cl)){
+ if(first) {
+ new (clusters) TRefArray(TProcessID::GetProcessWithUID(cl));
+ first=kFALSE;
+ }
+ clusters->Add(cl);
+ //printf("IsEMCal cluster %d, E %2.3f Size: %d \n",i,cl->E(),clusters->GetEntriesFast());
+ }
+ }
+ }
+ return clusters->GetEntriesFast();
+}
+
+
+//____________________________________________________________________________
+Bool_t AliAnalysisTaskEMCALPi0CalibSelection::IsEMCALCluster(AliVCluster* cluster) const {
+ // Check if it is a cluster from EMCAL. For old AODs cluster type has
+ // different number and need to patch here
+
+ if(fOldAOD)
+ {
+ if (cluster->GetType() == 2) return kTRUE;
+ else return kFALSE;
+ }
+ else
+ {
+ return cluster->IsEMCAL();
+ }
+
+}
+
+
//__________________________________________________
void AliAnalysisTaskEMCALPi0CalibSelection::UserCreateOutputObjects()
{
//Create output container, init geometry and calibration
fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
-
- fOutputContainer = new TList();
- char hname[128], htitl[128];
+ fOutputContainer = new TList();
+ const Int_t buffersize = 255;
+ char hname[buffersize], htitl[buffersize];
for(Int_t iMod=0; iMod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); iMod++) {
for(Int_t iRow=0; iRow<AliEMCALGeoParams::fgkEMCALRows; iRow++) {
for(Int_t iCol=0; iCol<AliEMCALGeoParams::fgkEMCALCols; iCol++) {
- sprintf(hname,"%d_%d_%d",iMod,iCol,iRow);
- sprintf(htitl,"Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
+ snprintf(hname,buffersize, "%d_%d_%d",iMod,iCol,iRow);
+ snprintf(htitl,buffersize, "Two-gamma inv. mass for super mod %d, cell(col,row)=(%d,%d)",iMod,iCol,iRow);
fHmpi0[iMod][iCol][iRow] = new TH1F(hname,htitl,fNbins,fMinBin,fMaxBin);
fOutputContainer->Add(fHmpi0[iMod][iCol][iRow]);
}
fHmgg = new TH2F("hmgg","2-cluster invariant mass",fNbins,fMinBin,fMaxBin,100,0,10);
fHmgg->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
fHmgg->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
-
fOutputContainer->Add(fHmgg);
+
+ fHmggDifferentSM = new TH2F("hmggDifferentSM","2-cluster invariant mass, different SM",fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggDifferentSM->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggDifferentSM);
+
+ fHOpeningAngle = new TH2F("hopang","2-cluster opening angle",100,0.,50.,100,0,10);
+ fHOpeningAngle->SetXTitle("#alpha_{#gamma #gamma}");
+ fHOpeningAngle->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAngle);
+
+ fHOpeningAngleDifferentSM = new TH2F("hopangDifferentSM","2-cluster opening angle, different SM",100,0,50.,100,0,10);
+ fHOpeningAngleDifferentSM->SetXTitle("#alpha_{#gamma #gamma}");
+ fHOpeningAngleDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAngleDifferentSM);
+
+ fHIncidentAngle = new TH2F("hinang","#gamma incident angle in SM",100,0.,20.,100,0,10);
+ fHIncidentAngle->SetXTitle("#alpha_{#gamma SM center}");
+ fHIncidentAngle->SetYTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fHIncidentAngle);
+
+ fHIncidentAngleDifferentSM = new TH2F("hinangDifferentSM","#gamma incident angle in SM, different SM pair",100,0,20.,100,0,10);
+ fHIncidentAngleDifferentSM->SetXTitle("#alpha_{#gamma - SM center}");
+ fHIncidentAngleDifferentSM->SetYTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fHIncidentAngleDifferentSM);
+
+ fHAsymmetry = new TH2F("hasym","2-cluster opening angle",100,0.,1.,100,0,10);
+ fHAsymmetry->SetXTitle("a");
+ fHAsymmetry->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetry);
+
+ fHAsymmetryDifferentSM = new TH2F("hasymDifferentSM","2-cluster opening angle, different SM",100,0,1.,100,0,10);
+ fHAsymmetryDifferentSM->SetXTitle("a");
+ fHAsymmetryDifferentSM->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetryDifferentSM);
+
+
+ TString pairname[] = {"A side (0-2)", "C side (1-3)","Row 0 (0-1)", "Row 1 (2-3)"};
+
+ for(Int_t iSM=0; iSM<4; iSM++) {
+
+ snprintf(hname, buffersize, "hmgg_SM%d",iSM);
+ snprintf(htitl, buffersize, "Two-gamma inv. mass for super mod %d",iSM);
+ fHmggSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggSM[iSM]);
+
+ snprintf(hname,buffersize, "hmgg_PairSM%d",iSM);
+ snprintf(htitl,buffersize, "Two-gamma inv. mass for SM pair: %s",pairname[iSM].Data());
+ fHmggPairSM[iSM] = new TH2F(hname,htitl,fNbins,fMinBin,fMaxBin,100,0,10);
+ fHmggPairSM[iSM]->SetXTitle("m_{#gamma #gamma} (MeV/c^{2})");
+ fHmggPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHmggPairSM[iSM]);
+
+
+ snprintf(hname, buffersize, "hopang_SM%d",iSM);
+ snprintf(htitl, buffersize, "Opening angle for super mod %d",iSM);
+ fHOpeningAngleSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
+ fHOpeningAngleSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
+ fHOpeningAngleSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAngleSM[iSM]);
+
+ snprintf(hname,buffersize, "hopang_PairSM%d",iSM);
+ snprintf(htitl,buffersize, "Opening angle for SM pair: %s",pairname[iSM].Data());
+ fHOpeningAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,50.,100,0,10);
+ fHOpeningAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma #gamma} (deg)");
+ fHOpeningAnglePairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHOpeningAnglePairSM[iSM]);
+
+ snprintf(hname, buffersize, "hinang_SM%d",iSM);
+ snprintf(htitl, buffersize, "Incident angle for super mod %d",iSM);
+ fHIncidentAngleSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
+ fHIncidentAngleSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
+ fHIncidentAngleSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fHIncidentAngleSM[iSM]);
+
+ snprintf(hname,buffersize, "hinang_PairSM%d",iSM);
+ snprintf(htitl,buffersize, "Incident angle for SM pair: %s",pairname[iSM].Data());
+ fHIncidentAnglePairSM[iSM] = new TH2F(hname,htitl,100,0.,20.,100,0,10);
+ fHIncidentAnglePairSM[iSM]->SetXTitle("#alpha_{#gamma - SM center} (deg)");
+ fHIncidentAnglePairSM[iSM]->SetYTitle("p_{T #gamma} (GeV/c)");
+ fOutputContainer->Add(fHIncidentAnglePairSM[iSM]);
+
+ snprintf(hname, buffersize, "hasym_SM%d",iSM);
+ snprintf(htitl, buffersize, "asymmetry for super mod %d",iSM);
+ fHAsymmetrySM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
+ fHAsymmetrySM[iSM]->SetXTitle("a");
+ fHAsymmetrySM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetrySM[iSM]);
+
+ snprintf(hname,buffersize, "hasym_PairSM%d",iSM);
+ snprintf(htitl,buffersize, "Asymmetry for SM pair: %s",pairname[iSM].Data());
+ fHAsymmetryPairSM[iSM] = new TH2F(hname,htitl,100,0.,1.,100,0,10);
+ fHAsymmetryPairSM[iSM]->SetXTitle("a");
+ fHAsymmetryPairSM[iSM]->SetYTitle("p_{T #gamma #gamma} (GeV/c)");
+ fOutputContainer->Add(fHAsymmetryPairSM[iSM]);
+
+
+ Int_t colmax = 48;
+ Int_t rowmax = 24;
+
+ fhTowerDecayPhotonHit[iSM] = new TH2F (Form("hTowerDecPhotonHit_Mod%d",iSM),Form("Entries in grid of cells in Module %d",iSM),
+ colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
+ fhTowerDecayPhotonHit[iSM]->SetYTitle("row (phi direction)");
+ fhTowerDecayPhotonHit[iSM]->SetXTitle("column (eta direction)");
+ fOutputContainer->Add(fhTowerDecayPhotonHit[iSM]);
+
+ fhTowerDecayPhotonEnergy[iSM] = new TH2F (Form("hTowerDecPhotonEnergy_Mod%d",iSM),Form("Accumulated energy in grid of cells in Module %d",iSM),
+ colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
+ fhTowerDecayPhotonEnergy[iSM]->SetYTitle("row (phi direction)");
+ fhTowerDecayPhotonEnergy[iSM]->SetXTitle("column (eta direction)");
+ fOutputContainer->Add(fhTowerDecayPhotonEnergy[iSM]);
+
+ fhTowerDecayPhotonAsymmetry[iSM] = new TH2F (Form("hTowerDecPhotonAsymmetry_Mod%d",iSM),Form("Accumulated asymmetry in grid of cells in Module %d",iSM),
+ colmax+2,-1.5,colmax+0.5, rowmax+2,-1.5,rowmax+0.5);
+ fhTowerDecayPhotonAsymmetry[iSM]->SetYTitle("row (phi direction)");
+ fhTowerDecayPhotonAsymmetry[iSM]->SetXTitle("column (eta direction)");
+ fOutputContainer->Add(fhTowerDecayPhotonAsymmetry[iSM]);
+
+ }
+
+
fhNEvents = new TH1I("hNEvents", "Number of analyzed events" , 1 , 0 , 1 ) ;
fOutputContainer->Add(fhNEvents);
if(kAOD){
//Input are ESDs
aod = dynamic_cast<AliAODEvent*>(InputEvent());
+ if(!aod) {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
+ return;
+ }
+
// Create new AOD with only CaloClusters and CaloCells
if(fCopyAOD) CreateAODFromAOD();
}
else if(kESD) {
//Input are ESDs
aod = AODEvent();
+ if(!aod) {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain AODs?");
+ return;
+ }
+
// Create AOD with CaloClusters and use it as input.
// If filtering task is already executed, this is not needed.
if(fCopyAOD) CreateAODFromESD();
Int_t runNum = aod->GetRunNumber();
if(DebugLevel() > 1) printf("Run number: %d\n",runNum);
- //FIXME Not neede the matrices for the moment MEFIX
+ //FIXME Not need the matrices for the moment MEFIX
//Get the matrix with geometry information
//Still not implemented in AOD, just a workaround to be able to work at least with ESDs
// if(!strcmp(InputEvent()->GetName(),"AliAODEvent")) {
// else{
// if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Load Misaligned matrices. \n");
// AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
+ // if(!esd) {
+ // printf("AliAnalysisTaskEMCALPi0CalibSelection::UserExec() - This event does not contain ESDs?");
+ // return;
+ // }
// for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
// if(esd->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
// }
TLorentzVector p1;
TLorentzVector p2;
+// TLorentzVector p11;
+// TLorentzVector p22;
+
TLorentzVector p12;
TRefArray * caloClustersArr = new TRefArray();
- aod->GetEMCALClusters(caloClustersArr);
-
+ if(!fOldAOD) aod->GetEMCALClusters(caloClustersArr);
+ else GetEMCALClusters(aod,caloClustersArr);
+
const Int_t kNumberOfEMCALClusters = caloClustersArr->GetEntries() ;
if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection - N CaloClusters: %d \n", kNumberOfEMCALClusters);
// EMCAL cells
AliAODCaloCells *emCells = aod->GetEMCALCells();
-
+
// loop over EMCAL clusters
for(Int_t iClu=0; iClu<kNumberOfEMCALClusters; iClu++) {
}
//clu1.GetMomentum(p1,v);
- c1->GetMomentum(p1,v);
-
- //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
- //MaxEnergyCellPos(emCells,&clu1,iSupMod1,ieta1,iphi1);
- MaxEnergyCellPos(emCells,c1,iSupMod1,ieta1,iphi1);
+ // Correct Non-Linearity
+ c1->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c1));
+ //printf("\t e1 org %2.3f, e1 cor %2.3f \n",e1ii,c1->E());
+
+ //c1->GetMomentum(p1,v);
+ //printf("\t cor: e %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f\n",p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
+
+ //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate
+ //Float_t pos[3];
+ //c1->GetPosition(pos);
+ //printf("Before Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
+ GetMaxEnergyCellPosAndClusterPos(emCells,c1,iSupMod1,ieta1,iphi1);
+ //printf("i1 %d, corr1 %2.3f, e1 %2.3f, , ecorr1 %2.3f, pt %2.3f, px %2.3f, py %2.3f, pz %2.3f,\n",iClu, 1./corrFac, e1, p1.E(), p1.Pt(),p1.Px(),p1.Py(),p1.Pz());
+ //c1->GetPosition(pos);
+ //printf("After Alignment: e %2.4f, x %2.4f, y %2.4f , z %2.4f\n",c1->E(),pos[0], pos[1],pos[2]);
+
+ c1->GetMomentum(p1,v);
+
for (Int_t jClu=iClu; jClu<kNumberOfEMCALClusters; jClu++) {
AliAODCaloCluster *c2 = (AliAODCaloCluster *) caloClustersArr->At(jClu);
if(c2->IsEqual(c1)) continue;
Float_t e2ii = c2->E();
- //clu2.GetMomentum(p2,v);
+ //Correct Non-Linearity
+ c2->SetE(fRecoUtils->CorrectClusterEnergyLinearity(c2));
+
+ //Get tower with maximum energy and fill in the end the pi0 histogram for this cell, recalculate cluster position and recalibrate
+ GetMaxEnergyCellPosAndClusterPos(emCells,c2,iSupMod2,ieta2,iphi2);
+
c2->GetMomentum(p2,v);
-
- //Get tower with maximum energy and fill in the end the pi0 histogram for this cell.
- //MaxEnergyCellPos(emCells,&clu2,iSupMod2,ieta2,iphi2);
- MaxEnergyCellPos(emCells,c2,iSupMod2,ieta2,iphi2);
-
+
p12 = p1+p2;
Float_t invmass = p12.M()*1000;
+ //printf("*** mass %f\n",invmass);
+ Float_t asym = TMath::Abs(p1.E()-p2.E())/(p1.E()+p2.E());
+ //printf("asymmetry %f\n",asym);
+
+ if(asym > fAsyCut) continue;
if(invmass < fMaxBin && invmass > fMinBin){
- fHmgg->Fill(invmass,p12.Pt());
+
+ //Check if cluster is in fidutial region, not too close to borders
+ Bool_t in1 = CheckCellFiducialRegion(c1, aod->GetEMCALCells());
+ Bool_t in2 = CheckCellFiducialRegion(c2, aod->GetEMCALCells());
+
+ if(in1 && in2){
+
+ fHmgg->Fill(invmass,p12.Pt());
+
+ if(iSupMod1==iSupMod2) fHmggSM[iSupMod1]->Fill(invmass,p12.Pt());
+ else fHmggDifferentSM ->Fill(invmass,p12.Pt());
+
+ if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) fHmggPairSM[0]->Fill(invmass,p12.Pt());
+ if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) fHmggPairSM[1]->Fill(invmass,p12.Pt());
+ if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) fHmggPairSM[2]->Fill(invmass,p12.Pt());
+ if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) fHmggPairSM[3]->Fill(invmass,p12.Pt());
+
+ if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+
+ //Opening angle of 2 photons
+ Float_t opangle = p1.Angle(p2.Vect())*TMath::RadToDeg();
+ //printf("*******>>>>>>>> In PEAK pt %f, angle %f \n",p12.Pt(),opangle);
+
+ //Inciden angle of each photon
+ //Float_t * posSM1cen = RecalculatePosition(11.5, 23.5, p1.E(),0, iSupMod1);
+ //Float_t * posSM2cen = RecalculatePosition(11.5, 23.5, p2.E(),0, iSupMod2);
+ Float_t posSM1cen[3]={0.,0.,0.};
+ fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, p1.E(),iSupMod1,0,
+ fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
+ Float_t posSM2cen[3]={0.,0.,0.};
+ fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, p2.E(),iSupMod2,0,
+ fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
+ //printf("SM1 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod1,posSM1cen[0],posSM1cen[1],posSM1cen[2]);
+ //printf("SM2 %d pos (%2.3f,%2.3f,%2.3f) \n",iSupMod2,posSM2cen[0],posSM2cen[1],posSM2cen[2]);
+
+ TVector3 vecSM1cen(posSM1cen[0]-v[0],posSM1cen[1]-v[1],posSM1cen[2]-v[2]);
+ TVector3 vecSM2cen(posSM2cen[0]-v[0],posSM2cen[1]-v[1],posSM2cen[2]-v[2]);
+ Float_t inangle1 = p1.Angle(vecSM1cen)*TMath::RadToDeg();
+ Float_t inangle2 = p2.Angle(vecSM2cen)*TMath::RadToDeg();
+ //printf("Incident angle: cluster 1 %2.3f; cluster 2 %2.3f\n",inangle1,inangle2);
+
+ fHOpeningAngle ->Fill(opangle,p12.Pt());
+ fHIncidentAngle->Fill(inangle1,p1.Pt());
+ fHIncidentAngle->Fill(inangle2,p2.Pt());
+ fHAsymmetry ->Fill(asym,p12.Pt());
+
+ if(iSupMod1==iSupMod2) {
+ fHOpeningAngleSM[iSupMod1] ->Fill(opangle,p12.Pt());
+ fHIncidentAngleSM[iSupMod1]->Fill(inangle1,p1.Pt());
+ fHIncidentAngleSM[iSupMod1]->Fill(inangle2,p2.Pt());
+ fHAsymmetrySM[iSupMod1] ->Fill(asym,p12.Pt());
+ }
+ else{
+ fHOpeningAngleDifferentSM ->Fill(opangle,p12.Pt());
+ fHIncidentAngleDifferentSM ->Fill(inangle1,p1.Pt());
+ fHIncidentAngleDifferentSM ->Fill(inangle2,p2.Pt());
+ fHAsymmetryDifferentSM ->Fill(asym,p12.Pt());
+ }
+
+ if((iSupMod1==0 && iSupMod2==2) || (iSupMod1==2 && iSupMod2==0)) {
+ fHOpeningAnglePairSM[0] ->Fill(opangle,p12.Pt());
+ fHIncidentAnglePairSM[0]->Fill(inangle1,p1.Pt());
+ fHIncidentAnglePairSM[0]->Fill(inangle2,p2.Pt());
+ fHAsymmetryPairSM[0] ->Fill(asym,p12.Pt());
+
+ }
+ if((iSupMod1==1 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==1)) {
+ fHOpeningAnglePairSM[1] ->Fill(opangle,p12.Pt());
+ fHIncidentAnglePairSM[1]->Fill(inangle1,p1.Pt());
+ fHIncidentAnglePairSM[1]->Fill(inangle2,p2.Pt());
+ fHAsymmetryPairSM[1] ->Fill(asym,p12.Pt());
+
+ }
+
+ if((iSupMod1==0 && iSupMod2==1) || (iSupMod1==1 && iSupMod2==0)) {
+ fHOpeningAnglePairSM[2] ->Fill(opangle,p12.Pt());
+ fHIncidentAnglePairSM[2]->Fill(inangle1,p1.Pt());
+ fHIncidentAnglePairSM[2]->Fill(inangle2,p2.Pt());
+ fHAsymmetryPairSM[2] ->Fill(asym,p12.Pt());
+
+
+ }
+ if((iSupMod1==2 && iSupMod2==3) || (iSupMod1==3 && iSupMod2==2)) {
+ fHOpeningAnglePairSM[3] ->Fill(opangle,p12.Pt());
+ fHIncidentAnglePairSM[3]->Fill(inangle1,p1.Pt());
+ fHIncidentAnglePairSM[3]->Fill(inangle2,p2.Pt());
+ fHAsymmetryPairSM[3] ->Fill(asym,p12.Pt());
+ }
+
+ }// pair in 100 < mass < 160
+
+ }//in acceptance cuts
+
+ //In case of filling only channels with second cluster in same SM
+ if(fSameSM && iSupMod1!=iSupMod2) continue;
if (fGroupNCells == 0){
- fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
- fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+ fHmpi0[iSupMod1][ieta1][iphi1]->Fill(invmass);
+ fHmpi0[iSupMod2][ieta2][iphi2]->Fill(invmass);
+
+ if(invmass > 100. && invmass < 160.){//restrict to clusters really close to pi0 peak
+ fhTowerDecayPhotonHit [iSupMod1]->Fill(ieta1,iphi1);
+ fhTowerDecayPhotonEnergy [iSupMod1]->Fill(ieta1,iphi1,p1.E());
+ fhTowerDecayPhotonAsymmetry[iSupMod1]->Fill(ieta1,iphi1,asym);
+
+ fhTowerDecayPhotonHit [iSupMod2]->Fill(ieta2,iphi2);
+ fhTowerDecayPhotonEnergy [iSupMod2]->Fill(ieta2,iphi2,p2.E());
+ fhTowerDecayPhotonAsymmetry[iSupMod2]->Fill(ieta2,iphi2,asym);
+
+ }// pair in mass of pi0
}
else {
//printf("Regroup N %d, eta1 %d, phi1 %d, eta2 %d, phi2 %d \n",fGroupNCells, ieta1, iphi1, ieta2, iphi2);
}
+
+//_______________________________________________________________
+Bool_t AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells)
+{
+
+ // Given the list of AbsId of the cluster, get the maximum cell and
+ // check if there are fNCellsFromBorder from the calorimeter border
+
+ //If the distance to the border is 0 or negative just exit accept all clusters
+ if(cells->GetType()==AliVCaloCells::kEMCALCell && fNCellsFromEMCALBorder <= 0 ) return kTRUE;
+
+ Int_t absIdMax = -1;
+ Float_t ampMax = -1;
+
+ for(Int_t i = 0; i < cluster->GetNCells() ; i++){
+ Int_t absId = cluster->GetCellAbsId(i) ;
+ Float_t amp = cells->GetCellAmplitude(absId);
+ if(amp > ampMax){
+ ampMax = amp;
+ absIdMax = absId;
+ }
+ }
+
+ if(DebugLevel() > 1)
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
+ absIdMax, ampMax, cluster->E());
+
+ if(absIdMax==-1) return kFALSE;
+
+ //Check if the cell is close to the borders:
+ Bool_t okrow = kFALSE;
+ Bool_t okcol = kFALSE;
+
+ Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
+ fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
+ fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
+ if(iSM < 0 || iphi < 0 || ieta < 0 ) {
+ Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
+ }
+
+ //Check rows/phi
+ if(iSM < 10){
+ if(iphi >= fNCellsFromEMCALBorder && iphi < 24-fNCellsFromEMCALBorder) okrow =kTRUE;
+ }
+ else{
+ if(iphi >= fNCellsFromEMCALBorder && iphi < 12-fNCellsFromEMCALBorder) okrow =kTRUE;
+ }
+
+ //Check collumns/eta
+ if(iSM%2==0){
+ if(ieta >= fNCellsFromEMCALBorder) okcol = kTRUE;
+ }
+ else {
+ if(ieta < 48-fNCellsFromEMCALBorder) okcol = kTRUE;
+ }
+
+ if(DebugLevel() > 1)
+ {
+ printf("AliAnalysisTaskEMCALPi0CalibSelection::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
+ fNCellsFromEMCALBorder, ieta, iphi, iSM);
+ if (okcol && okrow ) printf(" YES \n");
+ else printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
+ }
+
+ if (okcol && okrow) return kTRUE;
+ else return kFALSE;
+
+}
+
//__________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::MaxEnergyCellPos(AliAODCaloCells* const cells, AliAODCaloCluster* const clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
+void AliAnalysisTaskEMCALPi0CalibSelection::GetMaxEnergyCellPosAndClusterPos(AliVCaloCells* cells, AliVCluster* clu, Int_t& iSupMod, Int_t& ieta, Int_t& iphi)
{
//For a given CaloCluster calculates the absId of the cell
//with maximum energy deposit.
Int_t iIphi = -1;
Int_t iIeta = -1;
+ Float_t clEnergy = clu->E(); //Energy already recalibrated previously.
+ Float_t weight = 0., weightedCol = 0., weightedRow = 0., totalWeight=0.;
+ Bool_t areInSameSM = kTRUE; //exclude clusters with cells in different SMs for now
+ Int_t startingSM = -1;
+
for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
cellAbsId = clu->GetCellAbsId(iDig);
fraction = clu->GetCellAmplitudeFraction(iDig);
if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
+ Int_t imodrc = -1, iphirc = -1, ietarc =-1;
+ Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
+ fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc);
+ fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);
+ if (iDig==0) startingSM = imodrc;
+ else if(imodrc != startingSM) areInSameSM = kFALSE;
+
if(IsRecalibrationOn()) {
- Int_t imodrc = -1, iphirc = -1, ietarc =-1;
- Int_t iTowerrc = -1, iIphirc = -1, iIetarc =-1;
- fEMCALGeo->GetCellIndex(cellAbsId,imodrc,iTowerrc,iIphirc,iIetarc);
- fEMCALGeo->GetCellPhiEtaIndexInSModule(imodrc,iTowerrc,iIphirc, iIetarc,iphirc,ietarc);
recalFactor = GetEMCALChannelRecalibrationFactor(imodrc,ietarc,iphirc);
}
eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
+
+ weight = TMath::Log(eCell/clEnergy) + 4;
+ if(weight < 0) weight = 0;
+ totalWeight += weight;
+ weightedCol += ietarc*weight;
+ weightedRow += iphirc*weight;
+
//printf("Max cell? cell %d, amplitude org %f, fraction %f, recalibration %f, amplitude new %f \n",cellAbsId, cells->GetCellAmplitude(cellAbsId), fraction, recalFactor, eCell) ;
if(eCell > eMax) {
maxId = cellAbsId;
//printf("\t new max: cell %d, e %f, ecell %f\n",maxId, eMax,eCell);
}
- }
+ }// cell loop
//Get from the absid the supermodule, tower and eta/phi numbers
fEMCALGeo->GetCellIndex(maxId,iSupMod,iTower,iIphi,iIeta);
//Gives SuperModule and Tower numbers
fEMCALGeo->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
- iIphi, iIeta,iphi,ieta);
+ iIphi, iIeta,iphi,ieta);
+ Float_t xyzNew[3];
+ if(areInSameSM == kTRUE) {
+ //printf("In Same SM\n");
+ weightedCol = weightedCol/totalWeight;
+ weightedRow = weightedRow/totalWeight;
+
+ //Float_t *xyzNew = RecalculatePosition(weightedRow, weightedCol, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+ fEMCALGeo->RecalculateTowerPosition(weightedRow, weightedCol, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
+ fRecoUtils->GetMisalTransShiftArray(), fRecoUtils->GetMisalRotShiftArray(), xyzNew);
+ }
+ else {
+ //printf("In Different SM\n");
+ //Float_t *xyzNew = RecalculatePosition(iphi, ieta, clEnergy, 0, iSupMod); //1 = electrons, 0 photons
+ fEMCALGeo->RecalculateTowerPosition(iphi, ieta, iSupMod, clEnergy, 0, //1 = electrons, 0 photons
+ fRecoUtils->GetMisalTransShiftArray(), fRecoUtils->GetMisalRotShiftArray(), xyzNew);
+
+ }
+
+ clu->SetPosition(xyzNew);
+
//printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff