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   TFile* file = new TFile("tmpTree.root","recreate"); // Trick to avoid annoying messages FIXME
 515   TTree *digitsTree = new TTree("digitstree","digitstree");
 516   digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 517
 518   Int_t bc = InputEvent()->GetBunchCrossNumber();
 519   for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 520   {
 521
 522     // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
 523     id = cells->GetCellNumber(icell);
 524     Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
 525
 526     // Do not include cells with too low energy, nor exotic cell
 527     if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
 528
 529     // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
 530     if (time*1e9 < 1.)
 531     {
 532       time = fCellTime[id];
 533       //printf("cell %d time org %f - ",id, time*1.e9);
 534       fRecoUtils->RecalibrateCellTime(id,bc,time);
 535       //printf("recal %f\n",time*1.e9);
 536     }
 537
 538     if(  accept && fRecoUtils->IsExoticCell(id,cells,bc))
 539     {
 540       accept = kFALSE;
 541     }
 542
 543     if( !accept )
 544     {
 545       fCellLabels[id]      =-1; //reset the entry in the array for next event
 546       fCellSecondLabels[id]=-1; //reset the entry in the array for next event
 547       fCellTime[id]        = 0.;
 548       fCellMatchdEta[id]   =-999;
 549       fCellMatchdPhi[id]   =-999;
 550       if( DebugLevel() > 2 )
 551         printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
 552                                     id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
 553       continue;
 554     }
 555
 556     //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
 557     new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
 558
 559     fCellLabels[id]      =-1; //reset the entry in the array for next event
 560
 561     idigit++;
 562   }
 563
 564   //Fill the tree with digits
 565   digitsTree->Fill();
 566
 567   //-------------------------------------------------------------------------------------
 568   //Do the clusterization
 569   //-------------------------------------------------------------------------------------
 570   TTree *clustersTree = new TTree("clustertree","clustertree");
 571
 572   fClusterizer->SetInput(digitsTree);
 573   fClusterizer->SetOutput(clustersTree);
 574   fClusterizer->Digits2Clusters("");
 575
 576   //-------------------------------------------------------------------------------------
 577   //Transform the recpoints into AliVClusters
 578   //-------------------------------------------------------------------------------------
 579
 580   clustersTree->SetBranchStatus("*",0); //disable all branches
 581   clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 582
 583   TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 584   branch->SetAddress(&fClusterArr);
 585
 586   branch->GetEntry(0);
 587
 588   RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 589
 590   if(!fCaloClusterArr)
 591   {
 592     printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
 593     return;
 594   }
 595
 596   if( DebugLevel() > 0 )
 597   {
 598     printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
 599
 600     if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
 601       printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
 602              fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
 603              fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
 604     }
 605   }
 606
 607   //Reset the array with second labels for this event
 608   memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
 609
 610   //---CLEAN UP-----
 611   fClusterizer->Clear();
 612   fDigitsArr  ->Clear("C");
 613   fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 614
 615   clustersTree->Delete("all");
 616   digitsTree  ->Delete("all");
 617
 618   // Trick to avoid annoying messages FIXME
 619   file->Close();
 620   delete file;
 621 }
 622
 623 //_____________________________________________________
 624 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
 625 {
 626   // Take the event clusters and unfold them
 627
 628   AliVCaloCells *cells   = fEvent->GetEMCALCells();
 629   Double_t cellAmplitude = 0;
 630   Double_t cellTime      = 0;
 631   Short_t  cellNumber    = 0;
 632   Int_t    nClustersOrg  = 0;
 633
 634   // Fill the array with the EMCAL clusters, copy them
 635   for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
 636   {
 637     AliVCluster *clus = fEvent->GetCaloCluster(i);
 638     if(clus->IsEMCAL()){
 639
 640       //recalibrate/remove bad channels/etc if requested
 641       if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
 642         continue;
 643       }
 644
 645       if(fRecoUtils->IsRecalibrationOn()){
 646
 647         //Calibrate cluster
 648         fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
 649
 650         //CalibrateCells
 651         for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 652         {
 653           if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 654             break;
 655
 656           Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 657           fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
 658           fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 659
 660           //Do not include bad channels found in analysis?
 661           if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
 662              fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
 663             fCellLabels[cellNumber]      =-1; //reset the entry in the array for next event
 664             fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
 665             fCellTime[cellNumber]        = 0.;
 666             fCellMatchdEta[cellNumber]   =-999;
 667             fCellMatchdPhi[cellNumber]   =-999;
 668             continue;
 669           }
 670
 671           cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
 672
 673         }// cells loop
 674       }// recalibrate
 675
 676       //Cast to ESD or AOD, needed to create the cluster array
 677       AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 678       AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 679       if     (esdCluster){
 680         fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 681       }//ESD
 682       else if(aodCluster){
 683         fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 684       }//AOD
 685       else
 686         Warning("UserExec()"," - Wrong CaloCluster type?");
 687       nClustersOrg++;
 688     }
 689   }
 690
 691   //Do the unfolding
 692   fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 693
 694   //CLEAN-UP
 695   fUnfolder->Clear();
 696
 697 }
 698
 699 //_____________________________________________________
 700 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
 701 {
 702   // Put calo cells in standard branch
 703   AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
 704   Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
 705
 706   AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
 707   aodEMcells.CreateContainer(nEMcell);
 708   aodEMcells.SetType(AliVCaloCells::kEMCALCell);
 709   Double_t calibFactor = 1.;
 710   for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
 711     Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
 712     fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
 713     fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
 714
 715     if(fRecoUtils->IsRecalibrationOn()){
 716       calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
 717     }
 718
 719     if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
 720       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
 721     }
 722     else {
 723       aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
 724     }
 725   }
 726   aodEMcells.Sort();
 727
 728 }
 729
 730 //__________________________________________________
 731 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
 732 {
 733   //Put event header information in standard AOD branch
 734
 735   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
 736   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
 737
 738   Double_t pos[3]   ;
 739   Double_t covVtx[6];
 740   for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
 741
 742   AliAODHeader* header = AODEvent()->GetHeader();
 743   header->SetRunNumber(fEvent->GetRunNumber());
 744
 745   if(esdevent){
 746     TTree* tree = fInputHandler->GetTree();
 747     if (tree) {
 748       TFile* file = tree->GetCurrentFile();
 749       if (file) header->SetESDFileName(file->GetName());
 750     }
 751   }
 752   else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
 753
 754   header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
 755   header->SetOrbitNumber(fEvent->GetOrbitNumber());
 756   header->SetPeriodNumber(fEvent->GetPeriodNumber());
 757   header->SetEventType(fEvent->GetEventType());
 758
 759   //Centrality
 760   if(fEvent->GetCentrality()){
 761     header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
 762   }
 763   else{
 764     header->SetCentrality(0);
 765   }
 766
 767   //Trigger
 768   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
 769   if      (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
 770   else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
 771   header->SetTriggerMask(fEvent->GetTriggerMask());
 772   header->SetTriggerCluster(fEvent->GetTriggerCluster());
 773   if(esdevent){
 774     header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
 775     header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
 776     header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
 777   }
 778   else if (aodevent){
 779     header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
 780     header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
 781     header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
 782   }
 783
 784   header->SetMagneticField(fEvent->GetMagneticField());
 785   //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
 786
 787   header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
 788   header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
 789   header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
 790   header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
 791   header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
 792
 793   Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
 794   Float_t diamcov[3];
 795   fEvent->GetDiamondCovXY(diamcov);
 796   header->SetDiamond(diamxy,diamcov);
 797   if      (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
 798   else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
 799
 800   //
 801   //
 802   Int_t nVertices = 1 ;/* = prim. vtx*/;
 803   Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
 804
 805   AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
 806
 807   // Access to the AOD container of vertices
 808   TClonesArray &vertices = *(AODEvent()->GetVertices());
 809   Int_t jVertices=0;
 810
 811   // Add primary vertex. The primary tracks will be defined
 812   // after the loops on the composite objects (V0, cascades, kinks)
 813   fEvent->GetPrimaryVertex()->GetXYZ(pos);
 814   Float_t chi = 0;
 815   if      (esdevent){
 816     esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 817     chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
 818   }
 819   else if (aodevent){
 820     aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
 821     chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
 822   }
 823
 824   AliAODVertex * primary = new(vertices[jVertices++])
 825   AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
 826   primary->SetName(fEvent->GetPrimaryVertex()->GetName());
 827   primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
 828
 829 }
 830
 831 //_______________________________________________
 832 TString AliAnalysisTaskEMCALClusterize::GetPass()
 833 {
 834   // Get passx from filename.
 835
 836   if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
 837   {
 838     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
 839     return TString("");
 840   }
 841
 842   if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
 843   {
 844     AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
 845     return TString("");
 846   }
 847
 848   TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
 849   if      (pass.Contains("ass1")) return TString("pass1");
 850   else if (pass.Contains("ass2")) return TString("pass2");
 851   else if (pass.Contains("ass3")) return TString("pass3");
 852
 853   // No condition fullfilled, give a default value
 854   printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
 855   return TString("");
 856
 857 }
 858
 859 //_________________________________________
 860 void AliAnalysisTaskEMCALClusterize::Init()
 861 {
 862   //Init analysis with configuration macro if available
 863
 864   fOADBSet           = kFALSE;
 865   if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
 866
 867   fDigitsArr         = new TClonesArray("AliEMCALDigit",12000);
 868   fClusterArr        = new TObjArray(10000);
 869   fCaloClusterArr    = new TObjArray(10000);
 870   fBranchNames       = "ESD:AliESDHeader.,EMCALCells.";
 871
 872   if(!fRecParam)     fRecParam  = new AliEMCALRecParam;
 873   if(!fRecoUtils)    fRecoUtils = new AliEMCALRecoUtils();
 874
 875   if(fMaxEvent          <= 0) fMaxEvent          = 1000000000;
 876   if(fSelectCellMinE    <= 0) fSelectCellMinE    = 0.005;
 877   if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
 878
 879   if (fOCDBpath            == "") fOCDBpath            = "raw://" ;
 880   if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
 881
 882   if(gROOT->LoadMacro(fConfigName) >=0)
 883   {
 884     printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
 885     AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
 886     fGeomName         = clus->fGeomName;
 887     fLoadGeomMatrices = clus->fLoadGeomMatrices;
 888     fOCDBpath         = clus->fOCDBpath;
 889     fAccessOCDB       = clus->fAccessOCDB;
 890     fRecParam         = clus->fRecParam;
 891     fJustUnfold       = clus->fJustUnfold;
 892     fFillAODFile      = clus->fFillAODFile;
 893     fRecoUtils        = clus->fRecoUtils;
 894     fConfigName       = clus->fConfigName;
 895     fMaxEvent         = clus->fMaxEvent;
 896     fDoTrackMatching  = clus->fDoTrackMatching;
 897     fOutputAODBranchName = clus->fOutputAODBranchName;
 898     for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
 899   }
 900
 901   // Init geometry, I do not like much to do it like this ...
 902   if(fImportGeometryFromFile && !gGeoManager)
 903   {
 904     if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/EMCALTasks/macros/geometry.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
 905     printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
 906     TGeoManager::Import(fImportGeometryFilePath) ;
 907   }
 908
 909 }
 910
 911 //_______________________________________________________
 912 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 913 {
 914   //Select clusterization/unfolding algorithm and set all the needed parameters
 915
 916   if (fJustUnfold){
 917     // init the unfolding afterburner
 918     delete fUnfolder;
 919     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
 920     return;
 921   }
 922
 923   //First init the clusterizer
 924   delete fClusterizer;
 925   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 926     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 927   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 928     fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
 929   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){
 930     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 931     fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
 932     fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
 933   } else {
 934     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 935   }
 936
 937   //Now set the parameters
 938   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 939   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 940   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 941   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 942   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 943   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 944   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 945   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 946   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 947   fClusterizer->SetJustClusters          ( kTRUE                               );
 948   //In case of unfolding after clusterization is requested, set the corresponding parameters
 949   if(fRecParam->GetUnfold()){
 950     Int_t i=0;
 951     for (i = 0; i < 8; i++) {
 952       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 953     }//end of loop over parameters
 954     for (i = 0; i < 3; i++) {
 955       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 956       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 957     }//end of loop over parameters
 958
 959     fClusterizer->InitClusterUnfolding();
 960
 961   }// to unfold
 962 }
 963
 964 //_________________________________________________
 965 void AliAnalysisTaskEMCALClusterize::InitGeometry()
 966 {
 967   // Init geometry and set the geometry matrix, for the first event, skip the rest
 968   // Also set once the run dependent calibrations
 969
 970   if(fGeomMatrixSet) return;
 971
 972   Int_t runnumber = InputEvent()->GetRunNumber() ;
 973   if (!fGeom)
 974   {
 975     if(fGeomName=="")
 976     {
 977       if     (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
 978       else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
 979       else                        fGeomName = "EMCAL_COMPLETE12SMV1";
 980       printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
 981     }
 982
 983                 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 984
 985                 if(fDebug > 0)
 986     {
 987                         printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
 988                         if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
 989                         printf("\n");
 990                 }
 991         } // geometry pointer did not exist before
 992
 993   if(fLoadGeomMatrices)
 994   {
 995     printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
 996
 997     // OADB if available
 998     AliOADBContainer emcGeoMat("AliEMCALgeo");
 999     emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1000     TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1001
1002
1003     for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1004     {
1005
1006       if (!fGeomMatrix[mod]) // Get it from OADB
1007       {
1008         if(fDebug > 1 )
1009           printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1010                  mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1011         //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1012
1013         fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1014       }
1015
1016       if(fGeomMatrix[mod])
1017       {
1018         if(DebugLevel() > 1)
1019           fGeomMatrix[mod]->Print();
1020
1021         fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1022       }
1023
1024       fGeomMatrixSet=kTRUE;
1025
1026     }//SM loop
1027   }//Load matrices
1028   else if(!gGeoManager)
1029   {
1030     printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1031     //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1032     if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1033     {
1034       if(DebugLevel() > 1)
1035         Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1036     }//AOD
1037     else
1038     {
1039       if(DebugLevel() > 1)
1040         printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1041
1042       AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1043
1044       if(!esd)
1045       {
1046         Error("InitGeometry"," - This event does not contain ESDs?");
1047         AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1048         return;
1049       }
1050
1051       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1052       {
1053         if(DebugLevel() > 1)
1054           esd->GetEMCALMatrix(mod)->Print();
1055
1056         if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1057
1058       }
1059
1060       fGeomMatrixSet=kTRUE;
1061
1062     }//ESD
1063   }//Load matrices from Data
1064
1065
1066 }
1067
1068 //____________________________________________________
1069 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1070 {
1071
1072   // Check if event is exotic, get an exotic cell and compare with triggered patch
1073   // If there is a match remove event ... to be completed, filled with something provisional
1074
1075   if(!fRemoveExoticEvents) return kFALSE;
1076
1077   // Loop on cells
1078   AliVCaloCells * cells = fEvent->GetEMCALCells();
1079   Float_t totCellE = 0;
1080   Int_t bc = InputEvent()->GetBunchCrossNumber();
1081   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){
1082
1083     Float_t  ecell = 0 ;
1084     Double_t tcell = 0 ;
1085
1086     Int_t absID   = cells->GetCellNumber(icell);
1087     Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1088     if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1089   }
1090
1091   //  TString triggerclasses = "";
1092   //  if(esdevent) triggerclasses = esdevent             ->GetFiredTriggerClasses();
1093   //  else         triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
1094   //    //
1095   //    printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster  - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1096   //    AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1097   //    return;
1098   //
1099
1100   //printf("TotE cell %f\n",totCellE);
1101   if(totCellE < 1) return kTRUE;
1102
1103   return kFALSE;
1104
1105 }
1106
1107 //________________________________________________________________
1108 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1109 {
1110   //Check if event is LED
1111
1112   if(!fRemoveLEDEvents) return kFALSE;
1113
1114   //check events of LHC11a period
1115   if(run < 140000  || run > 146860) return kFALSE ;
1116
1117   // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1118   Int_t ncellsSM3 = 0;
1119   Int_t ncellsSM4 = 0;
1120   AliVCaloCells * cells = fEvent->GetEMCALCells();
1121   for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1122   {
1123     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1124     if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
1125   }
1126
1127   TString triggerclasses = "";
1128
1129   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
1130   if(esdevent) triggerclasses = esdevent              ->GetFiredTriggerClasses();
1131   else         triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
1132
1133   Int_t ncellcut = 21;
1134   if(triggerclasses.Contains("EMC")) ncellcut = 35;
1135
1136   if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
1137   {
1138     printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1139     AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1140     return kTRUE;
1141   }
1142
1143   return kFALSE;
1144
1145 }
1146
1147 //______________________________________________________________________________
1148 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr,
1149                                                         TObjArray *recPoints,
1150                                                         TObjArray *clusArray)
1151 {
1152   // Restore clusters from recPoints
1153   // Cluster energy, global position, cells and their amplitude fractions are restored
1154   Int_t j = 0;
1155   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
1156   {
1157     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
1158
1159     const Int_t ncells = recPoint->GetMultiplicity();
1160     Int_t ncellsTrue = 0;
1161
1162     if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1163
1164     // cells and their amplitude fractions
1165     UShort_t   absIds[ncells];
1166     Double32_t ratios[ncells];
1167
1168     //For later check embedding
1169     AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1170
1171     Float_t clusterE = 0;
1172     for (Int_t c = 0; c < ncells; c++)
1173     {
1174       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
1175
1176       absIds[ncellsTrue] = digit->GetId();
1177       ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1178
1179       // In case of unfolding, remove digits with unfolded energy too low
1180       if(fSelectCell){
1181         if     (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1182         {
1183           if(DebugLevel() > 1)
1184           {
1185             printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
1186                    recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1187           }
1188
1189           continue;
1190
1191         } // if cuts
1192       }// Select cells
1193
1194       //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1195       ncellsTrue++;
1196       clusterE  +=recPoint->GetEnergiesList()[c];
1197
1198       // In case of embedding, fill ratio with amount of signal,
1199       if (aodIH && aodIH->GetMergeEvents())
1200       {
1201         //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1202         AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1203         AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1204
1205         Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1206         //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1207         Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1208         //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1209
1210         if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1211         //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1212
1213       }//Embedding
1214
1215     }// cluster cell loop
1216
1217     if (ncellsTrue < 1)
1218     {
1219       if (DebugLevel() > 1)
1220         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1221                recPoint->GetEnergy(), ncells);
1222       continue;
1223     }
1224
1225     //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1226
1227     if(clusterE <  fRecParam->GetClusteringThreshold())
1228     {
1229       if (DebugLevel()>1)
1230         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1231       continue;
1232     }
1233
1234     TVector3 gpos;
1235     Float_t g[3];
1236
1237     // calculate new cluster position
1238     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
1239     recPoint->GetGlobalPosition(gpos);
1240     gpos.GetXYZ(g);
1241
1242     // create a new cluster
1243     (*clusArray)[j] = new AliAODCaloCluster() ;
1244     AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
1245     j++;
1246     clus->SetType(AliVCluster::kEMCALClusterv1);
1247     clus->SetE(clusterE);
1248     clus->SetPosition(g);
1249     clus->SetNCells(ncellsTrue);
1250     clus->SetCellsAbsId(absIds);
1251     clus->SetCellsAmplitudeFraction(ratios);
1252     clus->SetChi2(-1); //not yet implemented
1253     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1254     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1255     clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1256
1257     if(ncells == ncellsTrue)
1258     {
1259       Float_t elipAxis[2];
1260       recPoint->GetElipsAxis(elipAxis);
1261       clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1262       clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1263       clus->SetDispersion(recPoint->GetDispersion());
1264     }
1265     else if(fSelectCell)
1266     {
1267       // In case some cells rejected, in unfolding case, recalculate
1268       // shower shape parameters and position
1269       //if(DebugLevel() > 1)
1270         printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
1271       AliVCaloCells* cells = 0x0;
1272       if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent()  ->GetEMCALCells();
1273       else                                  cells = InputEvent()->GetEMCALCells();
1274       fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1275       fRecoUtils->RecalculateClusterPID(clus);
1276       fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1277
1278     }
1279
1280     //MC
1281     Int_t  parentMult = 0;
1282     Int_t *parentList = recPoint->GetParents(parentMult);
1283     clus->SetLabel(parentList, parentMult);
1284
1285     //Write the second major contributor to each MC cluster.
1286     Int_t iNewLabel ;
1287     for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1288     {
1289
1290       Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1291       if(idCell>=0){
1292         iNewLabel = 1 ; //iNewLabel makes sure we  don't write twice the same label.
1293         for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1294         {
1295           if ( fCellSecondLabels[idCell] == -1 )  iNewLabel = 0;  // -1 is never a good second label.
1296           if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) )  iNewLabel = 0;
1297         }
1298         if (iNewLabel == 1)
1299         {
1300           Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1301           for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1302           {
1303             newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1304           }
1305
1306           newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1307           clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1308           delete [] newLabelArray;
1309         }
1310       }//positive cell number
1311     }
1312
1313   } // recPoints loop
1314
1315 }
1316
1317
1318 //____________________________________________________________
1319 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1320 {
1321   // Init geometry, create list of output clusters
1322
1323   if(fOutputAODBranchName.Length()!=0)
1324   {
1325     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1326     fOutputAODBranch->SetName(fOutputAODBranchName);
1327     //fOutputAODBranch->SetOwner(kFALSE);
1328     AddAODBranch("TClonesArray", &fOutputAODBranch);
1329   }
1330   else
1331   {
1332     AliFatal("fOutputAODBranchName not set\n");
1333   }
1334
1335   //PostData(0,fOutputAODBranch);
1336
1337 }
1338
1339 //_______________________________________________________
1340 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1341 {
1342   // Main loop
1343   // Called for each event
1344
1345   //Remove the contents of output list set in the previous event
1346   fOutputAODBranch->Clear("C");
1347
1348   LoadBranches();
1349
1350   //Get the event, do some checks and settings
1351   CheckAndGetEvent() ;
1352
1353   if (!fEvent)
1354   {
1355     if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1356     return ;
1357   }
1358
1359   //Init pointers, geometry, clusterizer, ocdb, aodb
1360
1361   InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1362
1363   if(fAccessOCDB) AccessOCDB();
1364   if(fAccessOADB) AccessOADB(); // only once
1365
1366   InitClusterization();
1367
1368   // Make clusters
1369   if (fJustUnfold) ClusterUnfolding();
1370   else             ClusterizeCells() ;
1371
1372   //Recalculate track-matching for the new clusters
1373   if(fDoTrackMatching)
1374   {
1375     fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
1376   }
1377   //Put the new clusters in the AOD list
1378
1379   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntriesFast();
1380   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
1381     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1382
1383     newCluster->SetID(i);
1384
1385     //Add matched track
1386     if(fDoTrackMatching)
1387     {
1388       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1389       if(trackIndex >= 0)
1390       {
1391         newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1392         if(DebugLevel() > 1)
1393           printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
1394       }
1395
1396       Float_t dR = 999., dZ = 999.;
1397       fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
1398       newCluster->SetTrackDistance(dR,dZ);
1399
1400     }
1401     else
1402     {// Assign previously assigned matched track in reco, very very rough
1403       Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1404       newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1405     }
1406
1407     //printf("New cluster E %f, Time  %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1408     //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1409     //printf("\n");
1410
1411     // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1412     fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1413
1414     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
1415
1416     if(DebugLevel() > 1 )
1417       printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
1418
1419   } // cluster loop
1420
1421   fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1422
1423   // Clean up
1424   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
1425
1426 }
1427