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