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