ANALYSIS/TenderSupplies/AliEMCALTenderSupply.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 ///////////////////////////////////////////////////////////////////////////////
  18 //                                                                           //
  19 //  EMCAL tender, apply corrections to EMCAl clusters                        //
  20 //  and do track matching                                                    //
  21 //  Author : Deepa Thomas (Utrecht University)                               //
  22 //                                                                           //
  23 ///////////////////////////////////////////////////////////////////////////////
  24
  25 #include "TROOT.h"
  26 #include "TFile.h"
  27 #include "TTree.h"
  28 #include "TInterpreter.h"
  29 #include "TObjArray.h"
  30 #include "TClonesArray.h"
  31
  32 #include <AliLog.h>
  33 #include <AliESDEvent.h>
  34 #include <AliAnalysisManager.h>
  35 #include <AliTender.h>
  36 #include "AliOADBContainer.h"
  37 #include "AliCDBManager.h"
  38 #include "AliCDBStorage.h"
  39 #include "AliCDBEntry.h"
  40 #include "AliMagF.h"
  41 #include "TGeoGlobalMagField.h"
  42
  43 #include "AliESDCaloCluster.h"
  44 #include "AliEMCALTenderSupply.h"
  45
  46 //EMCAL
  47 #include "AliEMCALGeometry.h"
  48 #include "AliEMCALRecoUtils.h"
  49 #include "AliEMCALClusterizer.h"
  50 #include "AliEMCALRecParam.h"
  51 #include "AliEMCALRecPoint.h"
  52 #include "AliEMCALAfterBurnerUF.h"
  53 #include "AliEMCALClusterizerNxN.h"
  54 #include "AliEMCALClusterizerv1.h"
  55 #include "AliEMCALClusterizerv2.h"
  56 #include "AliEMCALDigit.h"
  57
  58 ClassImp(AliEMCALTenderSupply)
  59
  60 // mfasel:
  61 //  Remove all calls to TGrid::Connect - grid connection is global and better steer by the run macro
  62
  63 AliEMCALTenderSupply::AliEMCALTenderSupply() :
  64         AliTenderSupply()
  65         ,fEMCALGeo(0x0)
  66         ,fEMCALGeoName("EMCAL_FIRSTYEARV1")
  67 //      ,fEMCALRecoUtils(new AliEMCALRecoUtils)
  68         ,fEMCALRecoUtils(0)
  69         ,fConfigName("")
  70         ,fDebugLevel(0)
  71         ,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)
  72         ,fNonLinearThreshold(30)
  73         ,fReCalibCluster(kFALSE)
  74         ,fReCalibCell(kFALSE)
  75         ,fRecalClusPos(kFALSE)
  76         ,fFiducial(kFALSE)
  77         ,fNCellsFromEMCALBorder(1)
  78         ,fRecalDistToBadChannels(kFALSE)
  79         ,fInputTree(0)
  80         ,fInputFile(0)
  81         ,fFilepass(0)
  82         ,fMass(0.139)
  83         ,fStep(1)
  84         ,fRcut(0.05)
  85         ,fBasePath(".")
  86         ,fReClusterize(kFALSE)
  87         ,fClusterizer(0)
  88         ,fGeomMatrixSet(kFALSE)
  89         ,fLoadGeomMatrices(kFALSE)
  90         ,fRecParam(0)
  91         ,fOCDBpath(" ")
  92         ,fUnfolder(0)
  93         ,fDigitsArr(0)
  94         ,fClusterArr(0)
  95
  96 {
  97         //
  98         // default ctor
  99         //
 100         for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
 101         for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
 102         fRecParam        = new AliEMCALRecParam;
 103         fEMCALRecoUtils = new AliEMCALRecoUtils();
 104         fDigitsArr       = new TClonesArray("AliEMCALDigit",1000);
 105
 106 }
 107
 108 //_____________________________________________________
 109 AliEMCALTenderSupply::AliEMCALTenderSupply(const char *name, const AliTender *tender) :
 110         AliTenderSupply(name,tender)
 111         ,fEMCALGeo(0x0)
 112         ,fEMCALGeoName("EMCAL_FIRSTYEARV1")
 113         ,fEMCALRecoUtils(0)
 114         ,fConfigName("")
 115         ,fDebugLevel(0)
 116         ,fNonLinearFunc(AliEMCALRecoUtils::kNoCorrection)
 117         ,fNonLinearThreshold(30)
 118         ,fReCalibCluster(kFALSE)
 119         ,fReCalibCell(kFALSE)
 120         ,fRecalClusPos(kFALSE)
 121         ,fFiducial(kFALSE)
 122         ,fNCellsFromEMCALBorder(1)
 123         ,fRecalDistToBadChannels(kFALSE)
 124         ,fInputTree(0)
 125         ,fInputFile(0)
 126         ,fFilepass(0)
 127         ,fMass(0.139)
 128         ,fStep(1)
 129         ,fRcut(0.05)
 130         ,fBasePath(".")
 131         ,fReClusterize(kFALSE)
 132         ,fClusterizer(0)
 133         ,fGeomMatrixSet(kFALSE)
 134         ,fLoadGeomMatrices(kFALSE)
 135         ,fRecParam(0)
 136         ,fOCDBpath(" ")
 137         ,fUnfolder(0)
 138         ,fDigitsArr(0)
 139         ,fClusterArr(0)
 140
 141 {
 142         //
 143         // named ctor
 144         //
 145         for(Int_t i = 0; i < 10; i++) fEMCALMatrix[i] = 0 ;
 146         for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
 147         fRecParam        = new AliEMCALRecParam;
 148         fEMCALRecoUtils = new AliEMCALRecoUtils();
 149         fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
 150 }
 151
 152 //_____________________________________________________
 153 AliEMCALTenderSupply::~AliEMCALTenderSupply()
 154 {
 155         //Destructor
 156         delete fEMCALGeo;
 157         delete fEMCALRecoUtils;
 158         delete fInputTree;
 159         delete fInputFile;
 160         delete fClusterizer;
 161         delete fUnfolder;
 162
 163         if (fDigitsArr){
 164                 fDigitsArr->Clear("C");
 165                 delete fDigitsArr;
 166         }
 167
 168 }
 169
 170 //_____________________________________________________
 171 void AliEMCALTenderSupply::Init()
 172 {
 173         //
 174         // Initialise EMCAL tender
 175         //
 176
 177         if (fDebugLevel>0) AliInfo("Init EMCAL Tender supply\n");
 178
 179         AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
 180
 181         fInputTree = mgr->GetTree();
 182
 183         if(gROOT->LoadMacro(fConfigName) >=0){
 184                 AliDebug(1, Form("Loading settings from macro %s", fConfigName.Data()));
 185                 AliEMCALTenderSupply *tender = (AliEMCALTenderSupply*)gInterpreter->ProcessLine("ConfigEMCALTenderSupply()");
 186                 fDebugLevel         = tender->fDebugLevel;
 187                 fEMCALGeoName       = tender->fEMCALGeoName;
 188                 fEMCALRecoUtils     = tender->fEMCALRecoUtils;
 189                 fConfigName         = tender->fConfigName;
 190                 fNonLinearFunc      = tender->fNonLinearFunc;
 191                 fNonLinearThreshold = tender->fNonLinearThreshold;
 192                 fReCalibCluster     = tender->fReCalibCluster;
 193                 fReCalibCell        = tender->fReCalibCell;
 194                 fRecalClusPos       = tender->fRecalClusPos;
 195                 fFiducial           = tender->fFiducial;
 196                 fNCellsFromEMCALBorder = tender->fNCellsFromEMCALBorder;
 197                 fRecalDistToBadChannels = tender->fRecalDistToBadChannels;
 198                 fMass               = tender->fMass;
 199                 fStep               = tender->fStep;
 200                 fRcut               = tender->fRcut;
 201                 fReClusterize       = tender->fReClusterize;
 202                 fLoadGeomMatrices   = tender->fLoadGeomMatrices;
 203                 fRecParam           = tender->fRecParam;
 204                 fOCDBpath           = tender->fOCDBpath;
 205                 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = tender->fGeomMatrix[i] ;
 206         }
 207
 208         // Init goemetry
 209         fEMCALGeo =  AliEMCALGeometry::GetInstance(fEMCALGeoName) ;
 210
 211
 212         //Initialising Non linearity parameters
 213         fEMCALRecoUtils->SetNonLinearityThreshold(fNonLinearThreshold);
 214         fEMCALRecoUtils->SetNonLinearityFunction(fNonLinearFunc);
 215
 216         //Setting mass, step size and residual cut
 217         fEMCALRecoUtils->SwitchOnCutEtaPhiSum();
 218         fEMCALRecoUtils->SetCutR(fRcut);
 219         fEMCALRecoUtils->SetMass(fMass);
 220         fEMCALRecoUtils->SetStep(fStep);
 221
 222         //AliLog::SetGlobalDebugLevel(1);
 223
 224         if(fDebugLevel>1) fEMCALRecoUtils->Print("");
 225
 226
 227 }
 228
 229 //_____________________________________________________
 230 void AliEMCALTenderSupply::ProcessEvent()
 231 {
 232         //Event loop
 233         AliESDEvent *event=fTender->GetEvent();
 234         if (!event) return;
 235
 236         if(fTender->RunChanged()){
 237                 //Initialising parameters once per run number
 238                 if (!InitBadChannels()) return;
 239                 if (fRecalClusPos){ if (!InitMisalignMatrix()) return;}
 240                 if (fReCalibCluster || fReCalibCell){ if (!InitRecalib()) return;}
 241                 if (fReClusterize) InitClusterization();
 242         }
 243
 244         AliESDCaloCells *cells= event->GetEMCALCells();
 245
 246         //------------ EMCAL cells loop ------------
 247
 248         //Recalibrate cells
 249         if(fReCalibCell) RecalibrateCells();
 250
 251         //------------ Reclusterize ------------
 252         if(fReClusterize){
 253                 FillDigitsArray();
 254                 Clusterize();
 255                 UpdateClusters();
 256         }
 257
 258         //------------Good clusters------------
 259         TClonesArray *clusArr;
 260
 261         clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
 262         if(!clusArr) clusArr = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
 263         if(!clusArr) return;
 264
 265         Int_t nclusters = clusArr->GetEntriesFast(); //bug fix by Rongrong
 266         for (Int_t icluster=0; icluster < nclusters; icluster++ ){
 267                 AliVCluster *clust = static_cast<AliVCluster*>(clusArr->At(icluster));
 268                 if(!clust) continue;
 269                 if (!clust->IsEMCAL()) continue;
 270                 if(fEMCALRecoUtils->ClusterContainsBadChannel(fEMCALGeo, clust->GetCellsAbsId(), clust->GetNCells() )){
 271                         delete clusArr->RemoveAt(icluster);
 272                         continue;
 273                 }
 274                 if(fFiducial){
 275                         if(!fEMCALRecoUtils->CheckCellFiducialRegion(fEMCALGeo, clust, cells)){
 276                                 delete clusArr->RemoveAt(icluster);
 277                                 continue;
 278                         }
 279                 }
 280                 fEMCALRecoUtils->CorrectClusterEnergyLinearity(clust);
 281                 if(fRecalDistToBadChannels) fEMCALRecoUtils->RecalculateClusterDistanceToBadChannel(fEMCALGeo, cells, clust);
 282                 if(fReCalibCluster) fEMCALRecoUtils->RecalibrateClusterEnergy(fEMCALGeo, clust, cells);
 283                 if(fRecalClusPos) fEMCALRecoUtils->RecalculateClusterPosition(fEMCALGeo, cells, clust);
 284                 //fEMCALRecoUtils->SetTimeDependentCorrections(event->GetRunNumber());
 285         }
 286         clusArr->Compress();
 287
 288
 289         //-------- Track Matching ------------------
 290
 291         //set magnetic field
 292         Double_t  magF = event->GetMagneticField();
 293         Double_t magSign = 1.0;
 294         if(magF<0)magSign = -1.0;
 295
 296         if (!TGeoGlobalMagField::Instance()->GetField())
 297         {
 298                 AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG);
 299                 TGeoGlobalMagField::Instance()->SetField(field);
 300         }
 301
 302
 303         fEMCALRecoUtils->FindMatches(event,0x0,fEMCALGeo);
 304
 305         Int_t nTracks = event->GetNumberOfTracks();
 306
 307         if(nTracks>0){
 308                 SetClusterMatchedToTrack(event);
 309                 SetTracksMatchedToCluster(event);
 310         }
 311
 312 }
 313
 314 //_____________________________________________________
 315 void AliEMCALTenderSupply::SetClusterMatchedToTrack(AliESDEvent *event)
 316 {
 317         //checks if tracks are matched to EMC clusters and set the matched EMCAL cluster index to ESD track.
 318
 319         Int_t matchClusIndex = -1;
 320         Int_t nTracks = event->GetNumberOfTracks();
 321
 322         // Track loop
 323         for (Int_t iTrack = 0; iTrack < nTracks; iTrack++)
 324         {
 325                 AliESDtrack* track = event->GetTrack(iTrack);
 326                 if (!track) {
 327                         printf("ERROR: Could not receive track %d\n", iTrack);
 328                         continue;
 329                 }
 330
 331                 matchClusIndex = fEMCALRecoUtils->GetMatchedClusterIndex(iTrack);
 332                 track->SetEMCALcluster(matchClusIndex); //sets -1 if track not matched within residual
 333                 if(matchClusIndex != -1) track->SetStatus(AliESDtrack::kEMCALmatch);
 334         }
 335
 336         if (fDebugLevel>2) AliInfo("Track matched to closest cluster\n");
 337 }
 338
 339 //_____________________________________________________
 340 void AliEMCALTenderSupply::SetTracksMatchedToCluster(AliESDEvent *event)
 341 {
 342         //checks if EMC clusters are matched to ESD track.
 343         //Adds track indexes of all the tracks matched to a cluster withing residuals in ESDCalocluster
 344
 345         Int_t matchTrackIndex = -1;
 346         Int_t nTracks = event->GetNumberOfTracks();
 347
 348         // Cluster loop
 349         for (Int_t iClus=0; iClus < event->GetNumberOfCaloClusters(); iClus++)
 350         {
 351                 AliESDCaloCluster * cluster = event->GetCaloCluster(iClus);
 352                 if (!cluster->IsEMCAL()) continue;
 353
 354                 Int_t nMatched =0;
 355                 TArrayI arrayTrackMatched(nTracks);
 356
 357                 //get the closest track matched to the cluster
 358                 matchTrackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iClus);
 359                 arrayTrackMatched[nMatched] = matchTrackIndex;
 360                 nMatched++;
 361
 362                 //get all other tracks matched to the cluster
 363                 //track loop
 364                 for(Int_t iTrk=0; iTrk<nTracks; iTrk++){
 365                         AliESDtrack* track = event->GetTrack(iTrk);
 366                         if(iTrk == matchTrackIndex) continue;
 367
 368                         if(track->GetEMCALcluster() == iClus){
 369                                 arrayTrackMatched[nMatched] = iTrk;
 370                                 nMatched++;
 371                         }
 372                 }
 373
 374                 arrayTrackMatched.Set(nMatched);
 375                 cluster->AddTracksMatched(arrayTrackMatched);
 376
 377                 Float_t eta= -999, phi = -999;
 378                 if(matchTrackIndex != -1) fEMCALRecoUtils->GetMatchedResiduals(iClus, eta, phi);
 379                 cluster->SetTrackDistance(phi, eta);
 380         }
 381         if (fDebugLevel>2) AliInfo("Cluster matched to tracks \n");
 382 }
 383
 384 //_____________________________________________________
 385 Bool_t AliEMCALTenderSupply::InitMisalignMatrix()
 386 {
 387         //Initialising Misalignment matrix
 388
 389         AliESDEvent *event=fTender->GetEvent();
 390         if (!event) return kFALSE;
 391
 392         if (fDebugLevel>0) AliInfo("Initialising Misalignment matrix \n");
 393
 394         if(fInputTree){
 395                 fInputFile = fInputTree->GetCurrentFile();
 396                 if(fInputFile){
 397                         const char *fileName = fInputFile->GetName();
 398                         TString FileName = TString(fileName);
 399                         if     (FileName.Contains("pass1")) fFilepass = TString("pass1");
 400                         else if(FileName.Contains("pass2")) fFilepass = TString("pass2");
 401                         else if(FileName.Contains("pass3")) fFilepass = TString("pass3");
 402                         else AliError("pass number not found");
 403                 }
 404                 else AliError("File not found");
 405         }
 406         else AliError("Tree not found");
 407
 408         Int_t runGM = 0;
 409         runGM = event->GetRunNumber();
 410
 411         fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
 412
 413         if(runGM <=140000){ //2010 data
 414
 415                 Double_t rotationMatrix[4][9] = {{-0.014587, -0.999892, -0.002031, 0.999892, -0.014591,  0.001979, -0.002009, -0.002002,  0.999996},
 416                         {-0.014587,  0.999892,  0.002031, 0.999892,  0.014591, -0.001979, -0.002009,  0.002002, -0.999996},
 417                         {-0.345864, -0.938278, -0.003412, 0.938276, -0.345874,  0.003010, -0.004004, -0.002161,  0.999990},
 418                         {-0.345861,  0.938280,  0.003412, 0.938276,  0.345874, -0.003010, -0.004004,  0.002161, -0.999990}};
 419
 420                 Double_t translationMatrix[4][3] = {{0.351659,    447.576446,  176.269742},
 421                         {1.062577,    446.893974, -173.728870},
 422                         {-154.213287, 419.306156,  176.753692},
 423                         {-153.018950, 418.623681, -173.243605}};
 424
 425                 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
 426                 {
 427                         //if(DebugLevel() > 1)  fEMCALMatrix[mod]->Print();
 428                         fEMCALMatrix[mod] = new TGeoHMatrix();           // mfasel: prevent tender from crashing
 429                         fEMCALMatrix[mod]->SetRotation(rotationMatrix[mod]);
 430                         fEMCALMatrix[mod]->SetTranslation(translationMatrix[mod]);
 431                         fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
 432                 }
 433         }
 434
 435
 436         else if(runGM>140000 && runGM <148531 && (fFilepass = "pass1"))
 437         { // 2011 LHC11a pass1 data
 438                 AliOADBContainer emcalgeoCont(Form("emcal2011"));
 439                 emcalgeoCont.InitFromFile(Form("%s/BetaGood.root",fBasePath.Data()),Form("AliEMCALgeo"));
 440
 441                 TObjArray *mobj=(TObjArray*)emcalgeoCont.GetObject(100,"survey11byS");
 442
 443                 for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 444                         fEMCALMatrix[mod] = (TGeoHMatrix*) mobj->At(mod);
 445                         fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod);
 446
 447                         fEMCALMatrix[mod]->Print();
 448                 }
 449         }
 450
 451         else AliInfo("MISALLIGNMENT NOT APPLIED");
 452
 453         return kTRUE;
 454 }
 455
 456 //_____________________________________________________
 457 Bool_t AliEMCALTenderSupply::InitBadChannels()
 458 {
 459         //Initialising Bad channel maps
 460
 461         AliESDEvent *event=fTender->GetEvent();
 462         if (!event) return kFALSE;
 463
 464         Int_t fRunBC= 0;
 465         fRunBC = event->GetRunNumber();
 466
 467         if (fDebugLevel>0) AliInfo("Initialising Bad channel map \n");
 468
 469         if(fInputTree){
 470                 fInputFile = fInputTree->GetCurrentFile();
 471                 if(fInputFile){
 472                         const char *fileName = fInputFile->GetName();
 473                         TString FileName = TString(fileName);
 474                         if     (FileName.Contains("pass1")) fFilepass = TString("pass1");
 475                         else if(FileName.Contains("pass2")) fFilepass = TString("pass2");
 476                         else if(FileName.Contains("pass3")) fFilepass = TString("pass3");
 477                         else AliError("pass number not found");
 478                 }
 479                 else AliError("File not found");
 480         }
 481         else AliError("Tree not found");
 482
 483         if(fFiducial){
 484                 fEMCALRecoUtils->SetNumberOfCellsFromEMCALBorder(fNCellsFromEMCALBorder);
 485                 fEMCALRecoUtils->SwitchOnNoFiducialBorderInEMCALEta0();
 486         }
 487
 488         fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
 489         if(fRecalDistToBadChannels) fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
 490
 491         TFile *fbad;
 492         //2010
 493         if(fRunBC <=140000){
 494                 fbad = new TFile(Form("%s/BadChannels.root",fBasePath.Data()),"read");
 495                 if (fbad->IsZombie()){
 496                         TString fPath = TString(fBasePath.Data());
 497                         if(fPath.Contains("alien")) {
 498                                 if (fDebugLevel>1) AliInfo("Connecting to alien to get BadChannels.root \n");
 499                                 fbad = TFile::Open(Form("%s/BadChannelsDB/BadChannels.root",fBasePath.Data()));
 500                         }
 501                 }
 502
 503                 TH2I * hb0 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod0");
 504                 TH2I * hb1 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod1");
 505                 TH2I * hb2 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod2");
 506                 TH2I * hb3 = ( TH2I *)fbad->Get("EMCALBadChannelMap_Mod3");
 507                 fEMCALRecoUtils->SetEMCALChannelStatusMap(0,hb0);
 508                 fEMCALRecoUtils->SetEMCALChannelStatusMap(1,hb1);
 509                 fEMCALRecoUtils->SetEMCALChannelStatusMap(2,hb2);
 510                 fEMCALRecoUtils->SetEMCALChannelStatusMap(3,hb3);
 511
 512         }
 513
 514         //2011
 515         Int_t nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1;
 516
 517         if(fRunBC>=144871 && fRunBC<=146860){ //LHC11a 2.76 TeV pp
 518
 519                 if(fFilepass = "pass1"){ // pass1
 520
 521                         const Int_t nTowers=89;
 522                         Int_t hotChannels[nTowers]={74, 103, 152, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 368, 369, 370, 371, 372, 373, 374,375, 376, 377, 378, 379, 380, 381, 382, 383, 917, 1275, 1288, 1519, 1595, 1860, 1967, 2022, 2026, 2047, 2117, 2298, 2540, 2776, 3135, 3764, 6095, 6111, 6481, 6592, 6800, 6801, 6802, 6803, 6804, 6805, 6806, 6807, 6808, 6809, 6810, 6811, 6812, 6813, 6814, 6815, 7371, 7425, 7430, 7457, 7491, 7709, 8352, 8353, 8356, 8357, 8808, 8810, 8812, 8814, 9056, 9769, 9815, 9837};
 523                         for(Int_t i=0; i<nTowers; i++)
 524                         {
 525                                 fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
 526
 527                                 fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
 528                                 fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
 529                         }
 530                 }
 531
 532                 if(fFilepass = "pass2"){ // pass2
 533                         const Int_t nTowers=24;
 534                         Int_t hotChannels[nTowers]= {74, 103, 152, 917, 1059, 1175, 1276, 1288, 1376, 1382, 1595, 2022, 2026, 2210, 2540, 2778, 2793, 3135, 3764, 5767, 6481, 7371, 7878, 9769};
 535                         for(Int_t i=0; i<nTowers; i++)
 536                         {
 537                                 fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
 538
 539                                 fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
 540                                 fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
 541                         }
 542                 }
 543
 544         }
 545
 546         if(fRunBC>=151636 && fRunBC<=155384) //LHC11c : 7TeV by Rongrong
 547         {
 548                 const Int_t nTowers=8;
 549                 Int_t hotChannels[nTowers]={917, 2115, 2123, 2540, 6481, 9815, 10113, 10115};
 550
 551                 for(Int_t i=0; i<nTowers; i++)
 552                 {
 553                         fEMCALGeo->GetCellIndex(hotChannels[i],nSupMod,nModule,nIphi,nIeta);
 554                         fEMCALGeo->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta,iphi,ieta);
 555                         fEMCALRecoUtils->SetEMCALChannelStatus(nSupMod, ieta, iphi);
 556                 }
 557         }
 558
 559
 560         return kTRUE;
 561 }
 562
 563 //_____________________________________________________
 564 Bool_t AliEMCALTenderSupply::InitRecalib()
 565 {
 566         //Initialising Recalibration Factors
 567
 568         AliESDEvent *event=fTender->GetEvent();
 569         if (!event) return kFALSE;
 570
 571         if (fDebugLevel>0) AliInfo("Initialising Recalibration factors \n");
 572
 573         if(fInputTree){
 574                 fInputFile = fInputTree->GetCurrentFile();
 575                 if(fInputFile){
 576                         const char *fileName = fInputFile->GetName();
 577                         TString FileName = TString(fileName);
 578                         if     (FileName.Contains("pass1")) fFilepass = TString("pass1");
 579                         else if(FileName.Contains("pass2")) fFilepass = TString("pass2");
 580                         else if(FileName.Contains("pass3")) fFilepass = TString("pass3");
 581                         else AliError("pass number not found");
 582                 }
 583                 else AliError("File not found");
 584         }
 585         else AliError("Tree not found");
 586         //        else {cout << "Tree not found " <<endl; return kFALSE;}
 587
 588         Int_t runRC = event->GetRunNumber();
 589
 590         //if (event->GetRunNumber()==runRC)
 591         //     return kFALSE;
 592
 593         fEMCALRecoUtils->SwitchOnRecalibration();
 594
 595         TFile* fRecalib;
 596
 597         if(runRC <=140000){
 598                 fRecalib = new TFile(Form("%s/RecalibrationFactors.root",fBasePath.Data()),"read");
 599                 if (fRecalib->IsZombie()){
 600                         TString fPath = TString(fBasePath.Data());
 601                         if(fPath.Contains("alien")) {
 602
 603                                 if (fDebugLevel>1) AliInfo("Connecting to alien to get RecalibrationFactors.root \n");
 604
 605                                 if(     (runRC >= 114737 && runRC <= 117223 && (fFilepass = "pass2")) || //LHC10b pass2
 606                                                 (runRC >= 118503 && runRC <= 121040 && ((fFilepass = "pass2")||(fFilepass = "pass3"))) || //LHC10c pass2, LHC10c pass3
 607                                                 (runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass1")) || //LHC10d pass1
 608                                                 (runRC >= 127712 && runRC <= 130850 && (fFilepass = "pass1")) || //LHC10e pass1
 609                                                 (runRC >= 133004 && runRC < 134657  && (fFilepass = "pass1")))// LHC10f pass1 <134657
 610                                 {
 611                                         fRecalib = TFile::Open(Form("%s/RecalDB/summer_december_2010/RecalibrationFactors.root",fBasePath.Data()));
 612                                 }
 613
 614                                 else if((runRC >= 122195 && runRC <= 126437 && (fFilepass = "pass2")) || //LHC10d pass2
 615                                                 (runRC >= 134657 && runRC <= 135029 && (fFilepass = "pass1")) || //LHC10f pass1 >= 134657
 616                                                 (runRC >= 135654 && runRC <= 136377 && (fFilepass = "pass1")) || //LHC10g pass1
 617                                                 (runRC >= 136851 && runRC < 137231  && (fFilepass = "pass1"))) //LHC10h pass1 untill christmas
 618                                 {
 619                                         fRecalib = TFile::Open(Form("%s/RecalDB/december2010/RecalibrationFactors.root",fBasePath.Data()));
 620                                 }
 621
 622                                 else if(runRC >= 137231 && runRC <= 139517 && (fFilepass = "pass1")) //LHC10h pass1 from christmas
 623                                 {
 624                                         fRecalib = TFile::Open(Form("%s/RecalDB/summer2010/RecalibrationFactors.root",fBasePath.Data()));
 625                                 }
 626                                 else {
 627                                         AliError("Run number or pass number not found; RECALIBRATION NOT APPLIED");
 628                                 }
 629
 630                         }
 631                 }
 632
 633                 TH2F * r0 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM0");
 634                 TH2F * r1 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM1");
 635                 TH2F * r2 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM2");
 636                 TH2F * r3 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM3");
 637                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(0,r0);
 638                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(1,r1);
 639                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(2,r2);
 640                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(3,r3);
 641         }
 642
 643         if(runRC > 140000){
 644                 fRecalib = new TFile(Form("%s/RecalibrationFactors.root",fBasePath.Data()),"read");
 645                 if (fRecalib->IsZombie()){
 646                         TString fPath = TString(fBasePath.Data());
 647                         if(fPath.Contains("alien")) {
 648                                 if (fDebugLevel>1) AliInfo("Connecting to alien to get RecalibrationFactors.root \n");
 649
 650                                 fRecalib = TFile::Open(Form("%s/RecalDB/2011_v0/RecalibrationFactors.root",fBasePath.Data()));
 651                                 if(!fRecalib) AliError("Recalibration file not found");
 652                         }
 653                 }
 654
 655                 TH2F * r0 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM0");
 656                 TH2F * r1 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM1");
 657                 TH2F * r2 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM2");
 658                 TH2F * r3 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM3");
 659                 TH2F * r4 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM4");
 660                 TH2F * r5 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM5");
 661                 TH2F * r6 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM6");
 662                 TH2F * r7 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM7");
 663                 TH2F * r8 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM8");
 664                 TH2F * r9 = ( TH2F *)fRecalib->Get("EMCALRecalFactors_SM9");
 665
 666                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(0,r0);
 667                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(1,r1);
 668                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(2,r2);
 669                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(3,r3);
 670                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(4,r4);
 671                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(5,r5);
 672                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(6,r6);
 673                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(7,r7);
 674                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(8,r8);
 675                 fEMCALRecoUtils->SetEMCALChannelRecalibrationFactors(9,r9);
 676
 677         }
 678
 679         return kTRUE;
 680 }
 681
 682 //_____________________________________________________
 683 void AliEMCALTenderSupply::RecalibrateCells()
 684 {
 685         AliESDCaloCells *cells = fTender->GetEvent()->GetEMCALCells();
 686
 687         Int_t nEMCcell = cells->GetNumberOfCells();
 688         Double_t calibFactor = 1.;
 689
 690         for(Int_t icell=0; icell<nEMCcell; icell++){
 691                 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 692                 fEMCALGeo->GetCellIndex(cells->GetCellNumber(icell),imod,iTower,iIphi,iIeta);
 693                 fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 694                 calibFactor = fEMCALRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 695                 cells->SetCell(icell,cells->GetCellNumber(icell),cells->GetAmplitude(icell)*calibFactor,cells->GetTime(icell));
 696         }
 697
 698 }
 699
 700 //_____________________________________________________
 701 void AliEMCALTenderSupply::InitClusterization()
 702 {
 703         //Initialing clusterization/unfolding algorithm and set all the needed parameters
 704         AliESDEvent *event=fTender->GetEvent();
 705         if (!event){
 706                 AliWarning("Event not available!!!");
 707                 return;
 708         }
 709
 710         if (fDebugLevel>0) AliInfo(Form("Initialising Reclustering parameters: Clusterizer-%d \n",fRecParam->GetClusterizerFlag()));
 711
 712         //---set the geometry matrix
 713
 714         if (!fGeomMatrixSet) {
 715                 if (fLoadGeomMatrices) {
 716                         for(Int_t mod=0; mod < fEMCALGeo->GetNumberOfSuperModules(); ++mod) {
 717                                 if(fGeomMatrix[mod]){
 718                                         if(fDebugLevel > 2) fGeomMatrix[mod]->Print();
 719                                         fEMCALGeo->SetMisalMatrix(fGeomMatrix[mod],mod);
 720                                 }
 721                         }
 722                 } else { // get matrix from file (work around bug in aliroot)
 723                         for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 724                                 if(fDebugLevel > 2) event->GetEMCALMatrix(mod)->Print();
 725                                 if(event->GetEMCALMatrix(mod)) fEMCALGeo->SetMisalMatrix(event->GetEMCALMatrix(mod),mod) ;
 726                         }
 727                 }
 728                 fGeomMatrixSet=kTRUE;
 729         }
 730
 731         //---setup clusterizer
 732
 733         delete fClusterizer;
 734         if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 735                 fClusterizer = new AliEMCALClusterizerv1 (fEMCALGeo);
 736         else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 737                 fClusterizer = new AliEMCALClusterizerv2(fEMCALGeo);
 738         else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 739                 fClusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
 740         else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerv2) { //FIX this other way.
 741                 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fEMCALGeo);
 742                 clusterizer->SetNRowDiff(2);
 743                 clusterizer->SetNColDiff(2);
 744                 fClusterizer = clusterizer;
 745         } else{
 746                 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 747         }
 748
 749         //--- set the parameters
 750
 751         fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 752         fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 753         fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 754         fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 755         fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 756         fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 757         fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 758         fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 759         fClusterizer->SetInputCalibrated       ( kTRUE                               );
 760         fClusterizer->SetJustClusters          ( kTRUE                               );
 761
 762         //In case of unfolding after clusterization is requested, set the corresponding parameters
 763         if(fRecParam->GetUnfold()){
 764                 Int_t i=0;
 765                 for (i = 0; i < 8; i++) {
 766                         fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 767                 }//end of loop over parameters
 768                 for (i = 0; i < 3; i++) {
 769                         fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 770                         fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 771                 }//end of loop over parameters
 772
 773                 fClusterizer->InitClusterUnfolding();
 774         }// to unfold
 775
 776         fClusterizer->SetDigitsArr(fDigitsArr);
 777         fClusterizer->SetOutput(0);
 778         fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 779 }
 780 //_____________________________________________________
 781 void AliEMCALTenderSupply::FillDigitsArray()
 782 {
 783         //Fill digits from cells to a TClonesArray
 784         fDigitsArr->Clear("C");
 785
 786         AliESDEvent *event=fTender->GetEvent();
 787         if (!event){
 788                 AliWarning("Event not available!!!");
 789                 return;
 790         }
 791         AliESDCaloCells *cells = event->GetEMCALCells();
 792         Int_t ncells = cells->GetNumberOfCells();
 793         for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 794                 Double_t cellAmplitude=0, cellTime=0;
 795                 Short_t cellNumber=0;
 796                 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 797                         break;
 798                 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 799                 digit->SetId(cellNumber);
 800                 digit->SetTime(cellTime);
 801                 digit->SetTimeR(cellTime);
 802                 digit->SetIndexInList(idigit);
 803                 digit->SetType(AliEMCALDigit::kHG);
 804                 digit->SetAmplitude(cellAmplitude);
 805                 idigit++;
 806         }
 807 }
 808
 809 //_____________________________________________________
 810 void AliEMCALTenderSupply::Clusterize()
 811 {
 812         //Clusterize
 813         fClusterizer->Digits2Clusters("");
 814 }
 815
 816 //_____________________________________________________
 817 void AliEMCALTenderSupply::UpdateClusters()
 818 {
 819         //Update ESD cluster list
 820         AliESDEvent *event=fTender->GetEvent();
 821         if (!event){
 822                 AliWarning("Event not available!!!");
 823                 return;
 824         }
 825
 826         TClonesArray *clus;
 827         clus = dynamic_cast<TClonesArray*>(event->FindListObject("caloClusters"));
 828         if(!clus) clus = dynamic_cast<TClonesArray*>(event->FindListObject("CaloClusters"));
 829         if(!clus) return;
 830
 831         Int_t nents = clus->GetEntriesFast();
 832         for (Int_t i=0;i<nents;++i) {
 833                 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
 834                 if (!c)
 835                         continue;
 836                 if (c->IsEMCAL())
 837                         clus->RemoveAt(i);
 838         }
 839         clus->Compress();
 840         RecPoints2Clusters(clus);
 841 }
 842
 843 //_____________________________________________________
 844 void AliEMCALTenderSupply::RecPoints2Clusters(TClonesArray *clus)
 845 {
 846         // Convert AliEMCALRecoPoints to AliESDCaloClusters
 847         // Cluster energy, global position, cells and their amplitude fractions are restored.
 848         Bool_t esdobjects = 0;
 849         if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
 850                 esdobjects = 1;
 851
 852         AliESDEvent *event=fTender->GetEvent();
 853         if (!event){
 854                 AliWarning("Event not available!!!");
 855                 return;
 856         }
 857
 858         Int_t ncls = fClusterArr->GetEntriesFast();
 859         for(Int_t i=0, nout=clus->GetEntriesFast(); i < ncls; ++i) {
 860                 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 861
 862                 Int_t ncells_true = 0;
 863                 const Int_t ncells = recpoint->GetMultiplicity();
 864                 UShort_t   absIds[ncells];
 865                 Double32_t ratios[ncells];
 866                 Int_t *dlist = recpoint->GetDigitsList();
 867                 Float_t *elist = recpoint->GetEnergiesList();
 868                 for (Int_t c = 0; c < ncells; ++c) {
 869                         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 870                         absIds[ncells_true] = digit->GetId();
 871                         ratios[ncells_true] = elist[c]/digit->GetAmplitude();
 872                         if (ratios[ncells_true] < 0.001)
 873                                 continue;
 874                         ++ncells_true;
 875                 }
 876
 877                 if (ncells_true < 1) {
 878                         AliWarning("Skipping cluster with no cells");
 879                         continue;
 880                 }
 881
 882                 // calculate new cluster position
 883                 TVector3 gpos;
 884                 recpoint->GetGlobalPosition(gpos);
 885                 Float_t g[3];
 886                 gpos.GetXYZ(g);
 887
 888                 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->New(nout++));
 889                 c->SetType(AliVCluster::kEMCALClusterv1);
 890                 c->SetE(recpoint->GetEnergy());
 891                 c->SetPosition(g);
 892                 c->SetNCells(ncells_true);
 893                 c->SetDispersion(recpoint->GetDispersion());
 894                 c->SetEmcCpvDistance(-1);            //not yet implemented
 895                 c->SetChi2(-1);                      //not yet implemented
 896                 c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 897                 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 898                 Float_t elipAxis[2];
 899                 recpoint->GetElipsAxis(elipAxis);
 900                 c->SetM02(elipAxis[0]*elipAxis[0]) ;
 901                 c->SetM20(elipAxis[1]*elipAxis[1]) ;
 902
 903                 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
 904                 cesd->SetCellsAbsId(absIds);
 905                 cesd->SetCellsAmplitudeFraction(ratios);
 906         }
 907
 908 }
 909