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