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 {
  63   //Ctor
  64
  65   //Initialize parameters
  66   InitParameters();
  67   for(Int_t i = 0; i < 12; i++) fEMCALMatrix[i] = 0 ;
  68   for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
  69
  70 }
  71
  72 //_________________________________________
  73 AliCalorimeterUtils::~AliCalorimeterUtils()
  74 {
  75   //Dtor
  76
  77   //if(fPHOSGeo)  delete fPHOSGeo  ;
  78   if(fEMCALGeo) delete fEMCALGeo ;
  79
  80   if(fPHOSBadChannelMap) {
  81     fPHOSBadChannelMap->Clear();
  82     delete  fPHOSBadChannelMap;
  83   }
  84
  85   if(fPHOSRecalibrationFactors) {
  86     fPHOSRecalibrationFactors->Clear();
  87     delete  fPHOSRecalibrationFactors;
  88   }
  89
  90   if(fEMCALRecoUtils)   delete fEMCALRecoUtils ;
  91   if(fNMaskCellColumns) delete [] fMaskCellColumns;
  92
  93 }
  94
  95 //_____________________________________________________________________________________
  96 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
  97                                                     AliVCaloCells* cells,
  98                                                     AliVEvent * event, Int_t iev) const
  99 {
 100
 101         // Given the list of AbsId of the cluster, get the maximum cell and
 102         // check if there are fNCellsFromBorder from the calorimeter border
 103
 104   //If the distance to the border is 0 or negative just exit accept all clusters
 105         if(cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
 106         if(cells->GetType()==AliVCaloCells::kPHOSCell  && fNCellsFromPHOSBorder  <= 0 ) return kTRUE;
 107
 108   Int_t absIdMax        = -1;
 109         Float_t ampMax  = -1;
 110
 111   AliMixedEvent * mixEvent = dynamic_cast<AliMixedEvent*> (event);
 112   Int_t nMixedEvents = 0 ;
 113   Int_t * cellsCumul = NULL ;
 114   Int_t numberOfCells = 0 ;
 115   if (mixEvent){
 116     nMixedEvents = mixEvent->GetNumberOfEvents() ;
 117     if (cells->GetType()==AliVCaloCells::kEMCALCell) {
 118       cellsCumul =  mixEvent->GetEMCALCellsCumul() ;
 119       numberOfCells = mixEvent->GetNumberOfEMCALCells() ;
 120     }
 121
 122     else if (cells->GetType()==AliVCaloCells::kPHOSCell) {
 123       cellsCumul =  mixEvent->GetPHOSCellsCumul() ;
 124       numberOfCells = mixEvent->GetNumberOfPHOSCells() ;
 125     }
 126
 127     if(cellsCumul){
 128
 129       Int_t startCell = cellsCumul[iev] ;
 130       Int_t endCell   = (iev+1 < nMixedEvents)?cellsCumul[iev+1]:numberOfCells;
 131       //Find cells with maximum amplitude
 132       for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 133         Int_t absId = cluster->GetCellAbsId(i) ;
 134         for (Int_t j = startCell; j < endCell ;  j++) {
 135           Short_t cellNumber;
 136           Double_t amp ;
 137           Double_t time;
 138           cells->GetCell(j, cellNumber, amp, time) ;
 139           if (absId == cellNumber) {
 140             if(amp > ampMax){
 141               ampMax   = amp;
 142               absIdMax = absId;
 143             }
 144           }
 145         }
 146       }//loop on cluster cells
 147     }// cells cumul available
 148     else {
 149       printf("AliCalorimeterUtils::CheckCellFiducialRegion() - CellsCumul is NULL!!!\n");
 150       abort();
 151     }
 152   } else {//Normal SE Events
 153     for(Int_t i = 0; i < cluster->GetNCells() ; i++){
 154       Int_t absId = cluster->GetCellAbsId(i) ;
 155       Float_t amp       = cells->GetCellAmplitude(absId);
 156       if(amp > ampMax){
 157         ampMax   = amp;
 158         absIdMax = absId;
 159       }
 160     }
 161   }
 162
 163         if(fDebug > 1)
 164                 printf("AliCalorimeterUtils::CheckCellFiducialRegion() - Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f\n",
 165            absIdMax, ampMax, cluster->E());
 166
 167         if(absIdMax==-1) return kFALSE;
 168
 169         //Check if the cell is close to the borders:
 170         Bool_t okrow = kFALSE;
 171         Bool_t okcol = kFALSE;
 172
 173         if(cells->GetType()==AliVCaloCells::kEMCALCell){
 174
 175                 Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
 176                 fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
 177                 fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
 178                 if(iSM < 0 || iphi < 0 || ieta < 0 ) {
 179       Fatal("CheckCellFidutialRegion","Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name\n",iSM,ieta,iphi);
 180     }
 181
 182                 //Check rows/phi
 183     Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
 184                 if(iSM < 10){
 185                         if(iphi >= nborder && iphi < 24-nborder) okrow =kTRUE;
 186     }
 187                 else{
 188                         if(iphi >= nborder && iphi < 12-nborder) okrow =kTRUE;
 189                 }
 190
 191                 //Check columns/eta
 192                 if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0()){
 193                         if(ieta  > nborder && ieta < 48-nborder) okcol =kTRUE;
 194                 }
 195                 else{
 196                         if(iSM%2==0){
 197                                 if(ieta >= nborder)     okcol = kTRUE;
 198                         }
 199                         else {
 200                                 if(ieta <  48-nborder)  okcol = kTRUE;
 201                         }
 202                 }//eta 0 not checked
 203                 if(fDebug > 1)
 204                 {
 205                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ?",
 206              nborder, ieta, iphi, iSM);
 207                         if (okcol && okrow ) printf(" YES \n");
 208                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 209                 }
 210         }//EMCAL
 211         else if(cells->GetType()==AliVCaloCells::kPHOSCell){
 212                 Int_t relId[4];
 213                 Int_t irow = -1, icol = -1;
 214                 fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
 215                 irow = relId[2];
 216                 icol = relId[3];
 217                 //imod = relId[0]-1;
 218                 if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
 219                 if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
 220                 if(fDebug > 1)
 221                 {
 222                         printf("AliCalorimeterUtils::CheckCellFiducialRegion() - PHOS Cluster in %d cells fiducial volume: icol %d, irow %d, Module %d?",
 223              fNCellsFromPHOSBorder, icol, irow, relId[0]-1);
 224                         if (okcol && okrow ) printf(" YES \n");
 225                         else  printf(" NO: column ok? %d, row ok? %d \n",okcol,okrow);
 226                 }
 227         }//PHOS
 228
 229         if (okcol && okrow) return kTRUE;
 230         else                return kFALSE;
 231
 232 }
 233
 234 //_________________________________________________________________________________________________________
 235 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(TString calorimeter,UShort_t* cellList, Int_t nCells)
 236 {
 237         // Check that in the cluster cells, there is no bad channel of those stored
 238         // in fEMCALBadChannelMap or fPHOSBadChannelMap
 239
 240         if (!fRemoveBadChannels) return kFALSE;
 241         //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
 242         if(calorimeter == "EMCAL" && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
 243         if(calorimeter == "PHOS"  && !fPHOSBadChannelMap)  return kFALSE;
 244
 245         Int_t icol = -1;
 246         Int_t irow = -1;
 247         Int_t imod = -1;
 248         for(Int_t iCell = 0; iCell<nCells; iCell++){
 249
 250                 //Get the column and row
 251                 if(calorimeter == "EMCAL"){
 252       return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
 253                 }
 254                 else if(calorimeter=="PHOS"){
 255                         Int_t    relId[4];
 256                         fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
 257                         irow = relId[2];
 258                         icol = relId[3];
 259                         imod = relId[0]-1;
 260                         if(fPHOSBadChannelMap->GetEntries() <= imod)continue;
 261       //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
 262                         if(GetPHOSChannelStatus(imod, irow, icol)) return kTRUE;
 263                 }
 264                 else return kFALSE;
 265
 266         }// cell cluster loop
 267
 268         return kFALSE;
 269
 270 }
 271
 272 //_______________________________________________________________
 273 void AliCalorimeterUtils::CorrectClusterEnergy(AliVCluster *clus)
 274 {
 275   // Correct cluster energy non linearity
 276
 277   clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
 278
 279 }
 280
 281 //________________________________________________________________________________________
 282 Int_t  AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
 283                                              Float_t & clusterFraction) const
 284 {
 285
 286   //For a given CaloCluster gets the absId of the cell
 287   //with maximum energy deposit.
 288
 289   if( !clu || !cells ){
 290     AliInfo("Cluster or cells pointer is null!");
 291     return -1;
 292   }
 293
 294   Double_t eMax        =-1.;
 295   Double_t eTot        = 0.;
 296   Double_t eCell       =-1.;
 297   Float_t  fraction    = 1.;
 298   Float_t  recalFactor = 1.;
 299   Int_t    cellAbsId   =-1 , absId =-1 ;
 300   Int_t    iSupMod     =-1 , ieta  =-1 , iphi = -1, iRCU = -1;
 301
 302   TString           calo = "EMCAL";
 303   if(clu->IsPHOS()) calo = "PHOS";
 304
 305   for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++) {
 306
 307     cellAbsId = clu->GetCellAbsId(iDig);
 308
 309     fraction  = clu->GetCellAmplitudeFraction(iDig);
 310     if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
 311
 312     iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
 313
 314     if(IsRecalibrationOn()) {
 315       if(calo=="EMCAL") recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
 316       else              recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
 317     }
 318
 319     eCell  = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
 320
 321     if(eCell > eMax)  {
 322       eMax  = eCell;
 323       absId = cellAbsId;
 324     }
 325
 326     eTot+=eCell;
 327
 328   }// cell loop
 329
 330   clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
 331
 332   return absId;
 333
 334 }
 335
 336 //__________________________________________________________________________
 337 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(const AliVCluster* cluster,
 338                                                  const AliVEvent* event,
 339                                                  const Int_t index) const
 340 {
 341   // Get the matched track given its index, usually just the first match
 342   // Since it is different for ESDs and AODs here it is a wrap method to do it
 343
 344   AliVTrack *track = 0;
 345
 346   // EMCAL case only when matching is recalculated
 347   if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
 348   {
 349     Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
 350     //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
 351
 352     if(trackIndex < 0 )
 353     {
 354       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d, from recalculation\n", trackIndex);
 355     }
 356     else
 357     {
 358       track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
 359     }
 360
 361     return track ;
 362
 363   }
 364
 365   // Normal case, get info from ESD or AOD
 366   // ESDs
 367   if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 368   {
 369     Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 370
 371     if(iESDtrack < 0 )
 372     {
 373       printf("AliCalorimeterUtils::GetMatchedTrack() - Wrong track index %d\n", index);
 374       return 0x0;
 375     }
 376
 377     track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 378
 379   }
 380   else // AODs
 381   {
 382     if(cluster->GetNTracksMatched() > 0 )
 383       track = dynamic_cast<AliVTrack*>(cluster->GetTrackMatched(index));
 384   }
 385
 386   return track ;
 387
 388 }
 389
 390 //_____________________________________________________________________________________________________
 391 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
 392 {
 393         //Get the EMCAL/PHOS module number that corresponds to this particle
 394
 395         Int_t absId = -1;
 396         if(particle->GetDetector()=="EMCAL"){
 397                 fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
 398                 if(fDebug > 2)
 399                   printf("AliCalorimeterUtils::GetModuleNumber(PWG4AOD) - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 400              particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 401                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 402         }//EMCAL
 403         else if(particle->GetDetector()=="PHOS"){
 404     // In case we use the MC reader, the input are TParticles,
 405     // in this case use the corresponing method in PHOS Geometry to get the particle.
 406     if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
 407     {
 408       Int_t mod =-1;
 409       Double_t z = 0., x=0.;
 410       TParticle* primary = 0x0;
 411       AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
 412       if(stack) {
 413         primary = stack->Particle(particle->GetCaloLabel(0));
 414       }
 415       else {
 416         Fatal("GetModuleNumber(PWG4AOD)", "Stack not available, stop!");
 417       }
 418
 419       if(primary){
 420         fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
 421       }
 422       else{
 423         Fatal("GetModuleNumber(PWG4AOD)", "Primary not available, stop!");
 424       }
 425       return mod;
 426     }
 427     // Input are ESDs or AODs, get the PHOS module number like this.
 428     else{
 429       //FIXME
 430       //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
 431       //return GetModuleNumber(cluster);
 432       //MEFIX
 433       return -1;
 434     }
 435         }//PHOS
 436
 437         return -1;
 438 }
 439
 440 //_____________________________________________________________________
 441 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
 442 {
 443         //Get the EMCAL/PHOS module number that corresponds to this cluster
 444         TLorentzVector lv;
 445         Double_t v[]={0.,0.,0.}; //not necessary to pass the real vertex.
 446   if(!cluster){
 447     if(fDebug > 1) printf("AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
 448     return -1;
 449   }
 450         cluster->GetMomentum(lv,v);
 451         Float_t phi = lv.Phi();
 452         if(phi < 0) phi+=TMath::TwoPi();
 453         Int_t absId = -1;
 454         if(cluster->IsEMCAL()){
 455                 fEMCALGeo->GetAbsCellIdFromEtaPhi(lv.Eta(),phi, absId);
 456                 if(fDebug > 2)
 457                   printf("AliCalorimeterUtils::GetModuleNumber() - EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d\n",
 458              lv.Eta(), phi*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId));
 459                 return fEMCALGeo->GetSuperModuleNumber(absId) ;
 460         }//EMCAL
 461         else if(cluster->IsPHOS()) {
 462                 Int_t    relId[4];
 463                 if ( cluster->GetNCells() > 0) {
 464                         absId = cluster->GetCellAbsId(0);
 465                         if(fDebug > 2)
 466                                 printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: cluster eta %f, phi %f, e %f, absId %d\n",
 467                lv.Eta(), phi*TMath::RadToDeg(), lv.E(), absId);
 468                 }
 469                 else return -1;
 470
 471                 if ( absId >= 0) {
 472                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 473                         if(fDebug > 2)
 474                           printf("AliCalorimeterUtils::GetModuleNumber() - PHOS: Module %d\n",relId[0]-1);
 475                         return relId[0]-1;
 476                 }
 477                 else return -1;
 478         }//PHOS
 479
 480         return -1;
 481 }
 482
 483 //___________________________________________________________________________________________________
 484 Int_t AliCalorimeterUtils::GetModuleNumberCellIndexes(const Int_t absId, const TString calo,
 485                                                       Int_t & icol, Int_t & irow, Int_t & iRCU) const
 486 {
 487         //Get the EMCAL/PHOS module, columns, row and RCU number that corresponds to this absId
 488         Int_t imod = -1;
 489         if ( absId >= 0) {
 490                 if(calo=="EMCAL"){
 491                         Int_t iTower = -1, iIphi = -1, iIeta = -1;
 492                         fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
 493                         fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
 494       if(imod < 0 || irow < 0 || icol < 0 ) {
 495         Fatal("GetModuleNumberCellIndexes()","Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name\n",imod,icol,irow);
 496       }
 497
 498                         //RCU0
 499                         if (0<=irow&&irow<8) iRCU=0; // first cable row
 500                         else if (8<=irow&&irow<16 && 0<=icol&&icol<24) iRCU=0; // first half;
 501                         //second cable row
 502                         //RCU1
 503                         else if(8<=irow&&irow<16 && 24<=icol&&icol<48) iRCU=1; // second half;
 504                         //second cable row
 505                         else if(16<=irow&&irow<24) iRCU=1; // third cable row
 506
 507                         if (imod%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
 508                         if (iRCU<0) {
 509                                 Fatal("GetModuleNumberCellIndexes()","Wrong EMCAL RCU number = %d\n", iRCU);
 510                         }
 511
 512                         return imod ;
 513                 }//EMCAL
 514                 else{//PHOS
 515                         Int_t    relId[4];
 516                         fPHOSGeo->AbsToRelNumbering(absId,relId);
 517                         irow = relId[2];
 518                         icol = relId[3];
 519                         imod = relId[0]-1;
 520                         iRCU= (Int_t)(relId[2]-1)/16 ;
 521                         //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
 522                         if (iRCU >= 4) {
 523                                 Fatal("GetModuleNumberCellIndexes()","Wrong PHOS RCU number = %d\n", iRCU);
 524                         }
 525                         return imod;
 526                 }//PHOS
 527         }
 528
 529         return -1;
 530 }
 531
 532 //________________________________________
 533 void AliCalorimeterUtils::InitParameters()
 534 {
 535   //Initialize the parameters of the analysis.
 536   fEMCALGeoName         = "EMCAL_COMPLETEV1";
 537   fPHOSGeoName          = "PHOSgeo";
 538   fEMCALGeoMatrixSet    = kFALSE;
 539   fPHOSGeoMatrixSet     = kFALSE;
 540   fRemoveBadChannels    = kFALSE;
 541   fNCellsFromPHOSBorder = 0;
 542
 543   //  fMaskCellColumns = new Int_t[fNMaskCellColumns];
 544   //  fMaskCellColumns[0] = 6 ;  fMaskCellColumns[1] = 7 ;  fMaskCellColumns[2] = 8 ;
 545   //  fMaskCellColumns[3] = 35;  fMaskCellColumns[4] = 36;  fMaskCellColumns[5] = 37;
 546   //  fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
 547   //  fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
 548   //  fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
 549
 550 }
 551
 552
 553 //_____________________________________________________
 554 void AliCalorimeterUtils::InitPHOSBadChannelStatusMap()
 555 {
 556   //Init PHOS bad channels map
 557   if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSBadChannelStatusMap()\n");
 558   //In order to avoid rewriting the same histograms
 559   Bool_t oldStatus = TH1::AddDirectoryStatus();
 560   TH1::AddDirectory(kFALSE);
 561
 562   fPHOSBadChannelMap = new TObjArray(5);
 563   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));
 564
 565   fPHOSBadChannelMap->SetOwner(kTRUE);
 566   fPHOSBadChannelMap->Compress();
 567
 568   //In order to avoid rewriting the same histograms
 569   TH1::AddDirectory(oldStatus);
 570 }
 571
 572 //______________________________________________________
 573 void AliCalorimeterUtils::InitPHOSRecalibrationFactors()
 574 {
 575         //Init EMCAL recalibration factors
 576         if(fDebug > 0 )printf("AliCalorimeterUtils::InitPHOSRecalibrationFactors()\n");
 577         //In order to avoid rewriting the same histograms
 578         Bool_t oldStatus = TH1::AddDirectoryStatus();
 579         TH1::AddDirectory(kFALSE);
 580
 581         fPHOSRecalibrationFactors = new TObjArray(5);
 582         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));
 583         //Init the histograms with 1
 584         for (Int_t m = 0; m < 5; m++) {
 585                 for (Int_t i = 0; i < 56; i++) {
 586                         for (Int_t j = 0; j < 64; j++) {
 587                                 SetPHOSChannelRecalibrationFactor(m,j,i,1.);
 588                         }
 589                 }
 590         }
 591         fPHOSRecalibrationFactors->SetOwner(kTRUE);
 592         fPHOSRecalibrationFactors->Compress();
 593
 594         //In order to avoid rewriting the same histograms
 595         TH1::AddDirectory(oldStatus);
 596 }
 597
 598
 599 //___________________________________________
 600 void AliCalorimeterUtils::InitEMCALGeometry()
 601 {
 602         //Initialize EMCAL geometry if it did not exist previously
 603         if (!fEMCALGeo){
 604                 fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
 605                 if(fDebug > 0){
 606                         printf("AliCalorimeterUtils::InitEMCALGeometry()");
 607                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 608                         printf("\n");
 609                 }
 610         }
 611 }
 612
 613 //__________________________________________
 614 void AliCalorimeterUtils::InitPHOSGeometry()
 615 {
 616         //Initialize PHOS geometry if it did not exist previously
 617         if (!fPHOSGeo){
 618                 fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
 619                 if(fDebug > 0){
 620                         printf("AliCalorimeterUtils::InitPHOSGeometry()");
 621                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 622                         printf("\n");
 623                 }
 624         }
 625 }
 626
 627 //_________________________________________________________
 628 void AliCalorimeterUtils::Print(const Option_t * opt) const
 629 {
 630
 631   //Print some relevant parameters set for the analysis
 632   if(! opt)
 633     return;
 634
 635   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 636   printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
 637   printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
 638          fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
 639   if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
 640   printf("Recalibrate Clusters? %d\n",fRecalibration);
 641   printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
 642   printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
 643   printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
 644
 645   printf("    \n") ;
 646 }
 647
 648 //__________________________________________________________________
 649 Bool_t AliCalorimeterUtils::MaskFrameCluster(const Int_t iSM,
 650                                              const Int_t ieta) const
 651 {
 652   //Check if cell is in one of the regions where we have significant amount
 653   //of material in front. Only EMCAL
 654
 655   Int_t icol = ieta;
 656   if(iSM%2) icol+=48; // Impair SM, shift index [0-47] to [48-96]
 657
 658   if (fNMaskCellColumns && fMaskCellColumns) {
 659     for (Int_t imask = 0; imask < fNMaskCellColumns; imask++) {
 660       if(icol==fMaskCellColumns[imask]) return kTRUE;
 661     }
 662   }
 663
 664   return kFALSE;
 665
 666 }
 667
 668
 669 //__________________________________________________________________________
 670 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
 671                                                       AliVCaloCells * cells)
 672 {
 673         // Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that compose the cluster.
 674
 675   //Initialize some used variables
 676         Float_t energy   = 0;
 677         Int_t absId      = -1;
 678         Int_t icol = -1, irow = -1, iRCU = -1, module=1;
 679         Float_t factor = 1, frac = 0;
 680
 681         if(cells) {
 682
 683     //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
 684     UShort_t * index    = cluster->GetCellsAbsId() ;
 685     Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
 686     Int_t ncells     = cluster->GetNCells();
 687     TString calo     = "EMCAL";
 688     if(cluster->IsPHOS()) calo = "PHOS";
 689
 690     //Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
 691     for(Int_t icell = 0; icell < ncells; icell++){
 692       absId = index[icell];
 693       frac =  fraction[icell];
 694       if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
 695       module = GetModuleNumberCellIndexes(absId,calo,icol,irow,iRCU);
 696       if(cluster->IsPHOS()) factor = GetPHOSChannelRecalibrationFactor (module,irow,icol);
 697       else                  factor = GetEMCALChannelRecalibrationFactor(module,icol,irow);
 698       if(fDebug>2)
 699         printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - recalibrate cell: %s, module %d, col %d, row %d, cell fraction %f, recalibration factor %f, cell energy %f\n",
 700                calo.Data(),module,icol,irow,frac,factor,cells->GetCellAmplitude(absId));
 701
 702       energy += cells->GetCellAmplitude(absId)*factor*frac;
 703     }
 704
 705     if(fDebug>1)
 706       printf("AliCalorimeterUtils::RecalibrateClusterEnergy() - Energy before %f, after %f\n",cluster->E(),energy);
 707
 708         }// cells available
 709   else{
 710     Fatal("RecalibrateClusterEnergy()","Cells pointer does not exist!");
 711   }
 712
 713         return energy;
 714 }
 715
 716 //________________________________________________________________________________
 717 void AliCalorimeterUtils::SetGeometryTransformationMatrices(AliVEvent* inputEvent)
 718 {
 719   //Set the calorimeters transformation matrices
 720
 721   //Get the EMCAL transformation geometry matrices from ESD
 722   if(!fEMCALGeoMatrixSet && fEMCALGeo){
 723     if(fLoadEMCALMatrices){
 724       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined EMCAL geometry matrices\n");
 725       for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 726         if(fEMCALMatrix[mod]){
 727           if(fDebug > 1)
 728             fEMCALMatrix[mod]->Print();
 729           fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
 730         }
 731       }//SM loop
 732       fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 733
 734     }//Load matrices
 735     else if (!gGeoManager) {
 736
 737       if(fDebug > 1)
 738         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load EMCAL misalignment matrices. \n");
 739       if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 740         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 741           //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
 742           if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod)) {
 743             fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
 744           }
 745         }// loop over super modules
 746
 747         fEMCALGeoMatrixSet = kTRUE;//At least one, so good
 748
 749       }//ESD as input
 750       else {
 751         if(fDebug > 1)
 752           printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 753       }//AOD as input
 754     }//Get matrix from data
 755     else if(gGeoManager){
 756       fEMCALGeoMatrixSet = kTRUE;
 757     }
 758   }//EMCAL geo && no geoManager
 759
 760         //Get the PHOS transformation geometry matrices from ESD
 761   if(!fPHOSGeoMatrixSet && fPHOSGeo){
 762
 763     if(fLoadPHOSMatrices){
 764       printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load user defined PHOS geometry matrices\n");
 765       for(Int_t mod = 0 ; mod < 5 ; mod++){
 766         if(fPHOSMatrix[mod]){
 767           if(fDebug > 1 )
 768             fPHOSMatrix[mod]->Print();
 769           fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
 770         }
 771       }//SM loop
 772       fPHOSGeoMatrixSet = kTRUE;//At least one, so good
 773
 774     }//Load matrices
 775     else if (!gGeoManager) {
 776       if(fDebug > 1)
 777         printf(" AliCalorimeterUtils::SetGeometryTransformationMatrices() - Load PHOS misalignment matrices. \n");
 778                         if(!strcmp(inputEvent->GetName(),"AliESDEvent"))  {
 779                                 for(Int_t mod = 0; mod < 5; mod++){
 780                                         if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod)) {
 781                                                 //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
 782                                                 fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
 783                                         }
 784                                 }// loop over modules
 785         fPHOSGeoMatrixSet  = kTRUE; //At least one so good
 786                         }//ESD as input
 787                         else {
 788                                 if(fDebug > 1)
 789                                         printf("AliCalorimeterUtils::SetGeometryTransformationMatrices() - Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file\n");
 790       }//AOD as input
 791     }// get matrix from data
 792     else if(gGeoManager){
 793       fPHOSGeoMatrixSet = kTRUE;
 794     }
 795         }//PHOS geo     and  geoManager was not set
 796
 797 }
 798
 799 //__________________________________________________________________________________________
 800 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
 801 {
 802
 803   //Recalculate EMCAL cluster position
 804
 805   fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
 806
 807 }
 808
 809 //________________________________________________________________________________
 810 void AliCalorimeterUtils::RecalculateClusterTrackMatching(AliVEvent * event,
 811                                                           TObjArray* clusterArray)
 812 {
 813   //Recalculate track matching
 814
 815   if (fRecalculateMatching) {
 816     fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo)   ;
 817     //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
 818     //if(esdevent){
 819     //  fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent)  ;
 820     //  fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
 821     //}
 822   }
 823 }