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