PWG/CaloTrackCorrBase/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
  16 //_________________________________________________________________________
  17 // Class utility for Calorimeter specific selection methods                ///
  18 //
  19 //
  20 //
  21 //-- Author: Gustavo Conesa (LPSC-Grenoble)
  22 //////////////////////////////////////////////////////////////////////////////
  23
  24
  25 // --- ROOT system ---
  26 #include "TGeoManager.h"
  27
  28 // --- ANALYSIS system ---
  29 #include "AliCalorimeterUtils.h"
  30 #include "AliESDEvent.h"
  31 #include "AliMCEvent.h"
  32 #include "AliStack.h"
  33 #include "AliAODPWG4Particle.h"
  34 #include "AliVCluster.h"
  35 #include "AliVCaloCells.h"
  36 #include "AliMixedEvent.h"
  37
  38 // --- Detector ---
  39 #include "AliEMCALGeometry.h"
  40 #include "AliPHOSGeoUtils.h"
  41
  42 ClassImp(AliCalorimeterUtils)
  43
  44
  45 //____________________________________________
  46   AliCalorimeterUtils::AliCalorimeterUtils() :
  47     TObject(), fDebug(0),
  48     fEMCALGeoName("EMCAL_COMPLETE12SMV1"),
  49     fPHOSGeoName ("PHOSgeo"),
  50     fEMCALGeo(0x0),                   fPHOSGeo(0x0),
  51     fEMCALGeoMatrixSet(kFALSE),       fPHOSGeoMatrixSet(kFALSE),
  52     fLoadEMCALMatrices(kFALSE),       fLoadPHOSMatrices(kFALSE),
  53     fRemoveBadChannels(kFALSE),       fPHOSBadChannelMap(0x0),
  54     fNCellsFromPHOSBorder(0),
  55     fNMaskCellColumns(0),             fMaskCellColumns(0x0),
  56     fRecalibration(kFALSE),           fPHOSRecalibrationFactors(),
  57     fEMCALRecoUtils(new AliEMCALRecoUtils),
  58     fRecalculatePosition(kFALSE),     fCorrectELinearity(kFALSE),
  59     fRecalculateMatching(kFALSE),
  60     fCutR(20),                        fCutZ(20),
  61     fCutEta(20),                      fCutPhi(20),
  62     fLocMaxCutE(0),                   fLocMaxCutEDiff(0)
  63 {
  64   //Ctor
  65
  66   //Initialize parameters
  67   InitParameters();
  68   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
  69   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
  70
  71 }
  72
  73 //_________________________________________
  74 AliCalorimeterUtils::~AliCalorimeterUtils()
  75 {
  76   //Dtor
  77
  78   //if(fPHOSGeo)  delete fPHOSGeo  ;
  79   if(fEMCALGeo) delete fEMCALGeo ;
  80
  81   if(fPHOSBadChannelMap) {
  82     fPHOSBadChannelMap->Clear();
  83     delete  fPHOSBadChannelMap;
  84   }
  85
  86   if(fPHOSRecalibrationFactors) {
  87     fPHOSRecalibrationFactors->Clear();
  88     delete  fPHOSRecalibrationFactors;
  89   }
  90
  91   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
  92   if(fNMaskCellColumns) delete [] fMaskCellColumns;
  93
  94 }
  95
  96 //______________________________________________________________________________________
  97 Bool_t AliCalorimeterUtils::AreNeighbours(const TString calo,
  98                                           const Int_t absId1, const Int_t absId2 ) const
  99 {
 100   // Tells if (true) or not (false) two cells are neighbours
 101   // A neighbour is defined as being two cells which share a side or corner
 102
 103   Bool_t areNeighbours = kFALSE ;
 104
 105   Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
 106   Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
 107
 108   Int_t rowdiff =  0, coldiff =  0;
 109
 110   Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
 111   Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
 112
 113   if(calo=="EMCAL" && nSupMod1!=nSupMod2)
 114   {
 115     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
 116     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 117     if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
 118     else           icol2+=AliEMCALGeoParams::fgkEMCALCols;
 119         }
 120
 121   rowdiff = TMath::Abs( irow1 - irow2 ) ;
 122   coldiff = TMath::Abs( icol1 - icol2 ) ;
 123
 124   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 125     areNeighbours = kTRUE ;
 126
 127   return areNeighbours;
 128
 129 }
 130
 131
 132 //_____________________________________________________________________________________
 133 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
 134                                                     AliVCaloCells* cells,
 135                                                     AliVEvent * event, Int_t iev) const
 136 {
 137
 138         // Given the list of AbsId of the cluster, get the maximum cell and
 139         // check if there are fNCellsFromBorder from the calorimeter border
 140
 141   //If the distance to the border is 0 or negative just exit accept all clusters
 142         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
 143         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 144
 145   Int_t absIdMax        = -1;
 146         Float_t ampMax  = -1;
 147
 148   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 149   Int_t nMixedEvents = 0 ;
 150   Int_t * cellsCumul = NULL ;
 151   Int_t numberOfCells = 0 ;
 152   if (mixEvent){
 153     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 154     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 155       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 156       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 157     }
 158
 159     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 160       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 161       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 162     }
 163
 164     if(cellsCumul){
 165
 166       Int_t startCell = cellsCumul[iev] ;
 167       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 168       //Find cells with maximum amplitude
 169       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 170         Int_t absId = cluster->GetCellAbsId(i) ;
 171         for (Int_t j = startCell; j < endCell ;  j++) {
 172           Short_t cellNumber;
 173           Double_t amp ;
 174           Double_t time;
 175           cells->GetCell(j, cellNumber, amp, time) ;
 176           if (absId == cellNumber) {
 177             if(amp > ampMax){
 178               ampMax   = amp;
 179               absIdMax = absId;
 180             }
 181           }
 182         }
 183       }//loop on cluster cells
 184     }// cells cumul available
 185     else {
 186       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 187       abort();
 188     }
 189   } else {//Normal SE Events
 190     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 191       Int_t absId = cluster->GetCellAbsId(i) ;
 192       Float_t amp       = cells->GetCellAmplitude(absId);
 193       if(amp > ampMax){
 194         ampMax   = amp;
 195         absIdMax = absId;
 196       }
 197     }
 198   }
 199
 200         if(fDebug > 1)
 201                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 202            absIdMax, ampMax, cluster->E());
 203
 204         if(absIdMax==-1) return kFALSE;
 205
 206         //Check if the cell is close to the borders:
 207         Bool_t okrow = kFALSE;
 208         Bool_t okcol = kFALSE;
 209
 210         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 211
 212                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 213                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 214                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 215                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 216       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 217     }
 218
 219                 //Check rows/phi
 220     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 221                 if(iSM < 10){
 222                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 223     }
 224                 else{
 225                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
 226                 }
 227
 228                 //Check columns/eta
 229                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
 230                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 231                 }
 232                 else{
 233                         if(iSM%2==0){
 234                                 if(ieta >= nborder)     okcol = kTRUE;
 235                         }
 236                         else {
 237                                 if(ieta <  48-nborder)  okcol = kTRUE;
 238                         }
 239                 }//eta 0 not checked
 240                 if(fDebug > 1)
 241                 {
 242                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 243              nborder, ieta, iphi, iSM);
 244                         if (okcol && okrow ) printf(" YES \n");
 245                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 246                 }
 247         }//EMCAL
 248         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 249                 Int_t relId[4];
 250                 Int_t irow = -1, icol = -1;
 251                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 252                 irow = relId[2];
 253                 icol = relId[3];
 254                 //imod = relId[0]-1;
 255                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 256                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 257                 if(fDebug > 1)
 258                 {
 259                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 260              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 261                         if (okcol && okrow ) printf(" YES \n");
 262                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 263                 }
 264         }//PHOS
 265
 266         if (okcol && okrow) return kTRUE;
 267         else                return kFALSE;
 268
 269 }
 270
 271 //_________________________________________________________________________________________________________
 272 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
 273 {
 274         // Check that in the cluster cells, there is no bad channel of those stored
 275         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 276
 277         if (!fRemoveBadChannels) return kFALSE;
 278         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 279         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 280         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 281
 282         Int_t icol = -1;
 283         Int_t irow = -1;
 284         Int_t imod = -1;
 285         for(Int_t iCell = 0; iCell<nCells; iCell++){
 286
 287                 //Get the column and row
 288                 if(calorimeter == "EMCAL"){
 289       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 290                 }
 291                 else if(calorimeter=="PHOS"){
 292                         Int_t    relId[4];
 293                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 294                         irow = relId[2];
 295                         icol = relId[3];
 296                         imod = relId[0]-1;
 297                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 298       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
 299                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
 300                 }
 301                 else return kFALSE;
 302
 303         }// cell cluster loop
 304
 305         return kFALSE;
 306
 307 }
 308
 309 //_______________________________________________________________
 310 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
 311 {
 312   // Correct cluster energy non linearity
 313
 314   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 315
 316 }
 317
 318 //________________________________________________________________________________________
 319 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
 320                                              Float_t & clusterFraction) const
 321 {
 322
 323   //For a given CaloCluster gets the absId of the cell
 324   //with maximum energy deposit.
 325
 326   if( !clu || !cells ){
 327     AliInfo("Cluster or cells pointer is null!");
 328     return -1;
 329   }
 330
 331   Double_t eMax        =-1.;
 332   Double_t eTot        = 0.;
 333   Double_t eCell       =-1.;
 334   Float_t  fraction    = 1.;
 335   Float_t  recalFactor = 1.;
 336   Int_t    cellAbsId   =-1 , absId =-1 ;
 337   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
 338
 339   TString           calo = "EMCAL";
 340   if(clu->IsPHOS()) calo = "PHOS";
 341
 342   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 343
 344     cellAbsId = clu->GetCellAbsId(iDig);
 345
 346     fraction  = clu->GetCellAmplitudeFraction(iDig);
 347     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 348
 349     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
 350
 351     if(IsRecalibrationOn()) {
 352       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 353       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
 354     }
 355
 356     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 357
 358     if(eCell > eMax)  {
 359       eMax  = eCell;
 360       absId = cellAbsId;
 361     }
 362
 363     eTot+=eCell;
 364
 365   }// cell loop
 366
 367   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
 368
 369   return absId;
 370
 371 }
 372
 373 //__________________________________________________________________________
 374 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
 375                                                  const AliVEvent* event,
 376                                                  const Int_t index) const
 377 {
 378   // Get the matched track given its index, usually just the first match
 379   // Since it is different for ESDs and AODs here it is a wrap method to do it
 380
 381   AliVTrack *track = 0;
 382
 383   // EMCAL case only when matching is recalculated
 384   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
 385   {
 386     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
 387     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
 388
 389     if(trackIndex < 0 )
 390     {
 391       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
 392     }
 393     else
 394     {
 395       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
 396     }
 397
 398     return track ;
 399
 400   }
 401
 402   // Normal case, get info from ESD or AOD
 403   // ESDs
 404   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 405   {
 406     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 407
 408     if(iESDtrack < 0 )
 409     {
 410       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
 411       return 0x0;
 412     }
 413
 414     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 415
 416   }
 417   else // AODs
 418   {
 419     if(cluster->GetNTracksMatched() > 0 )
 420       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
 421   }
 422
 423   return track ;
 424
 425 }
 426
 427 //_____________________________________________________________________________________________________
 428 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 429 {
 430         //Get the EMCAL/PHOS module number that corresponds to this particle
 431
 432         Int_t absId = -1;
 433         if(particle->GetDetector()=="EMCAL"){
 434                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 435                 if(fDebug > 2)
 436                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 437              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 438                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 439         }//EMCAL
 440         else if(particle->GetDetector()=="PHOS"){
 441     // In case we use the MC reader, the input are TParticles,
 442     // in this case use the corresponing method in PHOS Geometry to get the particle.
 443     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 444     {
 445       Int_t mod =-1;
 446       Double_t z = 0., x=0.;
 447       TParticle* primary = 0x0;
 448       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 449       if(stack) {
 450         primary = stack->Particle(particle->GetCaloLabel(0));
 451       }
 452       else {
 453         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 454       }
 455
 456       if(primary){
 457         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 458       }
 459       else{
 460         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 461       }
 462       return mod;
 463     }
 464     // Input are ESDs or AODs, get the PHOS module number like this.
 465     else{
 466       //FIXME
 467       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 468       //return GetModuleNumber(cluster);
 469       //MEFIX
 470       return -1;
 471     }
 472         }//PHOS
 473
 474         return -1;
 475 }
 476
 477 //_____________________________________________________________________
 478 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 479 {
 480         //Get the EMCAL/PHOS module number that corresponds to this cluster
 481         TLorentzVector lv;
 482         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
 483   if(!cluster){
 484     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 485     return -1;
 486   }
 487         cluster->GetMomentum(lv,v);
 488         Float_t phi = lv.Phi();
 489         if(phi < 0) phi+=TMath::TwoPi();
 490         Int_t absId = -1;
 491         if(cluster->IsEMCAL()){
 492                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
 493                 if(fDebug > 2)
 494                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 495              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 496                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 497         }//EMCAL
 498         else if(cluster->IsPHOS()) {
 499                 Int_t    relId[4];
 500                 if ( cluster->GetNCells() > 0) {
 501                         absId = cluster->GetCellAbsId(0);
 502                         if(fDebug > 2)
 503                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
 504                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
 505                 }
 506                 else return -1;
 507
 508                 if ( absId >= 0) {
 509                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 510                         if(fDebug > 2)
 511                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
 512                         return relId[0]-1;
 513                 }
 514                 else return -1;
 515         }//PHOS
 516
 517         return -1;
 518 }
 519
 520 //___________________________________________________________________________________________________
 521 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
 522                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
 523 {
 524         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 525         Int_t imod = -1;
 526         if ( absId >= 0) {
 527                 if(calo=="EMCAL"){
 528                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 529                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 530                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 531       if(imod < 0 || irow < 0 || icol < 0 ) {
 532         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
 533       }
 534
 535                         //RCU0
 536                         if (0<=irow&&irow<8) iRCU=0; // first cable row
 537                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
 538                         //second cable row
 539                         //RCU1
 540                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
 541                         //second cable row
 542                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
 543
 544                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 545                         if (iRCU<0) {
 546                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
 547                         }
 548
 549                         return imod ;
 550                 }//EMCAL
 551                 else{//PHOS
 552                         Int_t    relId[4];
 553                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 554                         irow = relId[2];
 555                         icol = relId[3];
 556                         imod = relId[0]-1;
 557                         iRCU= (Int_t)(relId[2]-1)/16 ;
 558                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
 559                         if (iRCU >= 4) {
 560                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
 561                         }
 562                         return imod;
 563                 }//PHOS
 564         }
 565
 566         return -1;
 567 }
 568
 569 //___________________________________________________________________________________________
 570 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
 571 {
 572   // Find local maxima in cluster
 573
 574   const Int_t   nc = cluster->GetNCells();
 575   Int_t   *absIdList = new Int_t  [nc];
 576   Float_t *maxEList  = new Float_t[nc];
 577
 578   Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
 579
 580   delete [] absIdList;
 581   delete [] maxEList;
 582
 583   return nMax;
 584
 585 }
 586
 587 //___________________________________________________________________________________________
 588 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
 589                                                   Int_t *absIdList,     Float_t *maxEList)
 590 {
 591   // Find local maxima in cluster
 592
 593   Int_t iDigitN = 0 ;
 594   Int_t iDigit  = 0 ;
 595   Int_t absId1 = -1 ;
 596   Int_t absId2 = -1 ;
 597   const Int_t nCells = cluster->GetNCells();
 598
 599   TString calorimeter = "EMCAL";
 600   if(!cluster->IsEMCAL()) calorimeter = "PHOS";
 601
 602   //printf("cluster : ncells %d \n",nCells);
 603
 604   Float_t emax  = 0;
 605   Int_t   idmax =-1;
 606   for(iDigit = 0; iDigit < nCells ; iDigit++)
 607   {
 608     absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit]  ;
 609     Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 610     RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
 611     if( en > emax )
 612     {
 613       emax  = en ;
 614       idmax = absIdList[iDigit] ;
 615     }
 616     //Int_t icol = -1, irow = -1, iRCU = -1;
 617     //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
 618     //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
 619   }
 620
 621   for(iDigit = 0 ; iDigit < nCells; iDigit++)
 622   {
 623     if(absIdList[iDigit]>=0)
 624     {
 625       absId1 = cluster->GetCellsAbsId()[iDigit];
 626
 627       Float_t en1 = cells->GetCellAmplitude(absId1);
 628       RecalibrateCellAmplitude(en1,calorimeter,absId1);
 629
 630       //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
 631
 632       for(iDigitN = 0; iDigitN < nCells; iDigitN++)
 633       {
 634         absId2 = cluster->GetCellsAbsId()[iDigitN] ;
 635
 636         if(absId2==-1 || absId2==absId1) continue;
 637
 638         //printf("\t %d : absIDj %d\n",iDigitN, absId2);
 639
 640         Float_t en2 = cells->GetCellAmplitude(absId2);
 641         RecalibrateCellAmplitude(en2,calorimeter,absId2);
 642
 643         //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
 644
 645         if ( AreNeighbours(calorimeter, absId1, absId2) )
 646         {
 647           // printf("\t \t Neighbours \n");
 648           if (en1 > en2 )
 649           {
 650             absIdList[iDigitN] = -1 ;
 651             //printf("\t \t indexN %d not local max\n",iDigitN);
 652             // but may be digit too is not local max ?
 653             if(en1 < en2 + fLocMaxCutEDiff) {
 654               //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
 655               absIdList[iDigit] = -1 ;
 656             }
 657           }
 658           else
 659           {
 660             absIdList[iDigit] = -1 ;
 661             //printf("\t \t index %d not local max\n",iDigitN);
 662             // but may be digitN too is not local max ?
 663             if(en1 > en2 - fLocMaxCutEDiff)
 664             {
 665               absIdList[iDigitN] = -1 ;
 666               //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
 667             }
 668           }
 669         } // if Are neighbours
 670         //else printf("\t \t NOT Neighbours \n");
 671       } // while digitN
 672     } // slot not empty
 673   } // while digit
 674
 675   iDigitN = 0 ;
 676   for(iDigit = 0; iDigit < nCells; iDigit++)
 677   {
 678     if(absIdList[iDigit]>=0 )
 679     {
 680       absIdList[iDigitN] = absIdList[iDigit] ;
 681       Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
 682       RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
 683       if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
 684       maxEList[iDigitN] = en ;
 685       //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
 686       iDigitN++ ;
 687     }
 688   }
 689
 690   if(iDigitN == 0)
 691   {
 692     if(fDebug > 0)
 693       printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f\n",
 694              idmax,emax,cluster->E());
 695     iDigitN      = 1     ;
 696     maxEList[0]  = emax  ;
 697     absIdList[0] = idmax ;
 698   }
 699
 700   if(fDebug > 0)
 701   {
 702     printf("AliCalorimeterUtils::GetNumberOfLocalMaxima() - In cluster E %2.2f, M02 %2.2f, M20 %2.2f, N maxima %d \n",
 703            cluster->E(),cluster->GetM02(),cluster->GetM20(), iDigitN);
 704
 705     if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
 706     {
 707       printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
 708     }
 709   }
 710
 711   return iDigitN ;
 712
 713 }
 714
 715 //________________________________________
 716 void AliCalorimeterUtils::InitParameters()
 717 {
 718   //Initialize the parameters of the analysis.
 719
 720   fEMCALGeoName         = "EMCAL_COMPLETEV1";
 721   fPHOSGeoName          = "PHOSgeo";
 722
 723   fEMCALGeoMatrixSet    = kFALSE;
 724   fPHOSGeoMatrixSet     = kFALSE;
 725
 726   fRemoveBadChannels    = kFALSE;
 727
 728   fNCellsFromPHOSBorder = 0;
 729
 730   fLocMaxCutE           = 0.1 ;
 731   fLocMaxCutEDiff       = 0.0 ;
 732
 733   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
 734   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
 735   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
 736   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
 737   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 738   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 739
 740 }
 741
 742
 743 //_____________________________________________________
 744 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
 745 {
 746   //Init PHOS bad channels map
 747   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
 748   //In order to avoid rewriting the same histograms
 749   Bool_t oldStatus = TH1::AddDirectoryStatus();
 750   TH1::AddDirectory(kFALSE);
 751
 752   fPHOSBadChannelMap = new TObjArray(5);
 753   for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),Form("PHOS_BadMap_mod%d",i), 64, 0, 64, 56, 0, 56));
 754
 755   fPHOSBadChannelMap->SetOwner(kTRUE);
 756   fPHOSBadChannelMap->Compress();
 757
 758   //In order to avoid rewriting the same histograms
 759   TH1::AddDirectory(oldStatus);
 760 }
 761
 762 //______________________________________________________
 763 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
 764 {
 765         //Init EMCAL recalibration factors
 766         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
 767         //In order to avoid rewriting the same histograms
 768         Bool_t oldStatus = TH1::AddDirectoryStatus();
 769         TH1::AddDirectory(kFALSE);
 770
 771         fPHOSRecalibrationFactors = new TObjArray(5);
 772         for (int i = 0; i < 5; i++)fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),Form("PHOSRecalFactors_Mod%d",i), 64, 0, 64, 56, 0, 56));
 773         //Init the histograms with 1
 774         for (Int_t m = 0; m < 5; m++) {
 775                 for (Int_t i = 0; i < 56; i++) {
 776                         for (Int_t j = 0; j < 64; j++) {
 777                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
 778                         }
 779                 }
 780         }
 781         fPHOSRecalibrationFactors->SetOwner(kTRUE);
 782         fPHOSRecalibrationFactors->Compress();
 783
 784         //In order to avoid rewriting the same histograms
 785         TH1::AddDirectory(oldStatus);
 786 }
 787
 788
 789 //___________________________________________
 790 void AliCalorimeterUtils::InitEMCALGeometry()
 791 {
 792         //Initialize EMCAL geometry if it did not exist previously
 793         if (!fEMCALGeo){
 794                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 795                 if(fDebug > 0){
 796                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
 797                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 798                         printf("\n");
 799                 }
 800         }
 801 }
 802
 803 //__________________________________________
 804 void AliCalorimeterUtils::InitPHOSGeometry()
 805 {
 806         //Initialize PHOS geometry if it did not exist previously
 807         if (!fPHOSGeo){
 808                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
 809                 if(fDebug > 0){
 810                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
 811                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 812                         printf("\n");
 813                 }
 814         }
 815 }
 816
 817 //_________________________________________________________
 818 void AliCalorimeterUtils::Print(const Option_t * opt) const
 819 {
 820
 821   //Print some relevant parameters set for the analysis
 822   if(! opt)
 823     return;
 824
 825   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 826   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
 827   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
 828          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
 829   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
 830   printf("Recalibrate Clusters? %d\n",fRecalibration);
 831   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
 832   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 833   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
 834
 835   printf("Loc. Max. E > %2.2f\n",       fLocMaxCutE);
 836   printf("Loc. Max. E Diff > %2.2f\n",  fLocMaxCutEDiff);
 837
 838   printf("    \n") ;
 839 }
 840
 841 //__________________________________________________________________
 842 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
 843                                              const Int_t ieta) const
 844 {
 845   //Check if cell is in one of the regions where we have significant amount
 846   //of material in front. Only EMCAL
 847
 848   Int_t icol = ieta;
 849   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 850
 851   if (fNMaskCellColumns && fMaskCellColumns) {
 852     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 853       if(icol==fMaskCellColumns[imask]) return kTRUE;
 854     }
 855   }
 856
 857   return kFALSE;
 858
 859 }
 860
 861 //__________________________________________________________________________________________
 862 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp,
 863                                                    const TString calo, const Int_t id) const
 864 {
 865   //Recaculate cell energy if recalibration factor
 866
 867   Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
 868   Int_t nModule  = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
 869
 870   if (IsRecalibrationOn())
 871   {
 872     if(calo == "PHOS")
 873     {
 874       amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
 875     }
 876     else
 877     {
 878       amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
 879     }
 880   }
 881 }
 882
 883 //_________________________________________________________________________________
 884 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time,
 885                                               const TString calo,
 886                                               const Int_t id, const Int_t bc) const
 887 {
 888   // Recalculate time if time recalibration available for EMCAL
 889   // not ready for PHOS
 890
 891   if(calo == "EMCAL" && GetEMCALRecoUtils()->IsTimeRecalibrationOn())
 892   {
 893     GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
 894   }
 895
 896 }
 897
 898
 899 //__________________________________________________________________________
 900 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
 901                                                       AliVCaloCells * cells)
 902 {
 903         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 904
 905   //Initialize some used variables
 906         Float_t frac  = 0., energy = 0.;
 907
 908         if(cells)
 909   {
 910     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 911
 912     UShort_t * index    = cluster->GetCellsAbsId() ;
 913     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 914
 915     Int_t ncells     = cluster->GetNCells();
 916
 917     TString calo     = "EMCAL";
 918     if(cluster->IsPHOS())
 919       calo = "PHOS";
 920
 921     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 922     for(Int_t icell = 0; icell < ncells; icell++){
 923
 924       Int_t absId = index[icell];
 925
 926       frac =  fraction[icell];
 927       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
 928
 929       Float_t amp = cells->GetCellAmplitude(absId);
 930       RecalibrateCellAmplitude(amp,calo, absId);
 931
 932       if(fDebug>2)
 933         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, cell fraction %f, cell energy %f\n",
 934                calo.Data(),frac,cells->GetCellAmplitude(absId));
 935
 936       energy += amp*frac;
 937     }
 938
 939     if(fDebug>1)
 940       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
 941
 942         }// cells available
 943   else
 944   {
 945     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
 946   }
 947
 948         return energy;
 949 }
 950
 951 //________________________________________________________________________________
 952 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
 953 {
 954   //Set the calorimeters transformation matrices
 955
 956   //Get the EMCAL transformation geometry matrices from ESD
 957   if(!fEMCALGeoMatrixSet && fEMCALGeo){
 958     if(fLoadEMCALMatrices){
 959       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
 960       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 961         if(fEMCALMatrix[mod]){
 962           if(fDebug > 1)
 963             fEMCALMatrix[mod]->Print();
 964           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 965         }
 966       }//SM loop
 967       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 968
 969     }//Load matrices
 970     else if (!gGeoManager) {
 971
 972       if(fDebug > 1)
 973         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
 974       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 975         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 976           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 977           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
 978             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 979           }
 980         }// loop over super modules
 981
 982         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 983
 984       }//ESD as input
 985       else {
 986         if(fDebug > 1)
 987           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 988       }//AOD as input
 989     }//Get matrix from data
 990     else if(gGeoManager){
 991       fEMCALGeoMatrixSet = kTRUE;
 992     }
 993   }//EMCAL geo && no geoManager
 994
 995         //Get the PHOS transformation geometry matrices from ESD
 996   if(!fPHOSGeoMatrixSet && fPHOSGeo){
 997
 998     if(fLoadPHOSMatrices){
 999       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
1000       for(Int_t mod = 0 ; mod < 5 ; mod++){
1001         if(fPHOSMatrix[mod]){
1002           if(fDebug > 1 )
1003             fPHOSMatrix[mod]->Print();
1004           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
1005         }
1006       }//SM loop
1007       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
1008
1009     }//Load matrices
1010     else if (!gGeoManager) {
1011       if(fDebug > 1)
1012         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
1013                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
1014                                 for(Int_t mod = 0; mod < 5; mod++){
1015                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
1016                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
1017                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
1018                                         }
1019                                 }// loop over modules
1020         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
1021                         }//ESD as input
1022                         else {
1023                                 if(fDebug > 1)
1024                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
1025       }//AOD as input
1026     }// get matrix from data
1027     else if(gGeoManager){
1028       fPHOSGeoMatrixSet = kTRUE;
1029     }
1030         }//PHOS geo     and  geoManager was not set
1031
1032 }
1033
1034 //__________________________________________________________________________________________
1035 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1036 {
1037
1038   //Recalculate EMCAL cluster position
1039
1040   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1041
1042 }
1043
1044 //________________________________________________________________________________
1045 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
1046                                                           TObjArray* clusterArray)
1047 {
1048   //Recalculate track matching
1049
1050   if (fRecalculateMatching) {
1051     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ;
1052     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1053     //if(esdevent){
1054     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
1055     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1056     //}
1057   }
1058 }