PWG4/PartCorrBase/AliCalorimeterUtils.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id:  $ */
  16
  17 //_________________________________________________________________________
  18 // Class utility for Calorimeter specific selection methods                ///
  19 //
  20 //
  21 //
  22 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25
  26 // --- ROOT system ---
  27 #include "TGeoManager.h"
  28
  29 //---- ANALYSIS system ----
  30 #include "AliCalorimeterUtils.h"
  31 #include "AliESDEvent.h"
  32 #include "AliMCEvent.h"
  33 #include "AliStack.h"
  34 #include "AliAODPWG4Particle.h"
  35 #include "AliVCluster.h"
  36 #include "AliVCaloCells.h"
  37 #include "AliMixedEvent.h"
  38
  39 ClassImp(AliCalorimeterUtils)
  40
  41
  42 //____________________________________________________________________________
  43   AliCalorimeterUtils::AliCalorimeterUtils() :
  44     TObject(), fDebug(0),
  45     fEMCALGeoName("EMCAL_COMPLETEV1"),fPHOSGeoName("PHOSgeo"),
  46     fEMCALGeo(0x0), fPHOSGeo(0x0),
  47     fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
  48     fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
  49     fRemoveBadChannels(kFALSE),fPHOSBadChannelMap(0x0),
  50     fNCellsFromPHOSBorder(0),
  51     fNMaskCellColumns(0), fMaskCellColumns(0x0),
  52     fRecalibration(kFALSE), fPHOSRecalibrationFactors(),
  53     fEMCALRecoUtils(new AliEMCALRecoUtils),
  54     fRecalculatePosition(kFALSE),fCorrectELinearity(kFALSE),
  55     fRecalculateMatching(kFALSE),fCutR(20), fCutZ(20)
  56 {
  57   //Ctor
  58
  59   //Initialize parameters
  60   InitParameters();
  61   for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
  62   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix[i]  = 0 ;
  63
  64 }
  65
  66 //_________________________________
  67 AliCalorimeterUtils::~AliCalorimeterUtils() {
  68   //Dtor
  69
  70   //if(fPHOSGeo)  delete fPHOSGeo  ;
  71   if(fEMCALGeo) delete fEMCALGeo ;
  72
  73   if(fPHOSBadChannelMap) {
  74     fPHOSBadChannelMap->Clear();
  75     delete  fPHOSBadChannelMap;
  76   }
  77
  78   if(fPHOSRecalibrationFactors) {
  79     fPHOSRecalibrationFactors->Clear();
  80     delete  fPHOSRecalibrationFactors;
  81   }
  82
  83   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
  84   if(fNMaskCellColumns) delete [] fMaskCellColumns;
  85
  86 }
  87
  88 //_______________________________________________________________
  89 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster, AliVCaloCells* cells, AliVEvent * event, Int_t iev) const {
  90
  91         // Given the list of AbsId of the cluster, get the maximum cell and
  92         // check if there are fNCellsFromBorder from the calorimeter border
  93
  94   //If the distance to the border is 0 or negative just exit accept all clusters
  95         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
  96         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
  97
  98   Int_t absIdMax        = -1;
  99         Float_t ampMax  = -1;
 100
 101   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 102   Int_t nMixedEvents = 0 ;
 103   Int_t * cellsCumul = NULL ;
 104   Int_t numberOfCells = 0 ;
 105   if (mixEvent){
 106     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 107     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 108       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 109       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 110     }
 111
 112     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 113       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 114       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 115     }
 116
 117     if(cellsCumul){
 118
 119       Int_t startCell = cellsCumul[iev] ;
 120       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 121       //Find cells with maximum amplitude
 122       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 123         Int_t absId = cluster->GetCellAbsId(i) ;
 124         for (Int_t j = startCell; j < endCell ;  j++) {
 125           Short_t cellNumber;
 126           Double_t amp ;
 127           Double_t time;
 128           cells->GetCell(j, cellNumber, amp, time) ;
 129           if (absId == cellNumber) {
 130             if(amp > ampMax){
 131               ampMax   = amp;
 132               absIdMax = absId;
 133             }
 134           }
 135         }
 136       }//loop on cluster cells
 137     }// cells cumul available
 138     else {
 139       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 140       abort();
 141     }
 142   } else {//Normal SE Events
 143     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 144       Int_t absId = cluster->GetCellAbsId(i) ;
 145       Float_t amp       = cells->GetCellAmplitude(absId);
 146       if(amp > ampMax){
 147         ampMax   = amp;
 148         absIdMax = absId;
 149       }
 150     }
 151   }
 152
 153         if(fDebug > 1)
 154                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 155            absIdMax, ampMax, cluster->E());
 156
 157         if(absIdMax==-1) return kFALSE;
 158
 159         //Check if the cell is close to the borders:
 160         Bool_t okrow = kFALSE;
 161         Bool_t okcol = kFALSE;
 162
 163         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 164
 165                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 166                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 167                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 168                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 169       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 170     }
 171
 172                 //Check rows/phi
 173     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 174                 if(iSM < 10){
 175                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 176     }
 177                 else{
 178                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
 179                 }
 180
 181                 //Check columns/eta
 182                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
 183                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 184                 }
 185                 else{
 186                         if(iSM%2==0){
 187                                 if(ieta >= nborder)     okcol = kTRUE;
 188                         }
 189                         else {
 190                                 if(ieta <  48-nborder)  okcol = kTRUE;
 191                         }
 192                 }//eta 0 not checked
 193                 if(fDebug > 1)
 194                 {
 195                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 196              nborder, ieta, iphi, iSM);
 197                         if (okcol && okrow ) printf(" YES \n");
 198                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 199                 }
 200         }//EMCAL
 201         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 202                 Int_t relId[4];
 203                 Int_t irow = -1, icol = -1;
 204                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 205                 irow = relId[2];
 206                 icol = relId[3];
 207                 //imod = relId[0]-1;
 208                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 209                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 210                 if(fDebug > 1)
 211                 {
 212                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 213              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 214                         if (okcol && okrow ) printf(" YES \n");
 215                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 216                 }
 217         }//PHOS
 218
 219         if (okcol && okrow) return kTRUE;
 220         else                return kFALSE;
 221
 222 }
 223
 224 //_________________________________________________________________________________________________________
 225 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells){
 226         // Check that in the cluster cells, there is no bad channel of those stored
 227         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 228
 229         if (!fRemoveBadChannels) return kFALSE;
 230         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 231         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 232         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 233
 234         Int_t icol = -1;
 235         Int_t irow = -1;
 236         Int_t imod = -1;
 237         for(Int_t iCell = 0; iCell<nCells; iCell++){
 238
 239                 //Get the column and row
 240                 if(calorimeter == "EMCAL"){
 241       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 242                 }
 243                 else if(calorimeter=="PHOS"){
 244                         Int_t    relId[4];
 245                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 246                         irow = relId[2];
 247                         icol = relId[3];
 248                         imod = relId[0]-1;
 249                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 250                         if(GetPHOSChannelStatus(imod, icol, irow)) return kTRUE;
 251                 }
 252                 else return kFALSE;
 253
 254         }// cell cluster loop
 255
 256         return kFALSE;
 257
 258 }
 259
 260 //____________________________________________________________________________________________________________________________________________________
 261 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus){
 262   // Correct cluster energy non linearity
 263   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 264 }
 265
 266 //____________________________________________________________________________________________________________________________________________________
 267 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 268 {
 269         //Get the EMCAL/PHOS module number that corresponds to this particle
 270
 271         Int_t absId = -1;
 272         if(particle->GetDetector()=="EMCAL"){
 273                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 274                 if(fDebug > 2)
 275                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 276              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 277                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 278         }//EMCAL
 279         else if(particle->GetDetector()=="PHOS"){
 280     // In case we use the MC reader, the input are TParticles,
 281     // in this case use the corresponing method in PHOS Geometry to get the particle.
 282     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 283     {
 284       Int_t mod =-1;
 285       Double_t z = 0., x=0.;
 286       TParticle* primary = 0x0;
 287       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 288       if(stack) {
 289         primary = stack->Particle(particle->GetCaloLabel(0));
 290       }
 291       else {
 292         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 293       }
 294
 295       if(primary){
 296         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 297       }
 298       else{
 299         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 300       }
 301       return mod;
 302     }
 303     // Input are ESDs or AODs, get the PHOS module number like this.
 304     else{
 305       //FIXME
 306       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 307       //return GetModuleNumber(cluster);
 308       //MEFIX
 309       return -1;
 310     }
 311         }//PHOS
 312
 313         return -1;
 314 }
 315
 316 //____________________________________________________________________________________________________________________________________________________
 317 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 318 {
 319         //Get the EMCAL/PHOS module number that corresponds to this cluster
 320         TLorentzVector lv;
 321         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
 322   if(!cluster){
 323     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 324     return -1;
 325   }
 326         cluster->GetMomentum(lv,v);
 327         Float_t phi = lv.Phi();
 328         if(phi < 0) phi+=TMath::TwoPi();
 329         Int_t absId = -1;
 330         if(cluster->IsEMCAL()){
 331                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
 332                 if(fDebug > 2)
 333                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 334              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 335                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 336         }//EMCAL
 337         else if(cluster->IsPHOS()) {
 338                 Int_t    relId[4];
 339                 if ( cluster->GetNCells() > 0) {
 340                         absId = cluster->GetCellAbsId(0);
 341                         if(fDebug > 2)
 342                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
 343                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
 344                 }
 345                 else return -1;
 346
 347                 if ( absId >= 0) {
 348                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 349                         if(fDebug > 2)
 350                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
 351                         return relId[0]-1;
 352                 }
 353                 else return -1;
 354         }//PHOS
 355
 356         return -1;
 357 }
 358
 359 //_____________________________________________________________________________________________________________
 360 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo, Int_t & icol, Int_t & irow, Int_t & iRCU) const
 361 {
 362         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 363         Int_t imod = -1;
 364         if ( absId >= 0) {
 365                 if(calo=="EMCAL"){
 366                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 367                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 368                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 369       if(imod < 0 || irow < 0 || icol < 0 ) {
 370         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
 371       }
 372
 373                         //RCU0
 374                         if (0<=irow&&irow<8) iRCU=0; // first cable row
 375                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
 376                         //second cable row
 377                         //RCU1
 378                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
 379                         //second cable row
 380                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
 381
 382                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 383                         if (iRCU<0) {
 384                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
 385                         }
 386
 387                         return imod ;
 388                 }//EMCAL
 389                 else{//PHOS
 390                         Int_t    relId[4];
 391                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 392                         irow = relId[2];
 393                         icol = relId[3];
 394                         imod = relId[0]-1;
 395                         iRCU= (Int_t)(relId[2]-1)/16 ;
 396                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
 397                         if (iRCU >= 4) {
 398                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
 399                         }
 400                         return imod;
 401                 }//PHOS
 402         }
 403
 404         return -1;
 405 }
 406
 407 //_______________________________________________________________
 408 //void AliCalorimeterUtils::Init()
 409 //{
 410 //      //Init reader. Method to be called in AliAnaPartCorrMaker
 411 //
 412 //      fEMCALBadChannelMap->SetName(Form("EMCALBadMap_%s",fTaskName.Data()));
 413 //      fPHOSBadChannelMap->SetName(Form("PHOSBadMap_%s",fTaskName.Data()));
 414 //}
 415
 416 //_______________________________________________________________
 417 void AliCalorimeterUtils::InitParameters()
 418 {
 419   //Initialize the parameters of the analysis.
 420   fEMCALGeoName         = "EMCAL_COMPLETEV1";
 421   fPHOSGeoName          = "PHOSgeo";
 422   fEMCALGeoMatrixSet    = kFALSE;
 423   fPHOSGeoMatrixSet     = kFALSE;
 424   fRemoveBadChannels    = kFALSE;
 425   fNCellsFromPHOSBorder = 0;
 426
 427   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
 428   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
 429   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
 430   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
 431   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 432   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 433
 434 }
 435
 436
 437 //________________________________________________________________
 438 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap(){
 439   //Init PHOS bad channels map
 440   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
 441   //In order to avoid rewriting the same histograms
 442   Bool_t oldStatus = TH1::AddDirectoryStatus();
 443   TH1::AddDirectory(kFALSE);
 444
 445   fPHOSBadChannelMap = new TObjArray(5);
 446   for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOSBadChannelMap_Mod%d",i),Form("PHOSBadChannelMap_Mod%d",i), 56, 0, 56, 64, 0, 64));
 447
 448   fPHOSBadChannelMap->SetOwner(kTRUE);
 449   fPHOSBadChannelMap->Compress();
 450
 451   //In order to avoid rewriting the same histograms
 452   TH1::AddDirectory(oldStatus);
 453 }
 454
 455 //________________________________________________________________
 456 void AliCalorimeterUtils::InitPHOSRecalibrationFactors(){
 457         //Init EMCAL recalibration factors
 458         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
 459         //In order to avoid rewriting the same histograms
 460         Bool_t oldStatus = TH1::AddDirectoryStatus();
 461         TH1::AddDirectory(kFALSE);
 462
 463         fPHOSRecalibrationFactors = new TObjArray(5);
 464         for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 56, 0, 56, 64, 0, 64));
 465         //Init the histograms with 1
 466         for (Int_t m = 0; m < 5; m++) {
 467                 for (Int_t i = 0; i < 56; i++) {
 468                         for (Int_t j = 0; j < 64; j++) {
 469                                 SetPHOSChannelRecalibrationFactor(m,i,j,1.);
 470                         }
 471                 }
 472         }
 473         fPHOSRecalibrationFactors->SetOwner(kTRUE);
 474         fPHOSRecalibrationFactors->Compress();
 475
 476         //In order to avoid rewriting the same histograms
 477         TH1::AddDirectory(oldStatus);
 478 }
 479
 480
 481 //________________________________________________________________
 482 void AliCalorimeterUtils::InitEMCALGeometry()
 483 {
 484         //Initialize EMCAL geometry if it did not exist previously
 485         if (!fEMCALGeo){
 486                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 487                 if(fDebug > 0){
 488                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
 489                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 490                         printf("\n");
 491                 }
 492         }
 493 }
 494
 495 //________________________________________________________________
 496 void AliCalorimeterUtils::InitPHOSGeometry()
 497 {
 498         //Initialize PHOS geometry if it did not exist previously
 499         if (!fPHOSGeo){
 500                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
 501                 if(fDebug > 0){
 502                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
 503                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 504                         printf("\n");
 505                 }
 506         }
 507 }
 508
 509 //________________________________________________________________
 510 void AliCalorimeterUtils::Print(const Option_t * opt) const
 511 {
 512
 513   //Print some relevant parameters set for the analysis
 514   if(! opt)
 515     return;
 516
 517   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 518   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
 519   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
 520          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
 521   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
 522   printf("Recalibrate Clusters? %d\n",fRecalibration);
 523   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
 524   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 525   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
 526
 527   printf("    \n") ;
 528 }
 529
 530 //__________________________________________________
 531 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,  const Int_t ieta) const {
 532   //Check if cell is in one of the regions where we have significant amount
 533   //of material in front. Only EMCAL
 534
 535   Int_t icol = ieta;
 536   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 537
 538   if (fNMaskCellColumns && fMaskCellColumns) {
 539     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 540       if(icol==fMaskCellColumns[imask]) return kTRUE;
 541     }
 542   }
 543
 544   return kFALSE;
 545
 546 }
 547
 548
 549 //________________________________________________________________
 550 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster, AliVCaloCells * cells){
 551         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 552
 553   //Initialize some used variables
 554         Float_t energy   = 0;
 555         Int_t absId      = -1;
 556         Int_t icol = -1, irow = -1, iRCU = -1, module=1;
 557         Float_t factor = 1, frac = 0;
 558
 559         if(cells) {
 560
 561     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 562     UShort_t * index    = cluster->GetCellsAbsId() ;
 563     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 564     Int_t ncells     = cluster->GetNCells();
 565     TString calo     = "EMCAL";
 566     if(cluster->IsPHOS()) calo = "PHOS";
 567
 568     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 569     for(Int_t icell = 0; icell < ncells; icell++){
 570       absId = index[icell];
 571       frac =  fraction[icell];
 572       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
 573       module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
 574       if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,icol,irow);
 575       else                  factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
 576       if(fDebug>2)
 577         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n",
 578                calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
 579
 580       energy += cells->GetCellAmplitude(absId)*factor*frac;
 581     }
 582
 583     if(fDebug>1)
 584       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
 585
 586         }// cells available
 587   else{
 588     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
 589   }
 590
 591         return energy;
 592 }
 593
 594 //________________________________________________________________
 595 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
 596 {
 597   //Set the calorimeters transformation matrices
 598
 599   //Get the EMCAL transformation geometry matrices from ESD
 600   if(!fEMCALGeoMatrixSet && fEMCALGeo){
 601     if(fLoadEMCALMatrices){
 602       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
 603       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 604         if(fEMCALMatrix[mod]){
 605           if(fDebug > 1)
 606             fEMCALMatrix[mod]->Print();
 607           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 608         }
 609       }//SM loop
 610       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 611
 612     }//Load matrices
 613     else if (!gGeoManager) {
 614
 615       if(fDebug > 1)
 616         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
 617       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 618         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 619           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 620           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
 621             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 622           }
 623         }// loop over super modules
 624         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 625
 626       }//ESD as input
 627       else {
 628         if(fDebug > 1)
 629           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 630       }//AOD as input
 631     }//Get matrix from data
 632     else if(gGeoManager){
 633       fEMCALGeoMatrixSet = kTRUE;
 634     }
 635   }//EMCAL geo && no geoManager
 636
 637         //Get the PHOS transformation geometry matrices from ESD
 638   if(!fPHOSGeoMatrixSet && fPHOSGeo){
 639     if(fLoadPHOSMatrices){
 640       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined geometry matrices\n");
 641       for(Int_t mod=0; mod < 5; mod++){
 642         if(fPHOSMatrix[mod]){
 643           if(fDebug > 1)
 644             fPHOSMatrix[mod]->Print();
 645           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
 646         }
 647       }//SM loop
 648       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
 649     }//Load matrices
 650     else if (!gGeoManager) {
 651       if(fDebug > 1)
 652         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
 653                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 654                                 for(Int_t mod=0; mod < 5; mod++){
 655                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
 656                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
 657                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
 658                                         }
 659                                 }// loop over modules
 660         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
 661                         }//ESD as input
 662                         else {
 663                                 if(fDebug > 1)
 664                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 665       }//AOD as input
 666     }// get matrix from data
 667     else if(gGeoManager){
 668       fPHOSGeoMatrixSet = kTRUE;
 669     }
 670         }//PHOS geo     and  geoManager was not set
 671 }
 672
 673 //________________________________________________________________
 674 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu){
 675
 676   //Recalculate EMCAL cluster position
 677
 678   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
 679
 680 }