PWG4/CaloCalib/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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // This analysis provides a new list of clusters to be used in other analysis
  19 //
  20 // Author: Gustavo Conesa Balbastre,
  21 //         Adapted from analysis class from Deepa Thomas
  22 //
  23 //
  24 //_________________________________________________________________________
  25
  26 // --- Root ---
  27 #include "TString.h"
  28 #include "TRefArray.h"
  29 #include "TClonesArray.h"
  30 #include "TTree.h"
  31 #include "TGeoManager.h"
  32
  33 // --- AliRoot Analysis Steering
  34 #include "AliAnalysisTask.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliESDEvent.h"
  37 #include "AliGeomManager.h"
  38 #include "AliVCaloCells.h"
  39 #include "AliAODCaloCluster.h"
  40 #include "AliCDBManager.h"
  41 #include "AliCDBStorage.h"
  42 #include "AliCDBEntry.h"
  43 #include "AliLog.h"
  44 #include "AliVEventHandler.h"
  45
  46 // --- EMCAL
  47 #include "AliEMCALRecParam.h"
  48 #include "AliEMCALAfterBurnerUF.h"
  49 #include "AliEMCALGeometry.h"
  50 #include "AliEMCALClusterizerNxN.h"
  51 #include "AliEMCALClusterizerv1.h"
  52 #include "AliEMCALRecPoint.h"
  53 #include "AliEMCALDigit.h"
  54 #include "AliCaloCalibPedestal.h"
  55 #include "AliEMCALCalibData.h"
  56 #include "AliEMCALRecoUtils.h"
  57
  58 #include "AliAnalysisTaskEMCALClusterize.h"
  59
  60 ClassImp(AliAnalysisTaskEMCALClusterize)
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
  64   : AliAnalysisTaskSE(name)
  65   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  66   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  67   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  68   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  69   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kTRUE)
  70   , fRun(-1),            fRecoUtils(0)
  71
  72   {
  73   //ctor
  74   for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
  75   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  76   fClusterArr      = new TObjArray(100);
  77   fCaloClusterArr  = new TObjArray(100);
  78   fRecParam        = new AliEMCALRecParam;
  79   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
  80   fRecoUtils       = new AliEMCALRecoUtils();
  81 }
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
  85   : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
  86   , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"),  fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
  87   , fCalibData(0),       fPedestalData(0),     fOCDBpath("raw://")
  88   , fDigitsArr(0),       fClusterArr(0),       fCaloClusterArr(0)
  89   , fRecParam(0),        fClusterizer(0),      fUnfolder(0),           fJustUnfold(kFALSE)
  90   , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"),     fFillAODFile(kFALSE)
  91   , fRun(-1),            fRecoUtils(0)
  92 {
  93   // Constructor
  94   for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
  95   fDigitsArr       = new TClonesArray("AliEMCALDigit",200);
  96   fClusterArr      = new TObjArray(100);
  97   fCaloClusterArr  = new TObjArray(100);
  98   fRecParam        = new AliEMCALRecParam;
  99   fBranchNames     = "ESD:AliESDHeader.,EMCALCells.";
 100   fRecoUtils       = new AliEMCALRecoUtils();
 101 }
 102
 103 //________________________________________________________________________
 104 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
 105 {
 106   //dtor
 107
 108   if (fDigitsArr){
 109     fDigitsArr->Clear("C");
 110     delete fDigitsArr;
 111   }
 112
 113   if (fClusterArr){
 114     fClusterArr->Delete();
 115     delete fClusterArr;
 116   }
 117
 118   if (fCaloClusterArr){
 119     fCaloClusterArr->Delete();
 120     delete fCaloClusterArr;
 121   }
 122
 123   if(fClusterizer) delete fClusterizer;
 124   if(fUnfolder)    delete fUnfolder;
 125   if(fRecoUtils)   delete fRecoUtils;
 126
 127 }
 128
 129 //-------------------------------------------------------------------
 130 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
 131 {
 132   // Init geometry, create list of output clusters
 133
 134   fGeom =  AliEMCALGeometry::GetInstance(fGeomName) ;
 135
 136   fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 137   fOutputAODBranch->SetName(fOutputAODBranchName);
 138   AddAODBranch("TClonesArray", &fOutputAODBranch);
 139   Info("UserCreateOutputObjects","Create Branch: %s",fOutputAODBranchName.Data());
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
 144 {
 145   // Main loop
 146   // Called for each event
 147
 148   //Remove the contents of output list set in the previous event
 149   fOutputAODBranch->Clear("C");
 150
 151   AliVEvent   * event    = InputEvent();
 152   AliESDEvent * esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 153
 154   if (!event) {
 155     Error("UserExec","Event not available");
 156     return;
 157   }
 158
 159   //Magic line to write events to AOD file
 160   AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
 161   LoadBranches();
 162
 163   AccessOCDB();
 164
 165   //-------------------------------------------------------------------------------------
 166   //Set the geometry matrix, for the first event, skip the rest
 167   //-------------------------------------------------------------------------------------
 168   if(!fGeomMatrixSet){
 169     if(fLoadGeomMatrices){
 170       for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 171         if(fGeomMatrix[mod]){
 172           if(DebugLevel() > 1)
 173             fGeomMatrix[mod]->Print();
 174           fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
 175         }
 176         fGeomMatrixSet=kTRUE;
 177       }//SM loop
 178     }//Load matrices
 179     else if(!gGeoManager){
 180       Info("UserExec","Get geo matrices from data");
 181       //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
 182       if(!strcmp(event->GetName(),"AliAODEvent")) {
 183         if(DebugLevel() > 1)
 184           Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
 185       }//AOD
 186       else {
 187         if(DebugLevel() > 1)
 188           Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
 189         AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
 190         if(!esd) {
 191           Error("UserExec","This event does not contain ESDs?");
 192           return;
 193         }
 194         for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
 195           if(DebugLevel() > 1)
 196             esd->GetEMCALMatrix(mod)->Print();
 197           if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
 198         }
 199         fGeomMatrixSet=kTRUE;
 200       }//ESD
 201     }//Load matrices from Data
 202   }//first event
 203
 204   //Get the list of cells needed for unfolding and reclustering
 205   AliVCaloCells *cells= event->GetEMCALCells();
 206   Int_t nClustersOrg = 0;
 207   //-------------------------------------------
 208   //---------Unfolding clusters----------------
 209   //-------------------------------------------
 210   if (fJustUnfold) {
 211
 212     //Fill the array with the EMCAL clusters, copy them
 213     for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
 214     {
 215       AliVCluster *clus = event->GetCaloCluster(i);
 216       if(clus->IsEMCAL()){
 217         //printf("Org Cluster %d, E %f\n",i, clus->E());
 218         AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
 219         AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
 220         if     (esdCluster){
 221           fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
 222         }//ESD
 223         else if(aodCluster){
 224           fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
 225         }//AOD
 226         else
 227           Warning("UserExec()"," - Wrong CaloCluster type?");
 228         nClustersOrg++;
 229       }
 230     }
 231
 232     //Do the unfolding
 233     fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
 234
 235     //CLEAN-UP
 236     fUnfolder->Clear();
 237
 238     //printf("nClustersOrg %d\n",nClustersOrg);
 239   }
 240   //-------------------------------------------
 241   //---------- Recluster cells ----------------
 242   //-------------------------------------------
 243
 244   else{
 245     //-------------------------------------------------------------------------------------
 246     //Transform CaloCells into Digits
 247     //-------------------------------------------------------------------------------------
 248     Int_t idigit =  0;
 249     Int_t id     = -1;
 250     Float_t amp  = -1;
 251     Float_t time = -1;
 252
 253     Double_t cellAmplitude = 0;
 254     Double_t cellTime      = 0;
 255     Short_t  cellNumber    = 0;
 256
 257     TTree *digitsTree = new TTree("digitstree","digitstree");
 258     digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
 259
 260     for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
 261     {
 262       if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 263         break;
 264
 265       time = cellTime;
 266       amp  = cellAmplitude;
 267       id   = cellNumber;
 268
 269       AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
 270       digit->SetId(id);
 271       digit->SetAmplitude(amp);
 272       digit->SetTime(time);
 273       digit->SetTimeR(time);
 274       digit->SetIndexInList(idigit);
 275       digit->SetType(AliEMCALDigit::kHG);
 276       //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
 277       idigit++;
 278     }
 279
 280     //Fill the tree with digits
 281     digitsTree->Fill();
 282
 283     //-------------------------------------------------------------------------------------
 284     //Do the clusterization
 285     //-------------------------------------------------------------------------------------
 286     TTree *clustersTree = new TTree("clustertree","clustertree");
 287
 288     fClusterizer->SetInput(digitsTree);
 289     fClusterizer->SetOutput(clustersTree);
 290     fClusterizer->Digits2Clusters("");
 291
 292     //-------------------------------------------------------------------------------------
 293     //Transform the recpoints into AliVClusters
 294     //-------------------------------------------------------------------------------------
 295
 296     clustersTree->SetBranchStatus("*",0); //disable all branches
 297     clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
 298
 299     TBranch *branch = clustersTree->GetBranch("EMCALECARP");
 300     branch->SetAddress(&fClusterArr);
 301     branch->GetEntry(0);
 302
 303     RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
 304
 305     //---CLEAN UP-----
 306     fClusterizer->Clear();
 307     fDigitsArr  ->Clear("C");
 308     fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
 309
 310     clustersTree->Delete("all");
 311     digitsTree  ->Delete("all");
 312   }
 313
 314   //Recalculate track-matching for the new clusters, only with ESDs
 315   if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
 316
 317
 318   //-------------------------------------------------------------------------------------
 319   //Put the new clusters in the AOD list
 320   //-------------------------------------------------------------------------------------
 321
 322   Int_t kNumberOfCaloClusters   = fCaloClusterArr->GetEntries();
 323   //Info("UserExec","New clusters %d",kNumberOfCaloClusters);
 324   //if(nClustersOrg!=kNumberOfCaloClusters) Info("UserExec","Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
 325   for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
 326     AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
 327     //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
 328
 329     //Add matched track, if any, only with ESDs
 330     if(esdevent){
 331       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 332       if(trackIndex >= 0){
 333         newCluster->AddTrackMatched(event->GetTrack(trackIndex));
 334         if(DebugLevel() > 1)
 335           Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
 336       }
 337     }
 338     new((*fOutputAODBranch)[i])  AliAODCaloCluster(*newCluster);
 339   }
 340
 341   //---CLEAN UP-----
 342   fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
 343 }
 344
 345 //_____________________________________________________________________
 346 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
 347 {
 348   //Access to OCDB stuff
 349
 350   AliVEvent * event = InputEvent();
 351   if (!event)
 352   {
 353     Warning("AccessOCDB","Event not available!!!");
 354     return kFALSE;
 355   }
 356
 357   if (event->GetRunNumber()==fRun)
 358     return kTRUE;
 359   fRun = event->GetRunNumber();
 360
 361   if(DebugLevel() > 1 )
 362     Info("AccessODCD()"," Begin");
 363
 364   fGeom = AliEMCALGeometry::GetInstance(fGeomName);
 365
 366
 367   AliCDBManager *cdb = AliCDBManager::Instance();
 368
 369
 370   if (fOCDBpath.Length()){
 371     cdb->SetDefaultStorage(fOCDBpath.Data());
 372     Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
 373   }
 374
 375   cdb->SetRun(event->GetRunNumber());
 376
 377   //
 378   // EMCAL from RAW OCDB
 379   if (fOCDBpath.Contains("alien:"))
 380   {
 381     cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
 382     cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
 383   }
 384
 385   TString path = cdb->GetDefaultStorage()->GetBaseFolder();
 386
 387   // init parameters:
 388
 389   //Get calibration parameters
 390   if(!fCalibData)
 391   {
 392     AliCDBEntry *entry = (AliCDBEntry*)
 393       AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 394
 395     if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 396   }
 397
 398   if(!fCalibData)
 399     AliFatal("Calibration parameters not found in CDB!");
 400
 401   //Get calibration parameters
 402   if(!fPedestalData)
 403   {
 404     AliCDBEntry *entry = (AliCDBEntry*)
 405       AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
 406
 407     if (entry) fPedestalData =  (AliCaloCalibPedestal*) entry->GetObject();
 408   }
 409
 410   if(!fPedestalData)
 411     AliFatal("Dead map not found in CDB!");
 412
 413   InitClusterization();
 414
 415   return kTRUE;
 416 }
 417
 418 //________________________________________________________________________________________
 419 void AliAnalysisTaskEMCALClusterize::InitClusterization()
 420 {
 421   //Select clusterization/unfolding algorithm and set all the needed parameters
 422
 423   if (fJustUnfold){
 424     // init the unfolding afterburner
 425     delete fUnfolder;
 426     fUnfolder =  new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 427     return;
 428  }
 429
 430   //First init the clusterizer
 431   delete fClusterizer;
 432   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 433     fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
 434   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 435     fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 436   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
 437    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
 438    clusterizer->SetNRowDiff(2);
 439    clusterizer->SetNColDiff(2);
 440    fClusterizer = clusterizer;
 441   } else{
 442     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 443   }
 444
 445   //Now set the parameters
 446   fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
 447   fClusterizer->SetECALogWeight          ( fRecParam->GetW0()                  );
 448   fClusterizer->SetMinECut               ( fRecParam->GetMinECut()             );
 449   fClusterizer->SetUnfolding             ( fRecParam->GetUnfold()              );
 450   fClusterizer->SetECALocalMaxCut        ( fRecParam->GetLocMaxCut()           );
 451   fClusterizer->SetTimeCut               ( fRecParam->GetTimeCut()             );
 452   fClusterizer->SetTimeMin               ( fRecParam->GetTimeMin()             );
 453   fClusterizer->SetTimeMax               ( fRecParam->GetTimeMax()             );
 454   fClusterizer->SetInputCalibrated       ( kTRUE                               );
 455
 456   //In case of unfolding after clusterization is requested, set the corresponding parameters
 457   if(fRecParam->GetUnfold()){
 458     Int_t i=0;
 459     for (i = 0; i < 8; i++) {
 460       fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
 461     }//end of loop over parameters
 462     for (i = 0; i < 3; i++) {
 463       fClusterizer->SetPar5  (i, fRecParam->GetPar5(i));
 464       fClusterizer->SetPar6  (i, fRecParam->GetPar6(i));
 465     }//end of loop over parameters
 466
 467     fClusterizer->InitClusterUnfolding();
 468   }// to unfold
 469 }
 470
 471 //________________________________________________________________________________________
 472 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
 473 {
 474   // Restore clusters from recPoints
 475   // Cluster energy, global position, cells and their amplitude fractions are restored
 476
 477   for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
 478   {
 479     AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
 480
 481     const Int_t ncells = recPoint->GetMultiplicity();
 482     Int_t ncells_true = 0;
 483
 484     // cells and their amplitude fractions
 485     UShort_t   absIds[ncells];
 486     Double32_t ratios[ncells];
 487
 488     for (Int_t c = 0; c < ncells; c++) {
 489       AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
 490
 491       absIds[ncells_true] = digit->GetId();
 492       ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
 493
 494       if (ratios[ncells_true] > 0.001) ncells_true++;
 495     }
 496
 497     if (ncells_true < 1) {
 498       Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
 499       continue;
 500     }
 501
 502     TVector3 gpos;
 503     Float_t g[3];
 504
 505     // calculate new cluster position
 506     recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
 507     recPoint->GetGlobalPosition(gpos);
 508     gpos.GetXYZ(g);
 509
 510     // create a new cluster
 511     AliAODCaloCluster *clus = new AliAODCaloCluster();
 512     clus->SetType(AliVCluster::kEMCALClusterv1);
 513     clus->SetE(recPoint->GetEnergy());
 514     clus->SetPosition(g);
 515     clus->SetNCells(ncells_true);
 516     clus->SetCellsAbsId(absIds);
 517     clus->SetCellsAmplitudeFraction(ratios);
 518     clus->SetDispersion(recPoint->GetDispersion());
 519     clus->SetChi2(-1); //not yet implemented
 520     clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
 521     clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
 522     Float_t elipAxis[2];
 523     recPoint->GetElipsAxis(elipAxis);
 524     clus->SetM02(elipAxis[0]*elipAxis[0]) ;
 525     clus->SetM20(elipAxis[1]*elipAxis[1]) ;
 526     clus->SetDistToBadChannel(recPoint->GetDistanceToBadTower());
 527     clusArray->Add(clus);
 528   } // recPoints loop
 529 }