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