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