PWGGA/EMCALTasks/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
  16 //_________________________________________________________________________
  17 // This analysis provides a new list of clusters to be used in other analysis
  18 //
  19 // Author: Gustavo Conesa Balbastre,
  20 //         Adapted from analysis class from Deepa Thomas
  21 //
  22 // $Id$
  23 //_________________________________________________________________________
  24
  25 // --- Root ---
  26 #include <TString.h>
  27 #include <TRefArray.h>
  28 #include <TClonesArray.h>
  29 #include <TTree.h>
  30 #include <TGeoManager.h>
  31 #include <TROOT.h>
  32 #include <TInterpreter.h>
  33 #include <TFile.h>
  34
  35 // --- AliRoot Analysis Steering
  36 #include "AliAnalysisTask.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliESDEvent.h"
  39 #include "AliGeomManager.h"
  40 #include "AliVCaloCells.h"
  41 #include "AliAODCaloCluster.h"
  42 #include "AliCDBManager.h"
  43 #include "AliCDBStorage.h"
  44 #include "AliCDBEntry.h"
  45 #include "AliLog.h"
  46 #include "AliVEventHandler.h"
  47 #include "AliAODInputHandler.h"
  48 #include "AliOADBContainer.h"
  49
  50 // --- EMCAL
  51 #include "AliEMCALAfterBurnerUF.h"
  52 #include "AliEMCALGeometry.h"
  53 #include "AliEMCALClusterizerNxN.h"
  54 #include "AliEMCALClusterizerv1.h"
  55 #include "AliEMCALClusterizerv2.h"
  56 #include "AliEMCALRecPoint.h"
  57 #include "AliEMCALDigit.h"
  58 #include "AliCaloCalibPedestal.h"
  59 #include "AliEMCALCalibData.h"
  60
  61 #include "AliAnalysisTaskEMCALClusterize.h"
  62
  63 ClassImp(AliAnalysisTaskEMCALClusterize)
  64
  65 //______________________________________________________________________________
  66 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  67 : AliAnalysisTaskSE(name)
  68 , fEvent(0)
  69 , fGeom(0),               fGeomName("")
  70 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  71 , fCalibData(0),          fPedestalData(0)
  72 , fOCDBpath(""),          fAccessOCDB(kFALSE)
  73 , fDigitsArr(0),          fClusterArr(0),             fCaloClusterArr(0)
  74 , fRecParam(0),           fClusterizer(0)
  75 , fUnfolder(0),           fJustUnfold(kFALSE)
  76 , fOutputAODBranch(0),    fOutputAODBranchName(""),   fOutputAODBranchSet(0)
  77 , fFillAODFile(kFALSE),   fFillAODHeader(0)
  78 , fFillAODCaloCells(0),   fRun(-1)
  79 , fRecoUtils(0),          fConfigName("")
  80 , fCellLabels(),          fCellSecondLabels(),        fCellTime()
  81 , fCellMatchdEta(),       fCellMatchdPhi()
  82 , fMaxEvent(0),           fDoTrackMatching(kFALSE)
  83 , fSelectCell(kFALSE),    fSelectCellMinE(0),         fSelectCellMinFrac(0)
  84 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
  85 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
  86 , fOADBSet(kFALSE),       fAccessOADB(kTRUE),         fOADBFilePath("")
  87 , fCentralityClass("")
  88 {
  89   // Constructor
  90
  91   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
  92   for(Int_t j = 0; j < 24*48*11; j++)
  93   {
  94     fCellLabels[j]       = -1;
  95     fCellSecondLabels[j] = -1;
  96     fCellTime[j]         =  0.;
  97     fCellMatchdEta[j]    = -999;
  98     fCellMatchdPhi[j]    = -999;
  99   }
 100
 101   fCentralityBin[0] = fCentralityBin[1]=-1;
 102
 103 }
 104
 105 //______________________________________________________________
 106 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
 107 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
 108 , fEvent(0)
 109 , fGeom(0),                 fGeomName("")
 110 , fGeomMatrixSet(kFALSE),   fLoadGeomMatrices(kFALSE)
 111 , fCalibData(0),            fPedestalData(0)
 112 , fOCDBpath(""),            fAccessOCDB(kFALSE)
 113 , fDigitsArr(0),            fClusterArr(0),             fCaloClusterArr(0)
 114 , fRecParam(0),             fClusterizer(0)
 115 , fUnfolder(0),             fJustUnfold(kFALSE)
 116 , fOutputAODBranch(0),      fOutputAODBranchName(""),   fOutputAODBranchSet(0)
 117 , fFillAODFile(kFALSE),     fFillAODHeader(0)
 118 , fFillAODCaloCells(0),     fRun(-1)
 119 , fRecoUtils(0),            fConfigName("")
 120 , fCellLabels(),            fCellSecondLabels(),        fCellTime()
 121 , fCellMatchdEta(),         fCellMatchdPhi()
 122 , fMaxEvent(0),             fDoTrackMatching(kFALSE)
 123 , fSelectCell(kFALSE),      fSelectCellMinE(0),         fSelectCellMinFrac(0)
 124 , fRemoveLEDEvents(kTRUE),  fRemoveExoticEvents(kFALSE)
 125 , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
 126 , fOADBSet(kFALSE),         fAccessOADB(kTRUE),        fOADBFilePath("")
 127 , fCentralityClass("")
 128 {
 129   // Constructor
 130
 131   for(Int_t i = 0; i < 12;    i++)  fGeomMatrix[i] =  0;
 132   for(Int_t j = 0; j < 24*48*11; j++)
 133   {
 134     fCellLabels[j]       = -1;
 135     fCellSecondLabels[j] = -1;
 136     fCellTime[j]         =  0.;
 137     fCellMatchdEta[j]    = -999;
 138     fCellMatchdPhi[j]    = -999;
 139   }
 140
 141   fCentralityBin[0] = fCentralityBin[1]=-1;
 142
 143 }
 144
 145
 146 //_______________________________________________________________
 147 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 148 {
 149   //dtor
 150
 151   if (fDigitsArr)
 152   {
 153     fDigitsArr->Clear("C");
 154     delete fDigitsArr;
 155   }
 156
 157   if (fClusterArr)
 158   {
 159     fClusterArr->Delete();
 160     delete fClusterArr;
 161   }
 162
 163   if (fCaloClusterArr)
 164   {
 165     fCaloClusterArr->Delete();
 166     delete fCaloClusterArr;
 167   }
 168
 169   if(fClusterizer) delete fClusterizer;
 170   if(fUnfolder)    delete fUnfolder;
 171   if(fRecoUtils)   delete fRecoUtils;
 172
 173 }
 174
 175 //_______________________________________________
 176 void AliAnalysisTaskEMCALClusterize::AccessOADB()
 177 {
 178   // Set the AODB calibration, bad channels etc. parameters at least once
 179   // alignment matrices from OADB done in SetGeometryMatrices
 180
 181   //Set it only once
 182   if(fOADBSet) return ;
 183
 184   Int_t   runnumber = InputEvent()->GetRunNumber() ;
 185   TString pass      = GetPass();
 186
 187   printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
 188
 189   Int_t nSM = fGeom->GetNumberOfSuperModules();
 190
 191   // Bad map
 192   if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
 193   {
 194     AliOADBContainer *contBC=new AliOADBContainer("");
 195     contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
 196
 197     TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
 198
 199     if(arrayBC)
 200     {
 201         fRecoUtils->SwitchOnDistToBadChannelRecalculation();
 202         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
 203
 204         for (Int_t i=0; i<nSM; ++i)
 205         {
 206           TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
 207
 208           if (hbm)
 209             delete hbm;
 210
 211           hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
 212
 213           if (!hbm)
 214           {
 215             AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
 216             continue;
 217           }
 218
 219           hbm->SetDirectory(0);
 220           fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
 221
 222         } // loop
 223     } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
 224   }  // Remove bad
 225
 226   // Energy Recalibration
 227   if(fRecoUtils->IsRecalibrationOn())
 228   {
 229     AliOADBContainer *contRF=new AliOADBContainer("");
 230
 231     contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
 232
 233     TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
 234
 235     if(recal)
 236     {
 237       TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
 238
 239       if(recalpass)
 240       {
 241         TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
 242
 243         if(recalib)
 244         {
 245           printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
 246           for (Int_t i=0; i<nSM; ++i)
 247           {
 248             TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
 249
 250             if (h)
 251               delete h;
 252
 253             h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
 254
 255             if (!h)
 256             {
 257               AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
 258               continue;
 259             }
 260
 261             h->SetDirectory(0);
 262
 263             fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
 264           } // SM loop
 265         }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok
 266       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
 267     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n");  // run number array ok
 268
 269   } // Recalibration on
 270
 271   // Energy Recalibration, apply on top of previous calibration factors
 272   if(fRecoUtils->IsRunDepRecalibrationOn())
 273   {
 274     AliOADBContainer *contRFTD=new AliOADBContainer("");
 275
 276     contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
 277
 278     TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
 279
 280     if(htd)
 281     {
 282       printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
 283
 284       for (Int_t ism=0; ism<nSM; ++ism)
 285       {
 286         for (Int_t icol=0; icol<48; ++icol)
 287         {
 288           for (Int_t irow=0; irow<24; ++irow)
 289           {
 290             Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
 291
 292             Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
 293             factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
 294             //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
 295             //      GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
 296             fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
 297           } // columns
 298         } // rows
 299       } // SM loop
 300     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
 301   } // Run by Run T calibration
 302
 303   // Time Recalibration
 304   if(fRecoUtils->IsTimeRecalibrationOn())
 305   {
 306     AliOADBContainer *contTRF=new AliOADBContainer("");
 307
 308     contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
 309
 310     TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
 311
 312     if(trecal)
 313     {
 314       TString passTmp = pass;
 315       if(pass!="pass1" && pass!="pass2") passTmp = "pass2"; // TEMPORARY FIX FOR LHC11a analysis
 316       TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passTmp);
 317
 318       if(trecalpass)
 319       {
 320         printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n");
 321         for (Int_t ibc = 0; ibc < 4; ++ibc)
 322         {
 323           TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
 324
 325           if (h)
 326             delete h;
 327
 328           h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
 329
 330           if (!h)
 331           {
 332             AliError(Form("Could not load hAllTimeAvBC%d",ibc));
 333             continue;
 334           }
 335
 336           h->SetDirectory(0);
 337
 338           fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
 339         } // bunch crossing loop
 340       }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
 341     }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n");  // run number array ok
 342
 343   } // Time recalibration on
 344
 345   // Parameters already set once, so do not it again
 346   fOADBSet = kTRUE;
 347
 348 }
 349
 350 //_________________________________________________
 351 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 352 {
 353   //Access to OCDB stuff
 354
 355   fEvent = InputEvent();
 356   if (!fEvent)
 357   {
 358     Warning("AccessOCDB","Event not available!!!");
 359     return kFALSE;
 360   }
 361
 362   if (fEvent->GetRunNumber()==fRun)
 363     return kTRUE;
 364   fRun = fEvent->GetRunNumber();
 365
 366   if(DebugLevel() > 1 )
 367     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Begin");
 368
 369   AliCDBManager *cdb = AliCDBManager::Instance();
 370
 371
 372   if (fOCDBpath.Length())
 373   {
 374     cdb->SetDefaultStorage(fOCDBpath.Data());
 375     printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
 376   }
 377
 378   cdb->SetRun(fEvent->GetRunNumber());
 379
 380   //
 381   // EMCAL from RAW OCDB
 382   if (fOCDBpath.Contains("alien:"))
 383   {
 384     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 385     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 386   }
 387
 388   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 389
 390   // init parameters:
 391
 392   //Get calibration parameters
 393   if(!fCalibData)
 394   {
 395     AliCDBEntry *entry = (AliCDBEntry*)
 396     AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 397
 398     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 399   }
 400
 401   if(!fCalibData)
 402     AliFatal("Calibration parameters not found in CDB!");
 403
 404   //Get calibration parameters
 405   if(!fPedestalData)
 406   {
 407     AliCDBEntry *entry = (AliCDBEntry*)
 408     AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 409
 410     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 411   }
 412
 413   if(!fPedestalData)
 414     AliFatal("Dead map not found in CDB!");
 415
 416   return kTRUE;
 417 }
 418
 419 //_____________________________________________________
 420 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
 421 {
 422   // Get the input event, it can depend in embedded events what you want to get
 423   // Also check if the quality of the event is good if not reject it
 424
 425   fEvent = 0x0;
 426
 427   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 428   Int_t eventN = Entry();
 429   if(aodIH) eventN = aodIH->GetReadEntry();
 430
 431   if (eventN > fMaxEvent)
 432     return ;
 433
 434   //printf("Clusterizer --- Event %d-- \n",eventN);
 435
 436   //Check if input event are embedded events
 437   //If so, take output event
 438   if (aodIH && aodIH->GetMergeEvents())
 439   {
 440     fEvent  = AODEvent();
 441
 442     if(!aodIH->GetMergeEMCALCells())
 443       AliFatal("Events merged but not EMCAL cells, check analysis settings!");
 444
 445     if(DebugLevel() > 1)
 446     {
 447       printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
 448
 449       printf("\t InputEvent  N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
 450              InputEvent()->GetEMCALCells()->GetNumberOfCells());
 451
 452       printf("\t MergedEvent  N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
 453              aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
 454
 455       for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
 456       {
 457         AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
 458         if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E  %f\n",icl,sigCluster->E());
 459       }
 460
 461       printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
 462              AODEvent()->GetEMCALCells()->GetNumberOfCells());
 463     }
 464   }
 465   else
 466   {
 467     fEvent =  InputEvent();
 468     if(fFillAODCaloCells) FillAODCaloCells();
 469     if(fFillAODHeader)    FillAODHeader();
 470   }
 471
 472   if (!fEvent)
 473   {
 474     Error("UserExec","Event not available");
 475     return ;
 476   }
 477
 478   //-------------------------------------------------------------------------------------
 479   // Reject events if LED was firing, use only for LHC11a data
 480   // Reject event if triggered by exotic cell and remove exotic cells if not triggered
 481   //-------------------------------------------------------------------------------------
 482
 483   if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
 484
 485   if( IsExoticEvent() )                            { fEvent = 0x0 ; return ; }
 486
 487   //-------------------------------------------------------------------------------------
 488   // Set the cluster array in the event (output or input)
 489   //-------------------------------------------------------------------------------------
 490
 491   if     ( fFillAODFile )
 492   {
 493     //Magic line to write events to AOD filem put after event rejection
 494     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
 495   }
 496   else if( !fOutputAODBranchSet )
 497   {
 498     // Create array and put it in the input event, if output AOD not selected, only once
 499     InputEvent()->AddObject(fOutputAODBranch);
 500     fOutputAODBranchSet = kTRUE;
 501     printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data());
 502   }
 503
 504 }
 505
 506 //____________________________________________________
 507 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
 508 {
 509   // Recluster calocells, transform them into digits,
 510   // feed the clusterizer with them and get new list of clusters
 511
 512   //In case of MC, first loop on the clusters and fill MC label to array
 513   Int_t nClusters     = fEvent->GetNumberOfCaloClusters();
 514   Int_t nClustersOrg  = 0;
 515
 516   AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 517   if(aodIH && aodIH->GetEventToMerge())  //Embedding
 518     nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
 519
 520   for (Int_t i = 0; i < nClusters; i++)
 521   {
 522     AliVCluster *clus = 0;
 523     if(aodIH && aodIH->GetEventToMerge()) //Embedding
 524       clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
 525     else
 526       clus = fEvent->GetCaloCluster(i);
 527
 528     if(!clus) return;
 529
 530     if(clus->IsEMCAL())
 531     {
 532       Int_t label = clus->GetLabel();
 533       Int_t label2 = -1 ;
 534       //printf("Org cluster E %f, Time  %e, Id = ", clus->E(), clus->GetTOF() );
 535       if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
 536       UShort_t * index    = clus->GetCellsAbsId() ;
 537       for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
 538       {
 539         fCellLabels[index[icell]]       = label;
 540         fCellSecondLabels[index[icell]] = label2;
 541         fCellTime[index[icell]]         = clus->GetTOF();
 542         fCellMatchdEta[index[icell]]    = clus->GetTrackDz();
 543         fCellMatchdPhi[index[icell]]    = clus->GetTrackDx();
 544         //printf(" %d,", index[icell] );
 545       }
 546       nClustersOrg++;
 547     }
 548     // printf("\n");
 549   }
 550
 551   // Transform CaloCells into Digits
 552
 553   Int_t    idigit =  0;
 554   Int_t    id     = -1;
 555   Float_t  amp    = -1;
 556   Double_t time   = -1;
 557
 558   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 559
 560   Int_t bc = InputEvent()->GetBunchCrossNumber();
 561   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 562   {
 563     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
 564     id = cells->GetCellNumber(icell);
 565     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
 566
 567     // Do not include cells with too low energy, nor exotic cell
 568     if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
 569
 570     // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
 571     if (time*1e9 < 1.)
 572     {
 573       time = fCellTime[id];
 574       //printf("cell %d time org %f - ",id, time*1.e9);
 575       fRecoUtils->RecalibrateCellTime(id,bc,time);
 576       //printf("recal %f\n",time*1.e9);
 577     }
 578
 579     if(  accept && fRecoUtils->IsExoticCell(id,cells,bc))
 580     {
 581       accept = kFALSE;
 582     }
 583
 584     if( !accept )
 585     {
 586       fCellLabels[id]      =-1; //reset the entry in the array for next event
 587       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
 588       fCellTime[id]        = 0.;
 589       fCellMatchdEta[id]   =-999;
 590       fCellMatchdPhi[id]   =-999;
 591       if( DebugLevel() > 2 )
 592         printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
 593                id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
 594       continue;
 595     }
 596
 597     //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 598     new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 599
 600     fCellLabels[id] =-1; //reset the entry in the array for next event
 601
 602     idigit++;
 603   }
 604
 605   fDigitsArr->Sort();
 606
 607   //-------------------------------------------------------------------------------------
 608   //Do the clusterization
 609   //-------------------------------------------------------------------------------------
 610
 611   fClusterizer->Digits2Clusters("");
 612
 613   //-------------------------------------------------------------------------------------
 614   //Transform the recpoints into AliVClusters
 615   //-------------------------------------------------------------------------------------
 616
 617   RecPoints2Clusters();
 618
 619   if(!fCaloClusterArr)
 620   {
 621     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 622     return;
 623   }
 624
 625   if( DebugLevel() > 0 )
 626   {
 627     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
 628
 629     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
 630     {
 631       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
 632              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
 633              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 634     }
 635   }
 636
 637   //Reset the array with second labels for this event
 638   memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
 639
 640 }
 641
 642 //_____________________________________________________
 643 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
 644 {
 645   // Take the event clusters and unfold them
 646
 647   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 648   Double_t cellAmplitude = 0;
 649   Double_t cellTime      = 0;
 650   Short_t  cellNumber    = 0;
 651   Short_t  cellMCLabel   = 0;
 652   Double_t cellEFrac     = 0;
 653   Int_t    nClustersOrg  = 0;
 654
 655   // Fill the array with the EMCAL clusters, copy them
 656   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
 657   {
 658     AliVCluster *clus = fEvent->GetCaloCluster(i);
 659     if(clus->IsEMCAL())
 660     {
 661       //recalibrate/remove bad channels/etc if requested
 662       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
 663       {
 664         continue;
 665       }
 666
 667       if(fRecoUtils->IsRecalibrationOn())
 668       {
 669         //Calibrate cluster
 670         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 671
 672         //CalibrateCells
 673         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 674         {
 675           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
 676             break;
 677
 678           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 679           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 680           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 681
 682           //Do not include bad channels found in analysis?
 683           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 684              fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
 685           {
 686             fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
 687             fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
 688             fCellTime[cellNumber]        = 0.;
 689             fCellMatchdEta[cellNumber]   =-999;
 690             fCellMatchdPhi[cellNumber]   =-999;
 691             continue;
 692           }
 693
 694           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 695
 696         }// cells loop
 697       }// recalibrate
 698
 699       //Cast to ESD or AOD, needed to create the cluster array
 700       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 701       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 702
 703       if     (esdCluster)
 704       {
 705         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 706       }//ESD
 707       else if(aodCluster)
 708       {
 709         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 710       }//AOD
 711       else
 712         Warning("UserExec()"," - Wrong CaloCluster type?");
 713
 714       nClustersOrg++;
 715     }
 716   }
 717
 718   //Do the unfolding
 719   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 720
 721   //CLEAN-UP
 722   fUnfolder->Clear();
 723
 724 }
 725
 726 //_____________________________________________________
 727 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
 728 {
 729   // Put calo cells in standard branch
 730   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
 731   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 732
 733   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 734   aodEMcells.CreateContainer(nEMcell);
 735   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 736   Double_t calibFactor = 1.;
 737   for (Int_t iCell = 0; iCell < nEMcell; iCell++)
 738   {
 739     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 740     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 741     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 742
 743     if(fRecoUtils->IsRecalibrationOn())
 744     {
 745       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 746     }
 747
 748     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
 749     { //Channel is not declared as bad
 750       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
 751                          eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
 752     }
 753     else
 754     {
 755       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
 756     }
 757   }
 758   aodEMcells.Sort();
 759
 760 }
 761
 762 //__________________________________________________
 763 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
 764 {
 765   //Put event header information in standard AOD branch
 766
 767   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
 768   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
 769
 770   Double_t pos[3]   ;
 771   Double_t covVtx[6];
 772   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 773
 774   AliAODHeader* header = AODEvent()->GetHeader();
 775   header->SetRunNumber(fEvent->GetRunNumber());
 776
 777   if(esdevent)
 778   {
 779     TTree* tree = fInputHandler->GetTree();
 780     if (tree)
 781     {
 782       TFile* file = tree->GetCurrentFile();
 783       if (file) header->SetESDFileName(file->GetName());
 784     }
 785   }
 786   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 787
 788   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
 789   header->SetOrbitNumber(fEvent->GetOrbitNumber());
 790   header->SetPeriodNumber(fEvent->GetPeriodNumber());
 791   header->SetEventType(fEvent->GetEventType());
 792
 793   //Centrality
 794   if(fEvent->GetCentrality())
 795     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
 796   else
 797     header->SetCentrality(0);
 798
 799   //Trigger
 800   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 801   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 802   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 803
 804   header->SetTriggerMask(fEvent->GetTriggerMask());
 805   header->SetTriggerCluster(fEvent->GetTriggerCluster());
 806
 807   if      (esdevent)
 808   {
 809     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 810     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 811     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 812   }
 813   else if (aodevent)
 814   {
 815     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 816     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 817     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 818   }
 819
 820   header->SetMagneticField(fEvent->GetMagneticField());
 821   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 822
 823   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 824   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 825   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 826   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 827   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
 828
 829   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
 830   Float_t diamcov[3];
 831   fEvent->GetDiamondCovXY(diamcov);
 832   header->SetDiamond(diamxy,diamcov);
 833   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 834   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 835
 836   //
 837   Int_t nVertices = 1 ;/* = prim. vtx*/;
 838   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
 839
 840   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 841
 842   // Access to the AOD container of vertices
 843   TClonesArray &vertices = *(AODEvent()->GetVertices());
 844   Int_t jVertices=0;
 845
 846   // Add primary vertex. The primary tracks will be defined
 847   // after the loops on the composite objects (V0, cascades, kinks)
 848   fEvent->GetPrimaryVertex()->GetXYZ(pos);
 849   Float_t chi = 0;
 850   if      (esdevent)
 851   {
 852     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 853     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 854   }
 855   else if (aodevent)
 856   {
 857     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 858     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 859   }
 860
 861   AliAODVertex * primary = new(vertices[jVertices++])
 862   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 863   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
 864   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
 865
 866 }
 867
 868 //___________________________________________________________
 869 void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent()
 870 {
 871   // Get the CaloClusters array, do some final calculations
 872   // and put the clusters in the output or input event
 873   // as a separate branch
 874
 875   //First recalculate track-matching for the new clusters
 876   if(fDoTrackMatching)
 877   {
 878     fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
 879   }
 880   //Put the new clusters in the AOD list
 881
 882   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
 883
 884   for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
 885   {
 886     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 887
 888     newCluster->SetID(i);
 889
 890     // Correct cluster energy non linearity
 891     newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
 892
 893     //Add matched track
 894     if(fDoTrackMatching)
 895     {
 896       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 897       if(trackIndex >= 0)
 898       {
 899         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
 900         if(DebugLevel() > 1)
 901           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
 902       }
 903
 904       Float_t dR = 999., dZ = 999.;
 905       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
 906       newCluster->SetTrackDistance(dR,dZ);
 907
 908     }
 909     else
 910     {// Assign previously assigned matched track in reco, very very rough
 911       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
 912       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
 913     }
 914
 915     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
 916     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
 917     //printf("\n");
 918
 919     // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
 920     fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
 921
 922     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 923
 924     if(DebugLevel() > 1 )
 925       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
 926
 927   } // cluster loop
 928
 929   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
 930
 931
 932 }
 933
 934 //_______________________________________________
 935 TString AliAnalysisTaskEMCALClusterize::GetPass()
 936 {
 937   // Get passx from filename.
 938
 939   if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
 940   {
 941     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
 942     return TString("");
 943   }
 944
 945   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
 946   {
 947     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
 948     return TString("");
 949   }
 950
 951   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 952   if      (pass.Contains("ass1")) return TString("pass1");
 953   else if (pass.Contains("ass2")) return TString("pass2");
 954   else if (pass.Contains("ass3")) return TString("pass3");
 955   else if (pass.Contains("ass4")) return TString("pass4");
 956   else if (pass.Contains("ass5")) return TString("pass5");
 957
 958   // No condition fullfilled, give a default value
 959   printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
 960   return TString("");
 961
 962 }
 963
 964 //_________________________________________
 965 void AliAnalysisTaskEMCALClusterize::Init()
 966 {
 967   //Init analysis with configuration macro if available
 968
 969   fOADBSet           = kFALSE;
 970   if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
 971
 972   fBranchNames       = "ESD:AliESDHeader.,EMCALCells.";
 973
 974   if(!fRecParam)     fRecParam  = new AliEMCALRecParam;
 975   if(!fRecoUtils)    fRecoUtils = new AliEMCALRecoUtils();
 976
 977   if(fMaxEvent          <= 0) fMaxEvent          = 1000000000;
 978   if(fSelectCellMinE    <= 0) fSelectCellMinE    = 0.005;
 979   if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
 980
 981   //Centrality
 982   if(fCentralityClass  == "") fCentralityClass  = "V0M";
 983   fCentralityBin[0] = fCentralityBin[1]=-1;
 984
 985   if (fOCDBpath            == "") fOCDBpath            = "raw://" ;
 986   if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
 987
 988   if(gROOT->LoadMacro(fConfigName) >=0)
 989   {
 990     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
 991     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 992     fGeomName         = clus->fGeomName;
 993     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 994     fOCDBpath         = clus->fOCDBpath;
 995     fAccessOCDB       = clus->fAccessOCDB;
 996     fRecParam         = clus->fRecParam;
 997     fJustUnfold       = clus->fJustUnfold;
 998     fFillAODFile      = clus->fFillAODFile;
 999     fRecoUtils        = clus->fRecoUtils;
1000     fConfigName       = clus->fConfigName;
1001     fMaxEvent         = clus->fMaxEvent;
1002     fDoTrackMatching  = clus->fDoTrackMatching;
1003     fOutputAODBranchName = clus->fOutputAODBranchName;
1004     for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1005     fCentralityClass  = clus->fCentralityClass;
1006     fCentralityBin[0] = clus->fCentralityBin[0];
1007     fCentralityBin[1] = clus->fCentralityBin[1];
1008   }
1009
1010   // Init geometry, I do not like much to do it like this ...
1011   if(fImportGeometryFromFile && !gGeoManager)
1012   {
1013     if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/OADB/EMCAL/geometry_2011.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1014     printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
1015     TGeoManager::Import(fImportGeometryFilePath) ;
1016   }
1017
1018 }
1019
1020 //_______________________________________________________
1021 void AliAnalysisTaskEMCALClusterize::InitClusterization()
1022 {
1023   //Select clusterization/unfolding algorithm and set all the needed parameters
1024
1025   if (fJustUnfold)
1026   {
1027     // init the unfolding afterburner
1028     delete fUnfolder;
1029     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1030     return;
1031   }
1032
1033   //First init the clusterizer
1034   delete fClusterizer;
1035   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1036     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
1037   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1038     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
1039   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1040   {
1041     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
1042     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1043     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1044   }
1045   else
1046   {
1047     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1048   }
1049
1050   //Now set the parameters
1051   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1052   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
1053   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
1054   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
1055   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
1056   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
1057   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
1058   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
1059   fClusterizer->SetInputCalibrated       ( kTRUE                               );
1060   fClusterizer->SetJustClusters          ( kTRUE                               );
1061
1062   // Initialize the cluster rec points and digits arrays and get them.
1063   fClusterizer->SetOutput(0);
1064   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1065   fDigitsArr  = fClusterizer->GetDigits();
1066
1067   //In case of unfolding after clusterization is requested, set the corresponding parameters
1068   if(fRecParam->GetUnfold())
1069   {
1070     Int_t i=0;
1071     for (i = 0; i < 8; i++)
1072     {
1073       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1074     }//end of loop over parameters
1075
1076     for (i = 0; i < 3; i++)
1077     {
1078       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
1079       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
1080     }//end of loop over parameters
1081
1082     fClusterizer->InitClusterUnfolding();
1083
1084   }// to unfold
1085 }
1086
1087 //_________________________________________________
1088 void AliAnalysisTaskEMCALClusterize::InitGeometry()
1089 {
1090   // Init geometry and set the geometry matrix, for the first event, skip the rest
1091   // Also set once the run dependent calibrations
1092
1093   if(fGeomMatrixSet) return;
1094
1095   Int_t runnumber = InputEvent()->GetRunNumber() ;
1096   if (!fGeom)
1097   {
1098     if(fGeomName=="")
1099     {
1100       if     (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
1101       else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
1102       else                        fGeomName = "EMCAL_COMPLETE12SMV1";
1103       printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
1104     }
1105
1106                 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1107
1108                 if(fDebug > 0)
1109     {
1110                         printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
1111                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1112                         printf("\n");
1113                 }
1114         } // geometry pointer did not exist before
1115
1116   if(fLoadGeomMatrices)
1117   {
1118     printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
1119
1120     // OADB if available
1121     AliOADBContainer emcGeoMat("AliEMCALgeo");
1122     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1123     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1124
1125
1126     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1127     {
1128
1129       if (!fGeomMatrix[mod]) // Get it from OADB
1130       {
1131         if(fDebug > 1 )
1132           printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1133                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1134         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1135
1136         fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1137       }
1138
1139       if(fGeomMatrix[mod])
1140       {
1141         if(DebugLevel() > 1)
1142           fGeomMatrix[mod]->Print();
1143
1144         fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1145       }
1146
1147       fGeomMatrixSet=kTRUE;
1148
1149     }//SM loop
1150   }//Load matrices
1151   else if(!gGeoManager)
1152   {
1153     printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1154     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1155     if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1156     {
1157       if(DebugLevel() > 1)
1158         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1159     }//AOD
1160     else
1161     {
1162       if(DebugLevel() > 1)
1163         printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1164
1165       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1166
1167       if(!esd)
1168       {
1169         Error("InitGeometry"," - This event does not contain ESDs?");
1170         if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1171         return;
1172       }
1173
1174       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1175       {
1176         if(DebugLevel() > 1)
1177           esd->GetEMCALMatrix(mod)->Print();
1178
1179         if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1180
1181       }
1182
1183       fGeomMatrixSet=kTRUE;
1184
1185     }//ESD
1186   }//Load matrices from Data
1187
1188 }
1189
1190 //____________________________________________________
1191 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1192 {
1193
1194   // Check if event is exotic, get an exotic cell and compare with triggered patch
1195   // If there is a match remove event ... to be completed, filled with something provisional
1196
1197   if(!fRemoveExoticEvents) return kFALSE;
1198
1199   // Loop on cells
1200   AliVCaloCells * cells = fEvent->GetEMCALCells();
1201   Float_t totCellE = 0;
1202   Int_t bc = InputEvent()->GetBunchCrossNumber();
1203   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1204   {
1205
1206     Float_t  ecell = 0 ;
1207     Double_t tcell = 0 ;
1208
1209     Int_t absID   = cells->GetCellNumber(icell);
1210     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1211     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1212   }
1213
1214   //  TString triggerclasses = "";
1215   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
1216   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
1217   //    //
1218   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1219   //    if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1220   //    return;
1221   //
1222
1223   //printf("TotE cell %f\n",totCellE);
1224   if(totCellE < 1) return kTRUE;
1225
1226   return kFALSE;
1227
1228 }
1229
1230 //________________________________________________________________
1231 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1232 {
1233   //Check if event is LED
1234
1235   if(!fRemoveLEDEvents) return kFALSE;
1236
1237   //check events of LHC11a period
1238   if(run < 146858 || run > 146860) return kFALSE ;
1239
1240   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1241   Int_t ncellsSM3 = 0;
1242   AliVCaloCells * cells = fEvent->GetEMCALCells();
1243   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1244   {
1245     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1246   }
1247
1248   TString triggerclasses = "";
1249
1250   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
1251   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
1252   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
1253
1254   Int_t ncellcut = 21;
1255   if(triggerclasses.Contains("EMC")) ncellcut = 35;
1256
1257   if( ncellsSM3 >= ncellcut)
1258   {
1259     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d\n",(Int_t)Entry(),ncellsSM3);
1260     if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1261     return kTRUE;
1262   }
1263
1264   return kFALSE;
1265
1266 }
1267
1268 //_______________________________________________________
1269 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
1270 {
1271   // Restore clusters from recPoints
1272   // Cluster energy, global position, cells and their amplitude
1273   // fractions are restored
1274
1275   Int_t j = 0;
1276   for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1277   {
1278     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1279
1280     const Int_t ncells = recPoint->GetMultiplicity();
1281     Int_t ncellsTrue = 0;
1282
1283     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1284
1285     // cells and their amplitude fractions
1286     UShort_t   absIds[ncells];
1287     Double32_t ratios[ncells];
1288
1289     //For later check embedding
1290     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1291
1292     Float_t clusterE = 0;
1293     for (Int_t c = 0; c < ncells; c++)
1294     {
1295       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1296
1297       absIds[ncellsTrue] = digit->GetId();
1298       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1299
1300       // In case of unfolding, remove digits with unfolded energy too low
1301       if(fSelectCell)
1302       {
1303         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1304         {
1305           if(DebugLevel() > 1)
1306           {
1307             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1308                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1309           }
1310
1311           continue;
1312
1313         } // if cuts
1314       }// Select cells
1315
1316       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1317       ncellsTrue++;
1318       clusterE  +=recPoint->GetEnergiesList()[c];
1319
1320       // In case of embedding, fill ratio with amount of signal,
1321       if (aodIH && aodIH->GetMergeEvents())
1322       {
1323         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1324         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1325         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1326
1327         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1328         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1329         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1330         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1331
1332         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1333         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1334
1335       }//Embedding
1336
1337     }// cluster cell loop
1338
1339     if (ncellsTrue < 1)
1340     {
1341       if (DebugLevel() > 1)
1342         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1343                recPoint->GetEnergy(), ncells);
1344       continue;
1345     }
1346
1347     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1348
1349     if(clusterE <  fRecParam->GetClusteringThreshold())
1350     {
1351       if (DebugLevel()>1)
1352         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1353       continue;
1354     }
1355
1356     TVector3 gpos;
1357     Float_t g[3];
1358
1359     // calculate new cluster position
1360
1361     recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1362     recPoint->GetGlobalPosition(gpos);
1363     gpos.GetXYZ(g);
1364
1365     // create a new cluster
1366
1367     (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1368     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1369     j++;
1370     clus->SetType(AliVCluster::kEMCALClusterv1);
1371     clus->SetE(clusterE);
1372     clus->SetPosition(g);
1373     clus->SetNCells(ncellsTrue);
1374     clus->SetCellsAbsId(absIds);
1375     clus->SetCellsAmplitudeFraction(ratios);
1376     clus->SetChi2(-1); //not yet implemented
1377     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1378     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1379     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1380
1381     if(ncells == ncellsTrue)
1382     {
1383       Float_t elipAxis[2];
1384       recPoint->GetElipsAxis(elipAxis);
1385       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1386       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1387       clus->SetDispersion(recPoint->GetDispersion());
1388     }
1389     else if(fSelectCell)
1390     {
1391       // In case some cells rejected, in unfolding case, recalculate
1392       // shower shape parameters and position
1393       if(DebugLevel() > 1)
1394         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1395
1396       AliVCaloCells* cells = 0x0;
1397       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
1398       else                                  cells = InputEvent()->GetEMCALCells();
1399
1400       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1401       fRecoUtils->RecalculateClusterPID(clus);
1402       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1403
1404     }
1405
1406     //MC
1407     Int_t  parentMult = 0;
1408     Int_t *parentList = recPoint->GetParents(parentMult);
1409     clus->SetLabel(parentList, parentMult);
1410
1411     //Write the second major contributor to each MC cluster.
1412     Int_t iNewLabel ;
1413     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1414     {
1415
1416       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1417       if(idCell>=0)
1418       {
1419         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1420         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1421         {
1422           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1423           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1424         }
1425         if (iNewLabel == 1)
1426         {
1427           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1428           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1429           {
1430             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1431           }
1432
1433           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1434           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1435           delete [] newLabelArray;
1436         }
1437       }//positive cell number
1438     }
1439
1440   } // recPoints loop
1441
1442 }
1443
1444
1445 //____________________________________________________________
1446 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1447 {
1448   // Init geometry, create list of output clusters
1449
1450
1451   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1452
1453   if(fOutputAODBranchName.Length()==0)
1454   {
1455     fOutputAODBranchName = "newEMCALClustersArray";
1456     printf("Cluster branch name not set, set it to newEMCALClustersArray \n");
1457   }
1458
1459   fOutputAODBranch->SetName(fOutputAODBranchName);
1460
1461   if( fFillAODFile )
1462   {
1463     //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1464
1465
1466     //fOutputAODBranch->SetOwner(kFALSE);
1467
1468     AddAODBranch("TClonesArray", &fOutputAODBranch);
1469   }
1470
1471 }
1472
1473 //_______________________________________________________
1474 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1475 {
1476   // Do clusterization event by event, execute different steps
1477   // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed, initializations
1478   //    load things and clear arrays
1479   // 2) Clusterize a) just unfolding existing clusters (fJustUnfold)
1480   //               b) recluster cells
1481   //                   +  convert cells into digits (calibrating them)
1482   //                   +  recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN)
1483   //                      with methods in AliEMCALClusterizer
1484   //                   +  transform recPoints into CaloClusters
1485   // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch
1486
1487   //-------
1488   // Step 1
1489
1490   //Remove the contents of AOD branch output list set in the previous event
1491   fOutputAODBranch->Clear("C");
1492
1493   LoadBranches();
1494
1495   //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
1496   //If we need a centrality bin, we select only those events in the corresponding bin.
1497   if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1498   {
1499     Float_t cen = GetEventCentrality();
1500     if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
1501   }
1502
1503   // intermediate array with new clusters : init the array only once or clear from previous event
1504   if(!fCaloClusterArr) fCaloClusterArr    = new TObjArray(10000);
1505   else                 fCaloClusterArr->Delete();//Clear("C"); it leaks?
1506
1507   //Get the event, do some checks and settings
1508   CheckAndGetEvent() ;
1509
1510   if (!fEvent)
1511   {
1512     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1513     return ;
1514   }
1515
1516   //Init pointers, geometry, clusterizer, ocdb, aodb
1517
1518   InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1519
1520   if(fAccessOCDB) AccessOCDB();
1521   if(fAccessOADB) AccessOADB(); // only once
1522
1523   InitClusterization();
1524
1525   //-------
1526   // Step 2
1527
1528   // Make clusters
1529   if (fJustUnfold) ClusterUnfolding();
1530   else             ClusterizeCells() ;
1531
1532   //-------
1533   // Step 3
1534
1535   FillCaloClusterInEvent();
1536
1537 }
1538