fLogWeight(4.5), fCopyAOD(kFALSE), fSameSM(kFALSE), fOldAOD(kFALSE),
fEMCALGeoName("EMCAL_FIRSTYEAR"), fNCellsFromEMCALBorder(0),
fRemoveBadChannels(kFALSE),fEMCALBadChannelMap(0x0),
- fRecalibration(kFALSE),fEMCALRecalibrationFactors(),
fRecoUtils(new AliEMCALRecoUtils),
fNbins(300), fMinBin(0.), fMaxBin(300.),fOutputContainer(0x0),
fHmgg(0x0), fHmggDifferentSM(0x0),
delete fEMCALBadChannelMap;
}
-
- if(fEMCALRecalibrationFactors) {
- fEMCALRecalibrationFactors->Clear();
- delete fEMCALRecalibrationFactors;
- }
-
if(fRecoUtils) delete fRecoUtils ;
}
snprintf(onePar,buffersize, "Histograms: bins %d; energy range: %2.2f < E < %2.2f GeV;",fNbins,fMinBin,fMaxBin) ;
fCuts->Add(new TObjString(onePar));
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) ;
+ fRemoveBadChannels,fCopyAOD,fRecoUtils->IsRecalibrationOn(), fOldAOD, fSameSM) ;
fCuts->Add(new TObjString(onePar));
snprintf(onePar,buffersize, "EMCAL Geometry name: < %s >",fEMCALGeoName.Data()) ;
fCuts->Add(new TObjString(onePar));
//Analysis per event.
if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection <<< Event %d >>>\n",(Int_t)Entry());
+ if(fRecoUtils->GetParticleType()!=AliEMCALRecoUtils::kPhoton){
+ printf("Wrong particle type for cluster position recalculation! = %d\n", fRecoUtils->GetParticleType());
+ abort();
+ }
+
fhNEvents->Fill(0); //Event analyzed
AliAODEvent* aod = 0x0;
// }
if(DebugLevel() > 1) printf("AliAnalysisTaskEMCALPi0CalibSelection Will use fLogWeight %.3f .\n",fLogWeight);
-
+ Int_t absId1 = -1;
Int_t iSupMod1 = -1;
Int_t iphi1 = -1;
Int_t ieta1 = -1;
+ Int_t absId2 = -1;
Int_t iSupMod2 = -1;
Int_t iphi2 = -1;
Int_t ieta2 = -1;
//clu1.Recalibrate(fCalibData, emCells, fEMCALGeoName);
//clu1.EvalEnergy();
//clu1.EvalAll(fLogWeight, fEMCALGeoName);
- if(IsRecalibrationOn()) {
- Float_t energy = RecalibrateClusterEnergy(c1, emCells);
- //clu1.SetE(energy);
- c1->SetE(energy);
- }
+ if(fRecoUtils->IsRecalibrationOn()) fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c1, emCells);
//Float_t e1ii = clu1.E(); // cluster energy after correction
Float_t e1ii = c1->E(); // cluster energy after correction
-
+
if(DebugLevel() > 2)
- {
- //printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, clu1.GetDispersion(),clu1.GetM02(),clu1.GetM20());
printf("Recal: i %d, E %f, dispersion %f, m02 %f, m20 %f\n",iClu,e1ii, c1->GetDispersion(),c1->GetM02(),c1->GetM20());
- Float_t pos2[]={0,0,0};
- //clu1.GetPosition(pos2);
- c1->GetPosition(pos2);
- printf("Recal: i %d, x %f, y %f, z %f\n",iClu, pos2[0], pos2[1], pos2[2]);
- }
+
//clu1.GetMomentum(p1,v);
//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);
+ fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c1);
+ fRecoUtils->GetMaxEnergyCell (fEMCALGeo, emCells,c1,absId1,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]);
//clu2.Recalibrate(fCalibData, emCells,fEMCALGeoName);
//clu2.EvalEnergy();
//clu2.EvalAll(fLogWeight,fEMCALGeoName);
- if(IsRecalibrationOn()) {
- Float_t energy = RecalibrateClusterEnergy(c2, emCells);
- //clu2.SetE(energy);
- c2->SetE(energy);
- }
-
+ if(fRecoUtils->IsRecalibrationOn()) fRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, c2, emCells);
+
Float_t e2ii = c2->E();
//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);
-
+ fRecoUtils->RecalculateClusterPosition(fEMCALGeo, emCells,c2);
+ fRecoUtils->GetMaxEnergyCell (fEMCALGeo, emCells,c2,absId2,iSupMod2,ieta2,iphi2);
+
c2->GetMomentum(p2,v);
p12 = p1+p2;
//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, iSupMod1, p1.E(), 0,
- fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
+ Float_t depth = fRecoUtils->GetDepth(p1.Energy(),fRecoUtils->GetParticleType(),iSupMod1);
+ fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod1, depth, fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM1cen);
Float_t posSM2cen[3]={0.,0.,0.};
- fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, p2.E(), 0,
- fRecoUtils->GetMisalTransShiftArray(),fRecoUtils->GetMisalRotShiftArray(),posSM2cen);
+ depth = fRecoUtils->GetDepth(p2.Energy(),fRecoUtils->GetParticleType(),iSupMod2);
+ fEMCALGeo->RecalculateTowerPosition(11.5, 23.5, iSupMod2, depth, 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]);
}
-//__________________________________________________
-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.
-
- Double_t eMax = -1.;
- Double_t eCell = -1.;
- Float_t fraction = 1.;
- Int_t cellAbsId = -1;
- Float_t recalFactor = 1.;
-
- Int_t maxId = -1;
- Int_t iTower = -1;
- 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()) {
- 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) {
- eMax = eCell;
- 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);
-
- 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);
-
-}
//__________________________________________________
//void AliAnalysisTaskEMCALPi0CalibSelection::SetCalibCorrections(AliEMCALCalibData* const cdata)
}
-
-//________________________________________________________________
-void AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors(){
- //Init EMCAL recalibration factors
- if(DebugLevel() > 0 )printf("AliAnalysisTaskEMCALPi0CalibSelection::InitEMCALRecalibrationFactors()\n");
- //In order to avoid rewriting the same histograms
- Bool_t oldStatus = TH1::AddDirectoryStatus();
- TH1::AddDirectory(kFALSE);
-
- fEMCALRecalibrationFactors = new TObjArray(12);
- for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
- //Init the histograms with 1
- for (Int_t sm = 0; sm < 12; sm++) {
- for (Int_t i = 0; i < 48; i++) {
- for (Int_t j = 0; j < 24; j++) {
- SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
- }
- }
- }
- fEMCALRecalibrationFactors->SetOwner(kTRUE);
- fEMCALRecalibrationFactors->Compress();
-
- //In order to avoid rewriting the same histograms
- TH1::AddDirectory(oldStatus);
-}
-
-//________________________________________________________________
-Float_t AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AliAODCaloCluster * cluster, AliAODCaloCells * cells){
- // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
- // AOD case
-
- if(!cells) {
- printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy(AOD) - Cells pointer does not exist, stop!");
- abort();
- }
-
- //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
- UShort_t * index = cluster->GetCellsAbsId() ;
- Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
- Int_t ncells = cluster->GetNCells();
-
- //Initialize some used variables
- Float_t energy = 0;
- Int_t absId = -1;
- Int_t icol = -1, irow = -1, imod=1;
- Float_t factor = 1, frac = 0;
-
- //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
- for(Int_t icell = 0; icell < ncells; icell++){
- absId = index[icell];
- frac = fraction[icell];
- if(frac < 1e-5) frac = 1; //in case of EMCAL, this is set as 0 since unfolding is off
- Int_t iTower = -1, iIphi = -1, iIeta = -1;
- fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
- if(fEMCALRecalibrationFactors->GetEntries() <= imod) continue;
- fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
- factor = GetEMCALChannelRecalibrationFactor(imod,icol,irow);
- if(DebugLevel()>2)
- printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - recalibrate cell: module %d, col %d, row %d, cell fraction %f,recalibration factor %f, cell energy %f\n",
- imod,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
-
- energy += cells->GetCellAmplitude(absId)*factor*frac;
- }
-
- if(DebugLevel()>1)
- printf("AliAnalysisTaskEMCALPi0CalibSelection::RecalibrateClusterEnergy - Energy before %f, after %f\n",cluster->E(),energy);
-
- return energy;
-
-}
-
void SetEMCALChannelStatusMap(TObjArray *map) {fEMCALBadChannelMap = map;}
Bool_t ClusterContainsBadChannel(UShort_t* cellList, Int_t nCells);
-
- // Recalibration
- Bool_t IsRecalibrationOn() const { return fRecalibration ; }
- void SwitchOnRecalibration() {fRecalibration = kTRUE ; InitEMCALRecalibrationFactors();}
- void SwitchOffRecalibration() {fRecalibration = kFALSE ; }
-
- void InitEMCALRecalibrationFactors() ;
-
- Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const {
- if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow);
- else return 1;}
-
- void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
- if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
- ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
-
- void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);}
-
- TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);}
-
- void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;}
- Float_t RecalibrateClusterEnergy(AliAODCaloCluster* cluster, AliAODCaloCells * cells);
-
+
void SetEMCALRecoUtils(AliEMCALRecoUtils * ru) {fRecoUtils = ru;}
AliEMCALRecoUtils* GetEMCALRecoUtils() const {return fRecoUtils;}
Bool_t fRemoveBadChannels; // Check the channel status provided and remove clusters with bad channels
TObjArray *fEMCALBadChannelMap; // Array of histograms with map of bad channels, EMCAL
- Bool_t fRecalibration; // Switch on or off the recalibration
- TObjArray *fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
-
AliEMCALRecoUtils * fRecoUtils; // Access to reconstruction utilities
//Output histograms
TH1I* fhNEvents; //! Number of events counter histogram
TList * fCuts ; //! List with analysis cuts
- ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,8);
+ ClassDef(AliAnalysisTaskEMCALPi0CalibSelection,9);
};
//gSystem->Unload("libPWG4CaloCalib.so");
//Try to set the new library
//gSystem->Load("./PWG4CaloCalib/libPWG4CaloCalib.so");
- //gSystem->ListLibraries();
-
+ gSystem->ListLibraries();
+
//-------------------------------------------------------------------------------------------------
//Create chain from ESD and from cross sections files, look below for options.
//-------------------------------------------------------------------------------------------------
cout<<"Wrong data type "<<kInputData<<endl;
break;
}
-
+
TChain *chain = new TChain(kTreeName) ;
CreateChain(mode, chain);
-
+
if(chain){
AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
// Make the analysis manager
//-------------------------------------
AliAnalysisManager *mgr = new AliAnalysisManager("Manager", "Manager");
-
+
// AOD output handler
if(copy){
- AliAODHandler* aodoutHandler = new AliAODHandler();
- aodoutHandler->SetOutputFileName("aod.root");
- ////aodoutHandler->SetCreateNonStandardAOD();
- mgr->SetOutputEventHandler(aodoutHandler);
+ AliAODHandler* aodoutHandler = new AliAODHandler();
+ aodoutHandler->SetOutputFileName("aod.root");
+ ////aodoutHandler->SetCreateNonStandardAOD();
+ mgr->SetOutputEventHandler(aodoutHandler);
}
-
+
//input
if(kInputData == "ESD")
- {
- // ESD handler
- AliESDInputHandler *esdHandler = new AliESDInputHandler();
- mgr->SetInputEventHandler(esdHandler);
- esdHandler->SetReadFriends(kFALSE);
- cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
- }
+ {
+ // ESD handler
+ AliESDInputHandler *esdHandler = new AliESDInputHandler();
+ mgr->SetInputEventHandler(esdHandler);
+ esdHandler->SetReadFriends(kFALSE);
+ cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
+ }
if(kInputData == "AOD")
- {
- // AOD handler
- AliAODInputHandler *aodHandler = new AliAODInputHandler();
- mgr->SetInputEventHandler(aodHandler);
- cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
-
- }
+ {
+ // AOD handler
+ AliAODInputHandler *aodHandler = new AliAODInputHandler();
+ mgr->SetInputEventHandler(aodHandler);
+ cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
+
+ }
- //mgr->SetDebugLevel(1);
+ // mgr->SetDebugLevel(1);
//-------------------------------------------------------------------------
//Define task, put here any other task that you want to use.
//-------------------------------------------------------------------------
AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
-
+
// ESD filter task
if(kInputData == "ESD"){
-
+
gROOT->LoadMacro("AddTaskPhysicsSelection.C");
AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
if(!copy){
- gROOT->LoadMacro("AddTaskESDFilter.C");
- AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
+ gROOT->LoadMacro("AddTaskESDFilter.C");
+ AliAnalysisTaskESDfilter *esdfilter = AddTaskESDFilter(kFALSE);
}
}
-
-
+
+
AliAnalysisTaskEMCALPi0CalibSelection * pi0calib = new AliAnalysisTaskEMCALPi0CalibSelection ("EMCALPi0Calibration");
//pi0calib->SetDebugLevel(10);
pi0calib->CopyAOD(copy);
pi0calib->SetClusterMinEnergy(0.5);
pi0calib->SetClusterMaxEnergy(15.);
- pi0calib->SetAsymmetryCut(0.7);
+ //pi0calib->SetAsymmetryCut(0.7);
pi0calib->SetClusterMinNCells(0);
pi0calib->SetNCellsGroup(0);
pi0calib->SwitchOnBadChannelsRemoval();
pi0calib->SwitchOnOldAODs();
pi0calib->SetNumberOfCellsFromEMCALBorder(1);
AliEMCALRecoUtils * reco = pi0calib->GetEMCALRecoUtils();
- reco->SetMisalTransShift(0,1.08); reco->SetMisalTransShift(1,8.35); reco->SetMisalTransShift(2,1.12);
- reco->SetMisalRotShift(3,-8.05); reco->SetMisalRotShift(4,8.05);
- reco->SetMisalTransShift(3,-0.42); reco->SetMisalTransShift(5,1.55);
+
+ reco->SetParticleType(AliEMCALRecoUtils::kPhoton);
+ reco->SetW0(4.);
+
+ reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
+ reco->SetMisalTransShift(0,1.134); reco->SetMisalTransShift(1,8.2); reco->SetMisalTransShift(2,1.197);
+ reco->SetMisalTransShift(3,-3.093); reco->SetMisalTransShift(5,6.82);reco->SetMisalTransShift(5,1.635);
+
+ //reco->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerIndex);
+ //reco->SetMisalTransShift(0,1.08); reco->SetMisalTransShift(1,8.35); reco->SetMisalTransShift(2,1.12);
+ //reco->SetMisalRotShift(3,-8.05); reco->SetMisalRotShift(4,8.05);
+ //reco->SetMisalTransShift(3,-0.42); reco->SetMisalTransShift(5,1.55);
+
reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaGamma);
reco->SetNonLinearityParam(0,1.04); reco->SetNonLinearityParam(1,-0.1445);
reco->SetNonLinearityParam(2,1.046);
-
- // reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaConversion);
- // reco->SetNonLinearityParam(0,1.033); reco->SetNonLinearityParam(1,0.0566186);
- // reco->SetNonLinearityParam(2,0.982133);
-
-
- // reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
- // reco->SetNonLinearityParam(0,1.001); reco->SetNonLinearityParam(1,-0.01264);
- // reco->SetNonLinearityParam(2,-0.03632);
- // reco->SetNonLinearityParam(3,0.1798); reco->SetNonLinearityParam(4,-0.522);
-
- reco->Print("");
-
+
+// reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0GammaConversion);
+// reco->SetNonLinearityParam(0,1.033); reco->SetNonLinearityParam(1,0.0566186);
+// reco->SetNonLinearityParam(2,0.982133);
+
+
+// reco->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MC);
+// reco->SetNonLinearityParam(0,1.001); reco->SetNonLinearityParam(1,-0.01264);
+// reco->SetNonLinearityParam(2,-0.03632);
+// reco->SetNonLinearityParam(3,0.1798); reco->SetNonLinearityParam(4,-0.522);
+
+ reco->SwitchOnRecalibration();
+ TFile * f = new TFile("RecalibrationFactors.root","read");
+ TH2F * h0 = (TH2F*)f->Get("EMCALRecalFactors_SM0")->Clone();
+ TH2F * h1 = (TH2F*)f->Get("EMCALRecalFactors_SM1")->Clone();
+ TH2F * h2 = (TH2F*)f->Get("EMCALRecalFactors_SM2")->Clone();
+ TH2F * h3 = (TH2F*)f->Get("EMCALRecalFactors_SM3")->Clone();
+
+ reco->SetEMCALChannelRecalibrationFactors(0,h0);
+ reco->SetEMCALChannelRecalibrationFactors(1,h1);
+ reco->SetEMCALChannelRecalibrationFactors(2,h2);
+ reco->SetEMCALChannelRecalibrationFactors(3,h3);
+ reco->Print("");
+
// SM0
pi0calib->SetEMCALChannelStatus(0,3,13); pi0calib->SetEMCALChannelStatus(0,44,1); pi0calib->SetEMCALChannelStatus(0,3,13);
pi0calib->SetEMCALChannelStatus(0,20,7); pi0calib->SetEMCALChannelStatus(0,38,2);
pi0calib->SetEMCALChannelStatus(2,19,22);
//SM3
pi0calib->SetEMCALChannelStatus(3,4,7);
-
+
mgr->AddTask(pi0calib);
- pi0calib->SwitchOnRecalibration();
- TFile * f = new TFile("RecalibrationFactors.root","read");
- TH2F * h0 = (TH2F*)f->Get("EMCALRecalFactors_SM0")->Clone();
- TH2F * h1 = (TH2F*)f->Get("EMCALRecalFactors_SM1")->Clone();
- TH2F * h2 = (TH2F*)f->Get("EMCALRecalFactors_SM2")->Clone();
- TH2F * h3 = (TH2F*)f->Get("EMCALRecalFactors_SM3")->Clone();
-
- pi0calib->SetEMCALChannelRecalibrationFactors(0,h0);
- pi0calib->SetEMCALChannelRecalibrationFactors(1,h1);
- pi0calib->SetEMCALChannelRecalibrationFactors(2,h2);
- pi0calib->SetEMCALChannelRecalibrationFactors(3,h3);
-
// Create containers for input/output
AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
if(copy) AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
AliAnalysisDataContainer *coutput2 =
- mgr->CreateContainer("pi0calib", TList::Class(), AliAnalysisManager::kOutputContainer, "pi0calib.root");
+ mgr->CreateContainer("pi0calib", TList::Class(), AliAnalysisManager::kOutputContainer, "pi0calib.root");
AliAnalysisDataContainer *cout_cuts = mgr->CreateContainer("Cuts", TList::Class(),
- AliAnalysisManager::kOutputContainer, "pi0calib.root");
+ AliAnalysisManager::kOutputContainer, "pi0calib.root");
mgr->ConnectInput (pi0calib, 0, cinput1);
if(copy) mgr->ConnectOutput (pi0calib, 0, coutput1 );
mgr->ConnectOutput (pi0calib, 1, coutput2 );
mgr->ConnectOutput (pi0calib, 2, cout_cuts);
-
+
//-----------------------
// Run the analysis
//-----------------------
mgr->InitAnalysis();
mgr->PrintStatus();
mgr->StartAnalysis(smode.Data(),chain);
-
- cout <<" Analysis ended sucessfully "<< endl ;
-
+
+cout <<" Analysis ended sucessfully "<< endl ;
+
}
else cout << "Chain was not produced ! "<<endl;
gSystem->Load("libANALYSISalice.so");
TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
gSystem->Load("libPWG4CaloCalib.so");
-
+ //SetupPar("PWG4CaloCalib");
/*
// gSystem->Load("libPWG4omega3pi.so");
// gSystem->Load("libCORRFW.so");
printf("AliCaloTrackReader::FillInputEMCAL() - Selected clusters E %3.2f, pt %3.2f, phi %3.2f, eta %3.2f\n",
momentum.E(),momentum.Pt(),momentum.Phi()*TMath::RadToDeg(),momentum.Eta());
- Float_t pos[3];
- clus->GetPosition(pos);
+ //Float_t pos[3];
+ //clus->GetPosition(pos);
//printf("Before Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
//Recalibrate the cluster energy
//Recalculate cluster position
if(GetCaloUtils()->IsRecalculationOfClusterPositionOn()){
GetCaloUtils()->RecalculateClusterPosition(GetEMCALCells(),clus);
- clus->GetPosition(pos);
+ //clus->GetPosition(pos);
//printf("After Corrections: e %f, x %f, y %f, z %f\n",clus->E(),pos[0],pos[1],pos[2]);
}
#include "AliVCluster.h"
#include "AliVCaloCells.h"
#include "AliMixedEvent.h"
-#include "AliEMCALRecoUtils.h"
-
ClassImp(AliCalorimeterUtils)
fEMCALBadChannelMap(0x0),fPHOSBadChannelMap(0x0),
fNCellsFromEMCALBorder(0), fNCellsFromPHOSBorder(0),
fNoEMCALBorderAtEta0(kFALSE),fRecalibration(kFALSE),
- fEMCALRecalibrationFactors(), fPHOSRecalibrationFactors(),
+ fPHOSRecalibrationFactors(),
fEMCALRecoUtils(new AliEMCALRecoUtils),fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE)
{
delete fPHOSBadChannelMap;
}
- if(fEMCALRecalibrationFactors) {
- fEMCALRecalibrationFactors->Clear();
- delete fEMCALRecalibrationFactors;
- }
if(fPHOSRecalibrationFactors) {
fPHOSRecalibrationFactors->Clear();
delete fPHOSRecalibrationFactors;
TH1::AddDirectory(oldStatus);
}
-//________________________________________________________________
-void AliCalorimeterUtils::InitEMCALRecalibrationFactors(){
- //Init EMCAL recalibration factors
- if(fDebug > 0 )printf("AliCalorimeterUtils::InitEMCALRecalibrationFactors()\n");
- //In order to avoid rewriting the same histograms
- Bool_t oldStatus = TH1::AddDirectoryStatus();
- TH1::AddDirectory(kFALSE);
-
- fEMCALRecalibrationFactors = new TObjArray(12);
- for (int i = 0; i < 12; i++) fEMCALRecalibrationFactors->Add(new TH2F(Form("EMCALRecalFactors_SM%d",i),Form("EMCALRecalFactors_SM%d",i), 48, 0, 48, 24, 0, 24));
- //Init the histograms with 1
- for (Int_t sm = 0; sm < 12; sm++) {
- for (Int_t i = 0; i < 48; i++) {
- for (Int_t j = 0; j < 24; j++) {
- SetEMCALChannelRecalibrationFactor(sm,i,j,1.);
- }
- }
- }
- fEMCALRecalibrationFactors->SetOwner(kTRUE);
- fEMCALRecalibrationFactors->Compress();
-
- //In order to avoid rewriting the same histograms
- TH1::AddDirectory(oldStatus);
-}
-
//________________________________________________________________
void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){
//Init EMCAL recalibration factors
//Recalculate EMCAL cluster position
- Double_t eMax = -1.;
- Double_t eCell = -1.;
- Float_t fraction = 1.;
- Int_t cellAbsId = -1;
- Float_t recalFactor = 1.;
-
- Int_t maxId = -1;
- Int_t imod = -1, iphi = -1, ieta =-1;
- Int_t iTower = -1, iIphi = -1, 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
- fEMCALGeo->GetCellIndex(cellAbsId,imod,iTower,iIphi,iIeta);
- fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
- if (iDig==0) startingSM = imod;
- else if(imod != startingSM) areInSameSM = kFALSE;
-
- if(IsRecalibrationOn()) {
- recalFactor = GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
- }
- eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
-
- weight = TMath::Log(eCell/clEnergy) + 4;
- if(weight < 0) weight = 0;
- totalWeight += weight;
- weightedCol += ieta*weight;
- weightedRow += iphi*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) {
- eMax = eCell;
- 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,imod,iTower,iIphi,iIeta);
- //Gives SuperModule and Tower numbers
- fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,
- 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, imod, clEnergy, 0, //1 = electrons, 0 photons
- fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->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, imod, clEnergy, 0, //1 = electrons, 0 photons
- fEMCALRecoUtils->GetMisalTransShiftArray(), fEMCALRecoUtils->GetMisalRotShiftArray(), xyzNew);
-
- }
-
- clu->SetPosition(xyzNew);
-
- //printf("\t Max : cell %d, iSupMod %d, ieta %d, iphi %d \n",maxId,iSupMod, ieta,iphi);
-
+ fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
+
}
class AliVCaloCells;
#include "AliPHOSGeoUtils.h"
#include "AliEMCALGeoUtils.h"
-class AliEMCALRecoUtils;
+#include "AliEMCALRecoUtils.h"
class AliCalorimeterUtils : public TObject {
// Recalibration
Bool_t IsRecalibrationOn() const { return fRecalibration ; }
- void SwitchOnRecalibration() {fRecalibration = kTRUE ; InitEMCALRecalibrationFactors(); InitPHOSRecalibrationFactors();}
- void SwitchOffRecalibration() {fRecalibration = kFALSE ; }
+ void SwitchOnRecalibration() {fRecalibration = kTRUE ; InitPHOSRecalibrationFactors(); fEMCALRecoUtils->SwitchOnRecalibration();}
+ void SwitchOffRecalibration() {fRecalibration = kFALSE;fEMCALRecoUtils->SwitchOffRecalibration();}
- void InitEMCALRecalibrationFactors() ;
void InitPHOSRecalibrationFactors () ;
- Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const {
- if(fEMCALRecalibrationFactors) return (Float_t) ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->GetBinContent(iCol,iRow);
- else return 1;}
+ Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow) const { return fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(iSM , iCol, iRow);}
Float_t GetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow) const {
if(fPHOSRecalibrationFactors)return (Float_t) ((TH2F*)fPHOSRecalibrationFactors->At(imod))->GetBinContent(iCol,iRow);
else return 1;}
void SetEMCALChannelRecalibrationFactor(Int_t iSM , Int_t iCol, Int_t iRow, Double_t c = 1) {
- if(!fEMCALRecalibrationFactors) InitEMCALRecalibrationFactors();
- ((TH2F*)fEMCALRecalibrationFactors->At(iSM))->SetBinContent(iCol,iRow,c);}
+ fEMCALRecoUtils->SetEMCALChannelRecalibrationFactor(iSM,iCol,iRow,c);}
void SetPHOSChannelRecalibrationFactor (Int_t imod, Int_t iCol, Int_t iRow, Double_t c = 1) {
if(!fPHOSRecalibrationFactors) InitPHOSRecalibrationFactors();
((TH2F*)fPHOSRecalibrationFactors->At(imod))->SetBinContent(iCol,iRow,c);}
- void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecalibrationFactors->AddAt(h,iSM);}
+ void SetEMCALChannelRecalibrationFactors(Int_t iSM , TH2F* h) {fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(iSM,h);}
void SetPHOSChannelRecalibrationFactors(Int_t imod , TH2F* h) {fPHOSRecalibrationFactors ->AddAt(h,imod);}
- TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return (TH2F*)fEMCALRecalibrationFactors->At(iSM);}
+ TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const {return fEMCALRecoUtils->GetEMCALChannelRecalibrationFactors(iSM);}
TH2F * GetPHOSChannelRecalibrationFactors(Int_t imod) const {return (TH2F*)fPHOSRecalibrationFactors->At(imod);}
- void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecalibrationFactors = map;}
+ void SetEMCALChannelRecalibrationFactors(TObjArray *map) {fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(map);}
void SetPHOSChannelRecalibrationFactors (TObjArray *map) {fPHOSRecalibrationFactors = map;}
Float_t RecalibrateClusterEnergy(AliVCluster* cluster, AliVCaloCells * cells);
Int_t fNCellsFromPHOSBorder; // Number of cells from PHOS border the cell with maximum amplitude has to be.
Bool_t fNoEMCALBorderAtEta0; // Do fiducial cut in EMCAL region eta = 0?
Bool_t fRecalibration; // Switch on or off the recalibration
- TObjArray * fEMCALRecalibrationFactors; // Array of histograms with map of recalibration factors, EMCAL
TObjArray * fPHOSRecalibrationFactors; // Array of histograms with map of recalibration factors, PHOS
AliEMCALRecoUtils* fEMCALRecoUtils; // EMCAL utils for cluster rereconstruction
Bool_t fRecalculatePosition; // Recalculate cluster position
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff