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