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