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