PWG/EMCAL/AliAnalysisTaskEMCALClusterizeFast.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 // This analysis provides a new list of clusters to be used in other analysis
  18 //
  19 // Author: Constantin Loizides, Salvatore Aiola
  20 //         Adapted from analysis class from Deepa Thomas
  21 //
  22 // $Id$
  23 //_________________________________________________________________________
  24
  25 // --- Root ---
  26 #include <TClonesArray.h>
  27 #include <TGeoManager.h>
  28 #include <TObjArray.h>
  29 #include <TString.h>
  30 #include <TTree.h>
  31 #include <TArrayI.h>
  32
  33 // --- AliRoot ---
  34 #include "AliAODCaloCluster.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAnalysisManager.h"
  37 #include "AliCDBEntry.h"
  38 #include "AliCDBManager.h"
  39 #include "AliCaloCalibPedestal.h"
  40 #include "AliEMCALAfterBurnerUF.h"
  41 #include "AliEMCALCalibData.h"
  42 #include "AliEMCALClusterizerNxN.h"
  43 #include "AliEMCALClusterizerv1.h"
  44 #include "AliEMCALClusterizerv2.h"
  45 #include "AliEMCALClusterizerFixedWindow.h"
  46 #include "AliEMCALDigit.h"
  47 #include "AliEMCALGeometry.h"
  48 #include "AliEMCALRecParam.h"
  49 #include "AliEMCALRecPoint.h"
  50 #include "AliEMCALRecoUtils.h"
  51 #include "AliESDEvent.h"
  52 #include "AliInputEventHandler.h"
  53 #include "AliLog.h"
  54
  55 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  56
  57 ClassImp(AliAnalysisTaskEMCALClusterizeFast)
  58
  59 //________________________________________________________________________
  60 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast() :
  61   AliAnalysisTaskSE(),
  62   fRun(-1),
  63   fDigitsArr(0),
  64   fClusterArr(0),
  65   fRecParam(new AliEMCALRecParam),
  66   fClusterizer(0),
  67   fUnfolder(0),
  68   fJustUnfold(kFALSE),
  69   fGeomName(),
  70   fGeomMatrixSet(kFALSE),
  71   fLoadGeomMatrices(kFALSE),
  72   fOCDBpath(),
  73   fCalibData(0),
  74   fPedestalData(0),
  75   fOutputAODBranch(0),
  76   fOutputAODBrName(),
  77   fRecoUtils(0),
  78   fLoadCalib(kFALSE),
  79   fLoadPed(kFALSE),
  80   fAttachClusters(kTRUE),
  81   fSubBackground(kFALSE),
  82   fNPhi(4),
  83   fNEta(4),
  84   fShiftPhi(2),
  85   fShiftEta(2),
  86   fTRUShift(0),
  87   fInputCellType(kFEEData),
  88   fTrackName(),
  89   fCaloCellsName(),
  90   fCaloClustersName("newCaloClusters"),
  91   fDoUpdateCells(kFALSE),
  92   fDoClusterize(kTRUE),
  93   fClusterBadChannelCheck(kFALSE),
  94   fRejectExoticClusters(kFALSE),
  95   fRejectExoticCells(kFALSE),
  96   fFiducial(kFALSE),
  97   fDoNonLinearity(kFALSE),
  98   fRecalDistToBadChannels(kFALSE),
  99   fCaloCells(0),
 100   fCaloClusters(0),
 101   fEsd(0),
 102   fAod(0),
 103   fGeom(0)
 104 {
 105   // Constructor
 106
 107   for(Int_t i = 0; i < 12; ++i)
 108     fGeomMatrix[i] = 0;
 109 }
 110
 111 //________________________________________________________________________
 112 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name) :
 113   AliAnalysisTaskSE(name),
 114   fRun(-1),
 115   fDigitsArr(0),
 116   fClusterArr(0),
 117   fRecParam(new AliEMCALRecParam),
 118   fClusterizer(0),
 119   fUnfolder(0),
 120   fJustUnfold(kFALSE),
 121   fGeomName(),
 122   fGeomMatrixSet(kFALSE),
 123   fLoadGeomMatrices(kFALSE),
 124   fOCDBpath(),
 125   fCalibData(0),
 126   fPedestalData(0),
 127   fOutputAODBranch(0),
 128   fOutputAODBrName(),
 129   fRecoUtils(0),
 130   fLoadCalib(kFALSE),
 131   fLoadPed(kFALSE),
 132   fAttachClusters(kTRUE),
 133   fSubBackground(kFALSE),
 134   fNPhi(4),
 135   fNEta(4),
 136   fShiftPhi(2),
 137   fShiftEta(2),
 138   fTRUShift(0),
 139   fInputCellType(kFEEData),
 140   fTrackName(),
 141   fCaloCellsName(),
 142   fCaloClustersName("newCaloClusters"),
 143   fDoUpdateCells(kFALSE),
 144   fDoClusterize(kTRUE),
 145   fClusterBadChannelCheck(kFALSE),
 146   fRejectExoticClusters(kFALSE),
 147   fRejectExoticCells(kFALSE),
 148   fFiducial(kFALSE),
 149   fDoNonLinearity(kFALSE),
 150   fRecalDistToBadChannels(kFALSE),
 151   fCaloCells(0),
 152   fCaloClusters(0),
 153   fEsd(0),
 154   fAod(0),
 155   fGeom(0)
 156 {
 157   // Constructor
 158
 159   fBranchNames     = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
 160   for(Int_t i = 0; i < 12; ++i)
 161     fGeomMatrix[i] = 0;
 162 }
 163
 164 //________________________________________________________________________
 165 AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
 166 {
 167   // Destructor.
 168
 169   delete fClusterizer;
 170   delete fUnfolder;
 171   delete fRecoUtils;
 172   delete fRecParam;
 173 }
 174
 175 //________________________________________________________________________
 176 void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
 177 {
 178   // Create output objects.
 179
 180   if (!fOutputAODBrName.IsNull()) {
 181     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 182     fOutputAODBranch->SetName(fOutputAODBrName);
 183     AddAODBranch("TClonesArray", &fOutputAODBranch);
 184     AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
 185   }
 186 }
 187
 188 //________________________________________________________________________
 189 void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
 190 {
 191   // Main loop, called for each event
 192
 193   // if no operation is requested, return
 194   if (!fDoClusterize && !fDoUpdateCells)
 195     return;
 196
 197   // remove the contents of output list set in the previous event
 198   if (fOutputAODBranch)
 199     fOutputAODBranch->Clear("C");
 200
 201   fEsd = dynamic_cast<AliESDEvent*>(InputEvent());
 202   fAod = dynamic_cast<AliAODEvent*>(InputEvent());
 203
 204   if (!fEsd&&!fAod) {
 205     Error("UserExec","Event not available");
 206     return;
 207   }
 208
 209   LoadBranches();
 210
 211   UInt_t offtrigger = 0;
 212   if (fEsd) {
 213     UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
 214     UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
 215     Bool_t desc1 = (mask1 >> 18) & 0x1;
 216     Bool_t desc2 = (mask2 >> 18) & 0x1;
 217     if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
 218       AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
 219                     mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
 220                     mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
 221       return;
 222     }
 223     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 224     offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
 225   } else {
 226     offtrigger =  fAod->GetHeader()->GetOfflineTrigger();
 227   }
 228
 229   if (!MCEvent()) {
 230     if (offtrigger & AliVEvent::kFastOnly) {
 231       AliError(Form("EMCAL not in fast only partition"));
 232       return;
 233     }
 234   }
 235
 236   Init();
 237
 238   if (fJustUnfold) {
 239     AliWarning("Unfolding not implemented");
 240     return;
 241   }
 242
 243   FillDigitsArray();
 244
 245   if (fDoClusterize)
 246     Clusterize();
 247
 248   if (fDoUpdateCells)
 249     UpdateCells();
 250
 251   if (!fDoClusterize || (!fAttachClusters && !fOutputAODBranch) || !fCaloClusters)
 252     return;
 253
 254   UpdateClusters();
 255   CalibrateClusters();
 256
 257   if (fOutputAODBranch && fCaloClusters != fOutputAODBranch)
 258     CopyClusters(fCaloClusters, fOutputAODBranch);
 259 }
 260
 261 //________________________________________________________________________
 262 void AliAnalysisTaskEMCALClusterizeFast::CopyClusters(TClonesArray *orig, TClonesArray *dest)
 263 {
 264   const Int_t Ncls = orig->GetEntries();
 265
 266   for(Int_t i=0; i < Ncls; ++i) {
 267     AliVCluster *oc = static_cast<AliVCluster*>(orig->At(i));
 268
 269     if (!oc)
 270       continue;
 271
 272     if (!oc->IsEMCAL())
 273       continue;
 274
 275     AliVCluster *dc = static_cast<AliVCluster*>(dest->New(i));
 276     dc->SetType(AliVCluster::kEMCALClusterv1);
 277     dc->SetE(oc->E());
 278     Float_t pos[3] = {0};
 279     oc->GetPosition(pos);
 280     dc->SetPosition(pos);
 281     dc->SetNCells(oc->GetNCells());
 282     dc->SetCellsAbsId(oc->GetCellsAbsId());
 283     dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
 284     dc->SetID(oc->GetID());
 285     dc->SetDispersion(oc->GetDispersion());
 286     dc->SetEmcCpvDistance(-1);
 287     dc->SetChi2(-1);
 288     dc->SetTOF(oc->GetTOF());     //time-of-flight
 289     dc->SetNExMax(oc->GetNExMax()); //number of local maxima
 290     dc->SetM02(oc->GetM02());
 291     dc->SetM20(oc->GetM20());
 292     dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
 293     dc->SetMCEnergyFraction(oc->GetMCEnergyFraction());
 294
 295     //MC
 296     UInt_t nlabels = oc->GetNLabels();
 297     Int_t *labels = oc->GetLabels();
 298
 299     if (nlabels == 0 || !labels)
 300       continue;
 301
 302     AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
 303     if (esdClus) {
 304       TArrayI parents(nlabels, labels);
 305       esdClus->AddLabels(parents);
 306     }
 307     else {
 308       AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
 309       if (aodClus)
 310         aodClus->SetLabel(labels, nlabels);
 311     }
 312   }
 313 }
 314
 315 //________________________________________________________________________
 316 void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
 317 {
 318   // Clusterize
 319
 320   if (fSubBackground) {
 321     fClusterizer->SetInputCalibrated(kTRUE);
 322     fClusterizer->SetCalibrationParameters(0);
 323   }
 324
 325   fClusterizer->Digits2Clusters("");
 326
 327   if (fSubBackground) {
 328     if (fCalibData) {
 329       fClusterizer->SetInputCalibrated(kFALSE);
 330       fClusterizer->SetCalibrationParameters(fCalibData);
 331     }
 332   }
 333 }
 334
 335 //________________________________________________________________________
 336 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
 337 {
 338   // Fill digits array
 339
 340   fDigitsArr->Clear("C");
 341   switch (fInputCellType) {
 342
 343   case kFEEData :
 344   case kFEEDataMCOnly :
 345   case kFEEDataExcludeMC :
 346     {
 347       Double_t avgE        = 0; // for background subtraction
 348       const Int_t ncells   = fCaloCells->GetNumberOfCells();
 349       for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 350         Double_t cellAmplitude=0, cellTime=0, cellEFrac = 0;
 351         Short_t  cellNumber=0;
 352         Int_t cellMCLabel=-1;
 353         if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
 354           break;
 355
 356         if (cellMCLabel > 0 && cellEFrac < 1e-6)
 357           cellEFrac = 1;
 358
 359         if (cellAmplitude < 1e-6 || cellNumber < 0)
 360           continue;
 361
 362         if (fInputCellType == kFEEDataMCOnly) {
 363           if (cellMCLabel <= 0)
 364             continue;
 365           else {
 366             cellAmplitude *= cellEFrac;
 367             cellEFrac = 1;
 368           }
 369         }
 370         else if (fInputCellType == kFEEDataExcludeMC) {
 371           if (cellMCLabel > 0)
 372             continue;
 373           else {
 374             cellAmplitude *= 1 - cellEFrac;
 375             cellEFrac = 0;
 376           }
 377         }
 378
 379         if(!AcceptCell(cellNumber)) continue;
 380
 381         AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
 382                                                                         (Float_t)cellAmplitude, (Float_t)cellTime,
 383                                                                         AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
 384
 385         if (!fDoClusterize||fSubBackground) {
 386           Float_t energy = cellAmplitude;
 387           Float_t time   = cellTime;
 388           fClusterizer->Calibrate(energy,time,cellNumber);
 389           digit->SetAmplitude(energy);
 390           avgE += energy;
 391         }
 392         idigit++;
 393       }
 394
 395       if (fSubBackground) {
 396         avgE /= fGeom->GetNumberOfSuperModules()*48*24;
 397         Int_t ndigis = fDigitsArr->GetEntries();
 398         for (Int_t i = 0; i < ndigis; ++i) {
 399           AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
 400           Double_t energy = digit->GetAmplitude() - avgE;
 401           if (energy<=0.001) {
 402             digit->SetAmplitude(0);
 403           } else {
 404             digit->SetAmplitude(energy);
 405           }
 406         }
 407       }
 408     }
 409     break;
 410
 411   case kPattern :
 412     {
 413       // Fill digits from a pattern
 414       Int_t maxd = fGeom->GetNCells() / 4;
 415       for (Int_t idigit = 0; idigit < maxd; idigit++){
 416         if (idigit % 24 == 12) idigit += 12;
 417         AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 418         digit->SetId(idigit * 4);
 419         digit->SetTime(600);
 420         digit->SetTimeR(600);
 421         digit->SetIndexInList(idigit);
 422         digit->SetType(AliEMCALDigit::kHG);
 423         digit->SetAmplitude(0.1);
 424       }
 425     }
 426     break;
 427
 428   case kL0FastORs    :
 429   case kL0FastORsTC  :
 430   case kL1FastORs    :
 431     {
 432       // Fill digits from FastORs
 433
 434       AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
 435
 436       if (!triggers || !(triggers->GetEntries() > 0))
 437         return;
 438
 439       Int_t idigit = 0;
 440       triggers->Reset();
 441
 442       while ((triggers->Next())) {
 443         Float_t L0Amplitude = 0;
 444         triggers->GetAmplitude(L0Amplitude);
 445
 446         if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
 447           continue;
 448
 449         Int_t L1Amplitude = 0;
 450         triggers->GetL1TimeSum(L1Amplitude);
 451
 452         if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
 453           continue;
 454
 455         Int_t triggerTime = 0;
 456         Int_t ntimes = 0;
 457         triggers->GetNL0Times(ntimes);
 458
 459         if (ntimes < 1 && fInputCellType == kL0FastORsTC)
 460           continue;
 461
 462         if (ntimes > 0) {
 463           Int_t trgtimes[25];
 464           triggers->GetL0Times(trgtimes);
 465           triggerTime = trgtimes[0];
 466         }
 467
 468         Int_t triggerCol = 0, triggerRow = 0;
 469         triggers->GetPosition(triggerCol, triggerRow);
 470
 471         Int_t find = -1;
 472         fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
 473
 474         if (find < 0)
 475           continue;
 476
 477         Int_t cidx[4] = {-1};
 478         Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
 479
 480         if (!ret)
 481           continue;
 482
 483         Float_t triggerAmplitude = 0;
 484
 485         if (fInputCellType == kL1FastORs) {
 486           triggerAmplitude = 0.25 * L1Amplitude;  // it will add 4 cells for 1 amplitude
 487         }
 488         else {
 489           triggerAmplitude = L0Amplitude;      // 10 bit truncated, so it is already divided by 4
 490         }
 491
 492         for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
 493           Int_t triggerNumber = cidx[idxpos];
 494           AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 495           digit->SetId(triggerNumber);
 496           digit->SetTime(triggerTime);
 497           digit->SetTimeR(triggerTime);
 498           digit->SetIndexInList(idigit);
 499           digit->SetType(AliEMCALDigit::kHG);
 500           digit->SetAmplitude(triggerAmplitude);
 501           idigit++;
 502         }
 503       }
 504     }
 505     break;
 506   }
 507 }
 508
 509 //________________________________________________________________________________________
 510 Bool_t AliAnalysisTaskEMCALClusterizeFast::AcceptCell(Int_t cellNumber) {
 511
 512   Bool_t accept = kTRUE;
 513   if(fRejectExoticCells) {
 514     //Remove exotic cells before making digits
 515     Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
 516     fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
 517     Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
 518     Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
 519     if(isEx) accept = kFALSE;
 520     if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
 521   }
 522   return accept;
 523 }
 524
 525 //________________________________________________________________________________________
 526 void AliAnalysisTaskEMCALClusterizeFast::CalibrateClusters()
 527 {
 528   // Go through clusters one by one and process separate correction
 529   // as those were defined or not
 530
 531   Int_t nclusters = fCaloClusters->GetEntriesFast();
 532   for (Int_t icluster=0; icluster < nclusters; ++icluster) {
 533     AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
 534     if (!clust)
 535       continue;
 536     if (!clust->IsEMCAL())
 537       continue;
 538
 539     // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
 540     if (fClusterBadChannelCheck) {
 541       // careful, the the ClusterContainsBadChannel is dependent on
 542       // SwitchOnBadChannelsRemoval, switching it ON automatically
 543       // and returning to original value after processing
 544       Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
 545       fRecoUtils->SwitchOnBadChannelsRemoval();
 546
 547       Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
 548
 549       // switch the bad channels removal back
 550       if (!badRemoval)
 551         fRecoUtils->SwitchOffBadChannelsRemoval();
 552
 553       if (badResult) {
 554         delete fCaloClusters->RemoveAt(icluster);
 555         continue; //TODO is it really needed to remove it? Or should we flag it?
 556       }
 557     }
 558
 559     // REMOVE EXOTIC CLUSTERS -------------------------------------
 560     // does process local cell recalibration energy and time without replacing
 561     // the global cell values, in case of no cell recalib done yet
 562     if (fRejectExoticClusters) {
 563       // careful, the IsExoticCluster is dependent on
 564       // SwitchOnRejectExoticCell, switching it ON automatically
 565       // and returning to original value after processing
 566       Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
 567       fRecoUtils->SwitchOnRejectExoticCell();
 568
 569       // get bunch crossing
 570       Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
 571
 572       Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
 573
 574       // switch the exotic channels removal back
 575       if (!exRemoval)
 576         fRecoUtils->SwitchOffRejectExoticCell();
 577
 578       if (exResult) {
 579         delete fCaloClusters->RemoveAt(icluster);
 580         continue; //TODO is it really needed to remove it? Or should we flag it?
 581       }
 582     }
 583
 584     // FIDUCIAL CUT -----------------------------------------------
 585     if (fFiducial) {
 586       // depends on SetNumberOfCellsFromEMCALBorder
 587       // SwitchOnNoFiducialBorderInEMCALEta0
 588       if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
 589         delete fCaloClusters->RemoveAt(icluster);
 590         continue; //TODO it would be nice to store the distance
 591       }
 592     }
 593
 594     // NONLINEARITY -----------------------------------------------
 595     if (fDoNonLinearity) {
 596       Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
 597       clust->SetE(correctedEnergy);
 598     }
 599
 600     // DISTANCE TO BAD CHANNELS -----------------------------------
 601     if (fRecalDistToBadChannels)
 602       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
 603   }
 604
 605   fCaloClusters->Compress();
 606 }
 607
 608 //________________________________________________________________________________________
 609 void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
 610 {
 611   Float_t g[3]={0};
 612   c->GetPosition(g);
 613   TVector3 gpos(g);
 614
 615   Double_t dEtaMin  = 1e9;
 616   Double_t dPhiMin  = 1e9;
 617   Int_t    imin     = -1;
 618   Double_t ceta     = gpos.Eta();
 619   Double_t cphi     = gpos.Phi();
 620   const Int_t ntrks = tarr->GetEntries();
 621   for(Int_t t = 0; t<ntrks; ++t) {
 622     AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
 623     if (!track)
 624       continue;
 625     const AliExternalTrackParam *outp = track->GetOuterParam();
 626     if (!outp)
 627       continue;
 628     Double_t trkPos[3] = {0.,0.,0.};
 629     if (!outp->GetXYZ(trkPos))
 630       continue;
 631     TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
 632     Double_t veta = vec.Eta();
 633     Double_t vphi = vec.Phi();
 634     if(vphi<0)
 635       vphi += 2*TMath::Pi();
 636     if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
 637       continue;
 638     Double_t dR = vec.DeltaR(gpos);
 639     if(dR > 25)
 640       continue;
 641     Float_t tmpEta=0, tmpPhi=0;
 642     if (0) {
 643       AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
 644       Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
 645       if (!ret)
 646         continue;
 647     } else {
 648       tmpEta = ceta - veta;
 649       tmpPhi = cphi - vphi;
 650     }
 651     if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
 652       dEtaMin = tmpEta;
 653       dPhiMin = tmpPhi;
 654       imin = t;
 655     }
 656   }
 657   c->SetEmcCpvDistance(imin);
 658   c->SetTrackDistance(dPhiMin, dEtaMin);
 659 }
 660
 661 //________________________________________________________________________________________
 662 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
 663 {
 664   // Cluster energy, global position, cells and their amplitude fractions are restored.
 665
 666   // tracks array for track/cluster matching
 667   TClonesArray *tarr = 0;
 668   if (!fTrackName.IsNull()) {
 669     tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
 670     if (!tarr) {
 671       AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
 672     }
 673   }
 674
 675   const Int_t Ncls = fClusterArr->GetEntries();
 676   AliDebug(1, Form("total no of clusters %d", Ncls));
 677   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
 678     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 679     Int_t ncells_true = 0;
 680     const Int_t ncells = recpoint->GetMultiplicity();
 681     UShort_t   absIds[ncells];
 682     Double32_t ratios[ncells];
 683     Int_t *dlist = recpoint->GetDigitsList();
 684     Float_t *elist = recpoint->GetEnergiesList();
 685     Double_t mcEnergy = 0;
 686     for (Int_t c = 0; c < ncells; ++c) {
 687       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 688       absIds[ncells_true] = digit->GetId();
 689       ratios[ncells_true] = elist[c]/recpoint->GetEnergy();
 690       if (digit->GetIparent(1) > 0)
 691         mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
 692       ++ncells_true;
 693     }
 694
 695     if (ncells_true < 1) {
 696       AliWarning("Skipping cluster with no cells");
 697       continue;
 698     }
 699
 700     // calculate new cluster position
 701     TVector3 gpos;
 702     recpoint->GetGlobalPosition(gpos);
 703     Float_t g[3];
 704     gpos.GetXYZ(g);
 705
 706     AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
 707
 708     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
 709     c->SetType(AliVCluster::kEMCALClusterv1);
 710     c->SetE(recpoint->GetEnergy());
 711     c->SetPosition(g);
 712     c->SetNCells(ncells_true);
 713     c->SetCellsAbsId(absIds);
 714     c->SetCellsAmplitudeFraction(ratios);
 715     c->SetID(recpoint->GetUniqueID());
 716     c->SetDispersion(recpoint->GetDispersion());
 717     c->SetEmcCpvDistance(-1);
 718     c->SetChi2(-1);
 719     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 720     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 721     Float_t elipAxis[2];
 722     recpoint->GetElipsAxis(elipAxis);
 723     c->SetM02(elipAxis[0]*elipAxis[0]);
 724     c->SetM20(elipAxis[1]*elipAxis[1]);
 725     c->SetMCEnergyFraction(mcEnergy);
 726
 727     //MC
 728     AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(c);
 729     if (esdClus) {
 730       Int_t  parentMult = 0;
 731       Int_t *parentList = recpoint->GetParents(parentMult);
 732       if (parentMult > 0) {
 733         TArrayI parents(parentMult, parentList);
 734         esdClus->AddLabels(parents);
 735       }
 736     }
 737     else {
 738       AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(c);
 739       if (aodClus) {
 740         Int_t  parentMult = 0;
 741         Int_t *parentList = recpoint->GetParents(parentMult);
 742         aodClus->SetLabel(parentList, parentMult);
 743       }
 744     }
 745
 746     if (tarr)
 747       TrackClusterMatching(c, tarr);
 748   }
 749 }
 750
 751 //________________________________________________________________________
 752 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
 753 {
 754   // Update cells in case re-calibration was done.
 755   if (!fCalibData&&!fSubBackground)
 756     return;
 757
 758   const Int_t   ncells = fCaloCells->GetNumberOfCells();
 759   const Int_t   ndigis = fDigitsArr->GetEntries();
 760   if (ncells!=ndigis) {
 761     fCaloCells->DeleteContainer();
 762     fCaloCells->CreateContainer(ndigis);
 763   }
 764   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
 765     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
 766     Double_t cellAmplitude = digit->GetCalibAmp();
 767     Short_t cellNumber     = digit->GetId();
 768     Double_t cellTime      = digit->GetTime();
 769     fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
 770   }
 771 }
 772
 773 //________________________________________________________________________
 774 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
 775 {
 776   // Update cells in case re-calibration was done.
 777
 778   const Int_t nents = fCaloClusters->GetEntries();
 779   for (Int_t i=0;i<nents;++i) {
 780     AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
 781     if (!c)
 782       continue;
 783     if (c->IsEMCAL())
 784       delete fCaloClusters->RemoveAt(i);
 785   }
 786
 787   fCaloClusters->Compress();
 788
 789   RecPoints2Clusters(fCaloClusters);
 790 }
 791
 792 //________________________________________________________________________________________
 793 void AliAnalysisTaskEMCALClusterizeFast::Init()
 794 {
 795   // Select clusterization/unfolding algorithm and set all the needed parameters.
 796
 797   if (InputEvent()->GetRunNumber()==fRun)
 798     return;
 799   fRun = InputEvent()->GetRunNumber();
 800
 801   if (fJustUnfold){
 802     // init the unfolding afterburner
 803     delete fUnfolder;
 804     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
 805     return;
 806   }
 807
 808   if (fGeomName.Length()>0)
 809     fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 810   else
 811     fGeom = AliEMCALGeometry::GetInstance();
 812   if (!fGeom) {
 813     AliFatal("Geometry not available!!!");
 814     return;
 815   }
 816
 817   if (!fGeomMatrixSet) {
 818     if (fLoadGeomMatrices) {
 819       for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
 820         if (fGeomMatrix[mod]){
 821           if (DebugLevel() > 2)
 822             fGeomMatrix[mod]->Print();
 823           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
 824         }
 825       }
 826     } else { // get matrix from file (work around bug in aliroot)
 827       for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
 828         const TGeoHMatrix *gm = 0;
 829         if (fEsd) {
 830           gm = fEsd->GetEMCALMatrix(mod);
 831         } else {
 832           AliAODHeader *aodheader = fAod->GetHeader();
 833           if (aodheader) {
 834             gm = aodheader->GetEMCALMatrix(mod);
 835           }
 836         }
 837         if (gm) {
 838           if (DebugLevel() > 2)
 839             gm->Print();
 840           fGeom->SetMisalMatrix(gm,mod);
 841         }
 842       }
 843     }
 844     fGeomMatrixSet=kTRUE;
 845   }
 846
 847   // setup digit array if needed
 848   if (!fDigitsArr) {
 849     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
 850     fDigitsArr->SetOwner(1);
 851   }
 852
 853   // then setup clusterizer
 854   if (fClusterizer) {
 855     // avoid to delete digits array
 856     fClusterizer->SetDigitsArr(0);
 857     delete fClusterizer;
 858   }
 859   if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 860     fClusterizer = new AliEMCALClusterizerv1(fGeom);
 861   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
 862     AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
 863     clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
 864     clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
 865     fClusterizer = clusterizer;
 866   }
 867   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 868     fClusterizer = new AliEMCALClusterizerv2(fGeom);
 869   else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
 870     AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
 871     clusterizer->SetNphi(fNPhi);
 872     clusterizer->SetNeta(fNEta);
 873     clusterizer->SetShiftPhi(fShiftPhi);
 874     clusterizer->SetShiftEta(fShiftEta);
 875     clusterizer->SetTRUshift(fTRUShift);
 876     fClusterizer = clusterizer;
 877   }
 878   else {
 879     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 880   }
 881   fClusterizer->InitParameters(fRecParam);
 882
 883   if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
 884     AliCDBManager *cdb = AliCDBManager::Instance();
 885     if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
 886       cdb->SetDefaultStorage(fOCDBpath);
 887     if (fRun!=cdb->GetRun())
 888       cdb->SetRun(fRun);
 889   }
 890   if (!fCalibData&&fLoadCalib&&fRun>0) {
 891     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
 892     if (entry)
 893       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
 894     if (!fCalibData)
 895       AliFatal("Calibration parameters not found in CDB!");
 896   }
 897   if (!fPedestalData&&fLoadPed&&fRun>0) {
 898     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
 899     if (entry)
 900       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
 901   }
 902   if (fCalibData) {
 903     fClusterizer->SetInputCalibrated(kFALSE);
 904     fClusterizer->SetCalibrationParameters(fCalibData);
 905   } else {
 906     fClusterizer->SetInputCalibrated(kTRUE);
 907   }
 908   fClusterizer->SetCaloCalibPedestal(fPedestalData);
 909   fClusterizer->SetJustClusters(kTRUE);
 910   fClusterizer->SetDigitsArr(fDigitsArr);
 911   fClusterizer->SetOutput(0);
 912   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 913
 914   // Get the emcal cells
 915   if ((fInputCellType == kFEEData ||  fInputCellType == kFEEDataMCOnly || fInputCellType == kFEEDataExcludeMC) && !fCaloCells) {
 916     if (fCaloCellsName.IsNull()) {
 917       fCaloCells = InputEvent()->GetEMCALCells();
 918     }
 919     else {
 920       fCaloCells =  dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
 921       if (!fCaloCells)
 922         AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
 923     }
 924     if (!fCaloCells)
 925       AliFatal("Could not get EMCal cells!");
 926   }
 927
 928   // Set output clusters collection
 929   if (!fAttachClusters) {
 930     fCaloClusters = fOutputAODBranch;
 931     return;
 932   }
 933
 934   if (!fCaloClusters) {
 935     if (fCaloClustersName.IsNull()) { //overwrite mode
 936       if (fEsd)
 937         fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
 938       else if (fAod)
 939         fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
 940     }
 941     else {
 942       fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
 943
 944       if (!fCaloClusters) {
 945         if (fEsd)
 946           fCaloClusters = new TClonesArray("AliESDCaloCluster");
 947         else if (fAod)
 948           fCaloClusters = new TClonesArray("AliAODCaloCluster");
 949
 950         fCaloClusters->SetName(fCaloClustersName);
 951         InputEvent()->AddObject(fCaloClusters);
 952       }
 953     }
 954
 955     if (!fCaloClusters)
 956       AliFatal("Could not get cluster collection!");
 957
 958     TClass *cl = fCaloClusters->GetClass();
 959     if (!cl->GetBaseClass("AliVCluster")) {
 960       AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
 961       fCaloClusters = 0;
 962       return;
 963     }
 964   }
 965 }