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