PWG4/CaloCalib/AliAnalysisTaskEMCALClusterize.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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // This analysis provides a new list of clusters to be used in other analysis
  19 //
  20 // Author: Gustavo Conesa Balbastre,
  21 //         Adapted from analysis class from Deepa Thomas
  22 //
  23 //
  24 //_________________________________________________________________________
  25
  26 // --- Root ---
  27 #include "TString.h"
  28 #include "TRefArray.h"
  29 #include "TClonesArray.h"
  30 #include "TTree.h"
  31 #include "TGeoManager.h"
  32 #include "TROOT.h"
  33 #include "TInterpreter.h"
  34 #include "TFile.h"
  35 //#include "string.h"
  36
  37 // --- AliRoot Analysis Steering
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliESDEvent.h"
  41 #include "AliGeomManager.h"
  42 #include "AliVCaloCells.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliCDBManager.h"
  45 #include "AliCDBStorage.h"
  46 #include "AliCDBEntry.h"
  47 #include "AliLog.h"
  48 #include "AliVEventHandler.h"
  49 #include "AliAODInputHandler.h"
  50
  51 // --- EMCAL
  52 #include "AliEMCALRecParam.h"
  53 #include "AliEMCALAfterBurnerUF.h"
  54 #include "AliEMCALGeometry.h"
  55 #include "AliEMCALClusterizerNxN.h"
  56 #include "AliEMCALClusterizerv1.h"
  57 #include "AliEMCALClusterizerv2.h"
  58 #include "AliEMCALRecPoint.h"
  59 #include "AliEMCALDigit.h"
  60 #include "AliCaloCalibPedestal.h"
  61 #include "AliEMCALCalibData.h"
  62 #include "AliEMCALRecoUtils.h"
  63
  64 #include "AliAnalysisTaskEMCALClusterize.h"
  65
  66 ClassImp(AliAnalysisTaskEMCALClusterize)
  67
  68 //________________________________________________________________________
  69 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  70   : AliAnalysisTaskSE(name)
  71   , fGeom(0), fGeomName("EMCAL_COMPLETEV1"),   fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  72   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
  73   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  74   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  75   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
  76   , fFillAODFile(kTRUE), fFillAODHeader(0),    fFillAODCaloCells(0)
  77   , fRun(-1),            fRecoUtils(0),        fConfigName("")
  78   , fCellLabels(),       fCellSecondLabels(),  fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
  79
  80 {
  81   //ctor
  82   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
  83   for(Int_t j = 0; j < 24*48*11; j++)  {
  84     fCellLabels[j]       = -1;
  85     fCellSecondLabels[j] = -1;
  86   }
  87   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  88   fClusterArr      = new TObjArray(100);
  89   fCaloClusterArr  = new TObjArray(100);
  90   fRecParam        = new AliEMCALRecParam;
  91   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  92   fRecoUtils       = new AliEMCALRecoUtils();
  93 }
  94
  95 //________________________________________________________________________
  96 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
  97   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
  98   , fGeom(0), fGeomName("EMCAL_COMPLETEV1"),    fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  99   , fCalibData(0),        fPedestalData(0),     fOCDBpath("raw://"),    fAccessOCDB(kFALSE)
 100   , fDigitsArr(0),        fClusterArr(0),       fCaloClusterArr(0)
 101   , fRecParam(0),         fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
 102   , fOutputAODBranch(0),  fOutputAODBranchName("newEMCALClusters")
 103   , fFillAODFile(kFALSE), fFillAODHeader(0),    fFillAODCaloCells(0)
 104   , fRun(-1),             fRecoUtils(0),        fConfigName("")
 105   , fCellLabels(),        fCellSecondLabels(),  fMaxEvent(1000000000),  fDoTrackMatching(kFALSE)
 106 {
 107   // Constructor
 108   for(Int_t i = 0; i < 10;    i++)  fGeomMatrix[i] =  0;
 109   for(Int_t j = 0; j < 24*48*11; j++)  {
 110     fCellLabels[j]       = -1;
 111     fCellSecondLabels[j] = -1;
 112   }
 113   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
 114   fClusterArr      = new TObjArray(100);
 115   fCaloClusterArr  = new TObjArray(100);
 116   fRecParam        = new AliEMCALRecParam;
 117   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
 118   fRecoUtils       = new AliEMCALRecoUtils();
 119 }
 120
 121
 122 //________________________________________________________________________
 123 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 124 {
 125   //dtor
 126
 127   if (fDigitsArr){
 128     fDigitsArr->Clear("C");
 129     delete fDigitsArr;
 130   }
 131
 132   if (fClusterArr){
 133     fClusterArr->Delete();
 134     delete fClusterArr;
 135   }
 136
 137   if (fCaloClusterArr){
 138     fCaloClusterArr->Delete();
 139     delete fCaloClusterArr;
 140   }
 141
 142   if(fClusterizer) delete fClusterizer;
 143   if(fUnfolder)    delete fUnfolder;
 144   if(fRecoUtils)   delete fRecoUtils;
 145
 146 }
 147
 148 //-------------------------------------------------------------------
 149 void AliAnalysisTaskEMCALClusterize::Init()
 150 {
 151   //Init analysis with configuration macro if available
 152
 153   if(gROOT->LoadMacro(fConfigName) >=0){
 154     printf("Configure analysis with %s\n",fConfigName.Data());
 155     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 156     fGeomName         = clus->fGeomName;
 157     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 158     fOCDBpath         = clus->fOCDBpath;
 159     fAccessOCDB       = clus->fAccessOCDB;
 160     fRecParam         = clus->fRecParam;
 161     fJustUnfold       = clus->fJustUnfold;
 162     fFillAODFile      = clus->fFillAODFile;
 163     fRecoUtils        = clus->fRecoUtils;
 164     fConfigName       = clus->fConfigName;
 165     fMaxEvent         = clus->fMaxEvent;
 166     fDoTrackMatching  = clus->fDoTrackMatching;
 167     fOutputAODBranchName = clus->fOutputAODBranchName;
 168     for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
 169
 170   }
 171
 172 }
 173
 174 //-------------------------------------------------------------------
 175 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
 176 {
 177   // Init geometry, create list of output clusters
 178
 179   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;
 180   if(fOutputAODBranchName.Length()!=0){
 181     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 182     fOutputAODBranch->SetName(fOutputAODBranchName);
 183     AddAODBranch("TClonesArray", &fOutputAODBranch);
 184   }
 185   else {
 186     AliFatal("fOutputAODBranchName not set\n");
 187   }
 188 }
 189
 190 //________________________________________________________________________
 191 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
 192   // Put calo cells in standard branch
 193   AliVEvent * event = InputEvent();
 194   AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
 195   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 196
 197   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 198   aodEMcells.CreateContainer(nEMcell);
 199   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 200   Double_t calibFactor = 1.;
 201   for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 202     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 203     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 204     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 205
 206     if(fRecoUtils->IsRecalibrationOn()){
 207       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 208     }
 209
 210     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 211       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 212     }
 213     else {
 214       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 215     }
 216   }
 217   aodEMcells.Sort();
 218
 219 }
 220
 221 //________________________________________________________________________
 222 void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
 223   //Put event header information in standard AOD branch
 224
 225   AliVEvent* event      = InputEvent();
 226   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
 227   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
 228
 229   Double_t pos[3]   ;
 230   Double_t covVtx[6];
 231   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 232
 233   AliAODHeader* header = AODEvent()->GetHeader();
 234   header->SetRunNumber(event->GetRunNumber());
 235
 236   if(esdevent){
 237     TTree* tree = fInputHandler->GetTree();
 238     if (tree) {
 239       TFile* file = tree->GetCurrentFile();
 240       if (file) header->SetESDFileName(file->GetName());
 241     }
 242   }
 243   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 244
 245   header->SetBunchCrossNumber(event->GetBunchCrossNumber());
 246   header->SetOrbitNumber(event->GetOrbitNumber());
 247   header->SetPeriodNumber(event->GetPeriodNumber());
 248   header->SetEventType(event->GetEventType());
 249
 250   //Centrality
 251   if(event->GetCentrality()){
 252     header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
 253   }
 254   else{
 255     header->SetCentrality(0);
 256   }
 257
 258   //Trigger
 259   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 260   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 261   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 262   header->SetTriggerMask(event->GetTriggerMask());
 263   header->SetTriggerCluster(event->GetTriggerCluster());
 264   if(esdevent){
 265     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 266     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 267     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 268   }
 269   else if (aodevent){
 270     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 271     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 272     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 273   }
 274
 275   header->SetMagneticField(event->GetMagneticField());
 276   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 277
 278   header->SetZDCN1Energy(event->GetZDCN1Energy());
 279   header->SetZDCP1Energy(event->GetZDCP1Energy());
 280   header->SetZDCN2Energy(event->GetZDCN2Energy());
 281   header->SetZDCP2Energy(event->GetZDCP2Energy());
 282   header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
 283
 284   Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
 285   Float_t diamcov[3];
 286   event->GetDiamondCovXY(diamcov);
 287   header->SetDiamond(diamxy,diamcov);
 288   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 289   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 290
 291   //
 292   //
 293   Int_t nVertices = 1 ;/* = prim. vtx*/;
 294   Int_t nCaloClus = event->GetNumberOfCaloClusters();
 295
 296   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 297
 298   // Access to the AOD container of vertices
 299   TClonesArray &vertices = *(AODEvent()->GetVertices());
 300   Int_t jVertices=0;
 301
 302   // Add primary vertex. The primary tracks will be defined
 303   // after the loops on the composite objects (V0, cascades, kinks)
 304   event->GetPrimaryVertex()->GetXYZ(pos);
 305   Float_t chi = 0;
 306   if      (esdevent){
 307     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 308     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 309   }
 310   else if (aodevent){
 311     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 312     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 313   }
 314
 315   AliAODVertex * primary = new(vertices[jVertices++])
 316   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 317   primary->SetName(event->GetPrimaryVertex()->GetName());
 318   primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
 319
 320 }
 321
 322 //________________________________________________________________________
 323 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
 324 {
 325   // Main loop
 326   // Called for each event
 327
 328   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 329   Int_t eventN = Entry();
 330   if(aodIH) eventN = aodIH->GetReadEntry();
 331
 332   if (eventN > fMaxEvent) return;
 333   //printf("Clusterizer --- Event %d-- \n",eventN);
 334
 335   //Remove the contents of output list set in the previous event
 336   fOutputAODBranch->Clear("C");
 337
 338   //Magic line to write events to AOD file
 339   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
 340   LoadBranches();
 341
 342   //Init pointers, clusterizer, ocdb
 343   if(fAccessOCDB) AccessOCDB();
 344   InitClusterization();
 345
 346   //Get the event
 347   AliVEvent   * event    = 0;
 348   AliESDEvent * esdevent = 0;
 349
 350   //Fill output event with header
 351
 352   //Check if input event are embedded events
 353   //If so, take output event
 354   if (aodIH && aodIH->GetMergeEvents()) {
 355     //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
 356     event  = AODEvent();
 357
 358     if(!aodIH->GetMergeEMCALCells())
 359       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
 360
 361     //    printf("InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
 362     //           InputEvent()->GetEMCALCells()->GetNumberOfCells());
 363     //    printf("MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
 364     //           aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
 365     //    for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
 366     //        AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
 367     //        if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
 368     //    }
 369     //    printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
 370     //           AODEvent()->GetEMCALCells()->GetNumberOfCells());
 371
 372   }
 373   else {
 374     event =  InputEvent();
 375     esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 376     if(fFillAODCaloCells) FillAODCaloCells();
 377     if(fFillAODHeader)    FillAODHeader();
 378   }
 379
 380   if (!event) {
 381     Error("UserExec","Event not available");
 382     return;
 383   }
 384
 385   //-------------------------------------------------------------------------------------
 386   //Set the geometry matrix, for the first event, skip the rest
 387   //-------------------------------------------------------------------------------------
 388   if(!fGeomMatrixSet){
 389     if(fLoadGeomMatrices){
 390       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 391         if(fGeomMatrix[mod]){
 392           if(DebugLevel() > 1)
 393             fGeomMatrix[mod]->Print();
 394           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
 395         }
 396         fGeomMatrixSet=kTRUE;
 397       }//SM loop
 398     }//Load matrices
 399     else if(!gGeoManager){
 400       Info("UserExec","Get geo matrices from data");
 401       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 402       if(!strcmp(event->GetName(),"AliAODEvent")) {
 403         if(DebugLevel() > 1)
 404           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 405       }//AOD
 406       else {
 407         if(DebugLevel() > 1)
 408           Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 409         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 410         if(!esd) {
 411           Error("UserExec","This event does not contain ESDs?");
 412           return;
 413         }
 414         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 415           if(DebugLevel() > 1)
 416             esd->GetEMCALMatrix(mod)->Print();
 417           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 418         }
 419         fGeomMatrixSet=kTRUE;
 420       }//ESD
 421     }//Load matrices from Data
 422
 423     //Recover time dependent corrections, put then in recalibration histograms. Do it once
 424     fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
 425
 426   }//first event
 427
 428   //Get the list of cells needed for unfolding and reclustering
 429   AliVCaloCells *cells   = event->GetEMCALCells();
 430   Int_t    nClustersOrg  = 0;
 431   Double_t cellAmplitude = 0;
 432   Double_t cellTime      = 0;
 433   Short_t  cellNumber    = 0;
 434
 435   //-------------------------------------------
 436   //---------Unfolding clusters----------------
 437   //-------------------------------------------
 438   if (fJustUnfold) {
 439
 440     //Fill the array with the EMCAL clusters, copy them
 441     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
 442     {
 443       AliVCluster *clus = event->GetCaloCluster(i);
 444       if(clus->IsEMCAL()){
 445
 446         //recalibrate/remove bad channels/etc if requested
 447         if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
 448           continue;
 449         }
 450
 451         if(fRecoUtils->IsRecalibrationOn()){
 452
 453           //Calibrate cluster
 454           //printf("Energy before %f ",clus->E());
 455           fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 456           //printf("; after %f\n",clus->E());
 457
 458           //CalibrateCells
 459           for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 460           {
 461             if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 462               break;
 463
 464
 465             Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 466             fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 467             fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 468
 469             //Do not include bad channels found in analysis?
 470             if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 471                 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 472               fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
 473               fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
 474               continue;
 475             }
 476
 477             cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 478
 479           }// cells loop
 480         }// recalibrate
 481
 482         //Cast to ESD or AOD, needed to create the cluster array
 483         AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 484         AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 485         if     (esdCluster){
 486           fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 487         }//ESD
 488         else if(aodCluster){
 489           fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 490         }//AOD
 491         else
 492           Warning("UserExec()"," - Wrong CaloCluster type?");
 493         nClustersOrg++;
 494       }
 495     }
 496
 497     //Do the unfolding
 498     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 499
 500     //CLEAN-UP
 501     fUnfolder->Clear();
 502
 503   } // just unfold ESD/AOD cluster
 504
 505   //-------------------------------------------
 506   //---------- Recluster cells ----------------
 507   //-------------------------------------------
 508
 509   else{
 510     //-------------------------------------------------------------------------------------
 511     //Transform CaloCells into Digits
 512     //-------------------------------------------------------------------------------------
 513
 514     //In case of MC, first loop on the clusters and fill MC label to array
 515     //.....................................................................
 516
 517     //    for(Int_t j = 0; j < 24*48*11; j++)  {
 518     //      if(fCellLabels[j]      !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j]      ) ;
 519     //      if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
 520     //    }
 521
 522     Int_t nClusters = event->GetNumberOfCaloClusters();
 523     if(aodIH && aodIH->GetEventToMerge())  //Embedding
 524       nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
 525     for (Int_t i = 0; i < nClusters; i++)
 526     {
 527       AliVCluster *clus = 0;
 528       if(aodIH && aodIH->GetEventToMerge()) //Embedding
 529         clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
 530       else
 531         clus = event->GetCaloCluster(i);
 532
 533       if(!clus) {
 534         printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
 535         return;
 536       }
 537
 538       if(clus->IsEMCAL()){
 539         //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
 540         Int_t label = clus->GetLabel();
 541         Int_t label2 = -1 ;
 542         if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 543         UShort_t * index    = clus->GetCellsAbsId() ;
 544         for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
 545           fCellLabels[index[icell]]      =label;
 546           fCellSecondLabels[index[icell]]=label2;
 547           //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
 548         }
 549         nClustersOrg++;
 550       }
 551     }
 552
 553     // Create digits
 554     //.................
 555     Int_t idigit =  0;
 556     Int_t id     = -1;
 557     Float_t amp  = -1;
 558     Float_t time = -1;
 559
 560     TTree *digitsTree = new TTree("digitstree","digitstree");
 561     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 562
 563     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 564     {
 565       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 566         break;
 567
 568       time = cellTime;
 569       amp  = cellAmplitude;
 570       id   = cellNumber;
 571
 572       Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 573       fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
 574       fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 575
 576       //Do not include bad channels found in analysis?
 577       if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 578          fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 579         fCellLabels[id]      =-1; //reset the entry in the array for next event
 580         fCellSecondLabels[id]=-1; //reset the entry in the array for next event
 581         //printf("Remove channel %d\n",id);
 582         continue;
 583       }
 584
 585       //Recalibrate?
 586       if(fRecoUtils->IsRecalibrationOn()){
 587         //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
 588         amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 589       }
 590
 591       //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 592       new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 593       //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
 594       //else                  printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
 595       fCellLabels[id]      =-1; //reset the entry in the array for next event
 596
 597       //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
 598       idigit++;
 599     }
 600
 601     //Fill the tree with digits
 602     digitsTree->Fill();
 603
 604     //-------------------------------------------------------------------------------------
 605     //Do the clusterization
 606     //-------------------------------------------------------------------------------------
 607     TTree *clustersTree = new TTree("clustertree","clustertree");
 608
 609     fClusterizer->SetInput(digitsTree);
 610     fClusterizer->SetOutput(clustersTree);
 611     fClusterizer->Digits2Clusters("");
 612
 613     //-------------------------------------------------------------------------------------
 614     //Transform the recpoints into AliVClusters
 615     //-------------------------------------------------------------------------------------
 616
 617     clustersTree->SetBranchStatus("*",0); //disable all branches
 618     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 619
 620     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 621     branch->SetAddress(&fClusterArr);
 622     branch->GetEntry(0);
 623
 624     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 625
 626     if(!fCaloClusterArr){
 627       printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 628       return;
 629     }
 630
 631     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
 632       printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters: Input entries %d - Output entries %d - %d (not fast)\n",
 633              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 634     }
 635
 636     //Reset the array with second labels for this event
 637     memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
 638
 639     //---CLEAN UP-----
 640     fClusterizer->Clear();
 641     fDigitsArr  ->Clear("C");
 642     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 643
 644     clustersTree->Delete("all");
 645     digitsTree  ->Delete("all");
 646   }
 647
 648   //Recalculate track-matching for the new clusters, only with ESDs
 649   if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
 650
 651
 652   //-------------------------------------------------------------------------------------
 653   //Put the new clusters in the AOD list
 654   //-------------------------------------------------------------------------------------
 655
 656   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
 657   //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
 658   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
 659     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 660     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
 661
 662     //Add matched track, if any, only with ESDs
 663     if(esdevent && fDoTrackMatching){
 664       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 665       if(trackIndex >= 0){
 666         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
 667         if(DebugLevel() > 1)
 668           Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
 669       }
 670     }
 671
 672     //In case of new bad channels, recalculate distance to bad channels
 673     if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
 674       //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
 675       fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
 676       //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
 677     }
 678
 679     //    if(newCluster->GetNLabels()>0){
 680     //      printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
 681     //      UShort_t * newindex    = newCluster->GetCellsAbsId() ;
 682     //      for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
 683     //       printf("\t absID %d\n",newindex[icell]);
 684     //     }
 685     //    }
 686
 687     //    printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
 688     //    printf("labels: ");
 689     //    for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
 690     //      printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
 691     //     printf("\n ");
 692
 693     newCluster->SetID(i);
 694     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 695   }
 696
 697   //if(fOutputAODBranchName.Length()!=0)
 698   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
 699
 700   //---CLEAN UP-----
 701   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
 702 }
 703
 704 //_____________________________________________________________________
 705 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 706 {
 707   //Access to OCDB stuff
 708
 709   AliVEvent * event = InputEvent();
 710   if (!event)
 711   {
 712     Warning("AccessOCDB","Event not available!!!");
 713     return kFALSE;
 714   }
 715
 716   if (event->GetRunNumber()==fRun)
 717     return kTRUE;
 718   fRun = event->GetRunNumber();
 719
 720   if(DebugLevel() > 1 )
 721     Info("AccessODCD()"," Begin");
 722
 723   //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 724
 725   AliCDBManager *cdb = AliCDBManager::Instance();
 726
 727
 728   if (fOCDBpath.Length()){
 729     cdb->SetDefaultStorage(fOCDBpath.Data());
 730     Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
 731   }
 732
 733   cdb->SetRun(event->GetRunNumber());
 734
 735   //
 736   // EMCAL from RAW OCDB
 737   if (fOCDBpath.Contains("alien:"))
 738   {
 739     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 740     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 741   }
 742
 743   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 744
 745   // init parameters:
 746
 747   //Get calibration parameters
 748   if(!fCalibData)
 749   {
 750     AliCDBEntry *entry = (AliCDBEntry*)
 751     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 752
 753     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 754   }
 755
 756   if(!fCalibData)
 757     AliFatal("Calibration parameters not found in CDB!");
 758
 759   //Get calibration parameters
 760   if(!fPedestalData)
 761   {
 762     AliCDBEntry *entry = (AliCDBEntry*)
 763     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 764
 765     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 766   }
 767
 768   if(!fPedestalData)
 769     AliFatal("Dead map not found in CDB!");
 770
 771   return kTRUE;
 772 }
 773
 774 //________________________________________________________________________________________
 775 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 776 {
 777   //Select clusterization/unfolding algorithm and set all the needed parameters
 778
 779   if (fJustUnfold){
 780     // init the unfolding afterburner
 781     delete fUnfolder;
 782     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 783     return;
 784   }
 785
 786   //First init the clusterizer
 787   delete fClusterizer;
 788   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 789     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 790   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 791     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
 792   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 793     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 794   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerv2) { //FIX this other way.
 795     AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 796     clusterizer->SetNRowDiff(2);
 797     clusterizer->SetNColDiff(2);
 798     fClusterizer = clusterizer;
 799   } else{
 800     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 801   }
 802
 803   //Now set the parameters
 804   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 805   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 806   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 807   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 808   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 809   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 810   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 811   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 812   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 813   fClusterizer->SetJustClusters          ( kTRUE                               );
 814   //In case of unfolding after clusterization is requested, set the corresponding parameters
 815   if(fRecParam->GetUnfold()){
 816     Int_t i=0;
 817     for (i = 0; i < 8; i++) {
 818       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 819     }//end of loop over parameters
 820     for (i = 0; i < 3; i++) {
 821       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 822       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 823     }//end of loop over parameters
 824
 825     fClusterizer->InitClusterUnfolding();
 826   }// to unfold
 827 }
 828
 829 //________________________________________________________________________________________
 830 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
 831 {
 832   // Restore clusters from recPoints
 833   // Cluster energy, global position, cells and their amplitude fractions are restored
 834   Int_t j = 0;
 835   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
 836   {
 837     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
 838
 839     const Int_t ncells = recPoint->GetMultiplicity();
 840     Int_t ncells_true = 0;
 841
 842     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
 843
 844     // cells and their amplitude fractions
 845     UShort_t   absIds[ncells];
 846     Double32_t ratios[ncells];
 847
 848     for (Int_t c = 0; c < ncells; c++) {
 849       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
 850
 851       absIds[ncells_true] = digit->GetId();
 852       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 853
 854       if (ratios[ncells_true] > 0.001) ncells_true++;
 855     }
 856
 857     if (ncells_true < 1) {
 858       printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",recPoint->GetEnergy(), ncells);
 859       continue;
 860     }
 861
 862     TVector3 gpos;
 863     Float_t g[3];
 864
 865     // calculate new cluster position
 866     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
 867     recPoint->GetGlobalPosition(gpos);
 868     gpos.GetXYZ(g);
 869
 870     // create a new cluster
 871     //AliAODCaloCluster *clus = new AliAODCaloCluster();
 872     (*clusArray)[j] = new AliAODCaloCluster() ;
 873     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
 874     j++;
 875     clus->SetType(AliVCluster::kEMCALClusterv1);
 876     clus->SetE(recPoint->GetEnergy());
 877     clus->SetPosition(g);
 878     clus->SetNCells(ncells_true);
 879     clus->SetCellsAbsId(absIds);
 880     clus->SetCellsAmplitudeFraction(ratios);
 881     clus->SetDispersion(recPoint->GetDispersion());
 882     clus->SetChi2(-1); //not yet implemented
 883     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
 884     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
 885     Float_t elipAxis[2];
 886     recPoint->GetElipsAxis(elipAxis);
 887     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 888     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 889     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
 890
 891     //MC
 892     Int_t  parentMult  = 0;
 893     Int_t *parentList = recPoint->GetParents(parentMult);
 894     clus->SetLabel(parentList, parentMult);
 895
 896     //    if(parentMult)printf("RecToESD: Labels %d",parentMult);
 897     //    for (Int_t iii = 0; iii < parentMult; iii++) {
 898     //      printf("\t label %d\n",parentList[iii]);
 899     //    }
 900
 901     //Write the second major contributor to each MC cluster.
 902     Int_t iNewLabel ;
 903     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
 904
 905       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
 906       if(idCell>=0){
 907         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
 908         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
 909         {
 910           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
 911           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
 912         }
 913         if (iNewLabel == 1)
 914         {
 915           //Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
 916           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
 917           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
 918             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
 919           }
 920
 921           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
 922           //if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
 923           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
 924           delete [] newLabelArray;
 925         }
 926       }//positive cell number
 927     }
 928
 929   } // recPoints loop
 930
 931 }