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