PWG4/CaloCalib/AliAnalysisTaskEMCALClusterizeFast.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 /* $Id$ */
  17
  18 // --- Root ---
  19 #include <TClonesArray.h>
  20 #include <TGeoManager.h>
  21 #include <TObjArray.h>
  22 #include <TString.h>
  23 #include <TTree.h>
  24 #include "AliAODCaloCluster.h"
  25 #include "AliAODEvent.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliCDBEntry.h"
  28 #include "AliCDBManager.h"
  29 #include "AliCaloCalibPedestal.h"
  30 #include "AliEMCALAfterBurnerUF.h"
  31 #include "AliEMCALCalibData.h"
  32 #include "AliEMCALClusterizerNxN.h"
  33 #include "AliEMCALClusterizerv1.h"
  34 #include "AliEMCALDigit.h"
  35 #include "AliEMCALGeometry.h"
  36 #include "AliEMCALRecParam.h"
  37 #include "AliEMCALRecPoint.h"
  38 #include "AliEMCALRecoUtils.h"
  39 #include "AliESDEvent.h"
  40 #include "AliLog.h"
  41
  42 #include "AliAnalysisTaskEMCALClusterizeFast.h"
  43
  44 ClassImp(AliAnalysisTaskEMCALClusterizeFast)
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast()
  48   : AliAnalysisTaskSE(),
  49     fRun(0),
  50     fDigitsArr(0),
  51     fClusterArr(0),
  52     fRecParam(0),
  53     fClusterizer(0),
  54     fUnfolder(0),
  55     fJustUnfold(kFALSE),
  56     fGeomName(),
  57     fGeomMatrixSet(kFALSE),
  58     fLoadGeomMatrices(kFALSE),
  59     fOCDBpath(),
  60     fCalibData(0),
  61     fPedestalData(0),
  62     fOutputAODBranch(0),
  63     fOutputAODBrName(),
  64     fRecoUtils(0),
  65     fLoadCalib(0),
  66     fLoadPed(0),
  67     fAttachClusters(0),
  68     fRecalibOnly(0),
  69     fSubBackground(0)
  70 {
  71   // Constructor
  72 }
  73
  74 //________________________________________________________________________
  75 AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
  76   : AliAnalysisTaskSE(name),
  77     fRun(0),
  78     fDigitsArr(0),
  79     fClusterArr(0),
  80     fRecParam(new AliEMCALRecParam),
  81     fClusterizer(0),
  82     fUnfolder(0),
  83     fJustUnfold(kFALSE),
  84     fGeomName("EMCAL_FIRSTYEARV1"),
  85     fGeomMatrixSet(kFALSE),
  86     fLoadGeomMatrices(kFALSE),
  87     fOCDBpath(),
  88     fCalibData(0),
  89     fPedestalData(0),
  90     fOutputAODBranch(0),
  91     fOutputAODBrName(),
  92     fRecoUtils(0),
  93     fLoadCalib(0),
  94     fLoadPed(0),
  95     fAttachClusters(0),
  96     fRecalibOnly(0),
  97     fSubBackground(0)
  98 {
  99   // Constructor
 100
 101   fBranchNames     = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
 102   for(Int_t i = 0; i < 12; ++i)
 103     fGeomMatrix[i] = 0;
 104 }
 105
 106 //________________________________________________________________________
 107 AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
 108 {
 109   // Destructor.
 110
 111   delete fDigitsArr;
 112   delete fClusterizer;
 113   delete fUnfolder;
 114   delete fRecoUtils;
 115 }
 116
 117 //-------------------------------------------------------------------
 118 void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
 119 {
 120   // Create output objects.
 121
 122   if (!fOutputAODBrName.IsNull()) {
 123     fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
 124     fOutputAODBranch->SetName(fOutputAODBrName);
 125     AddAODBranch("TClonesArray", &fOutputAODBranch);
 126     AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
 127   }
 128 }
 129
 130 //________________________________________________________________________
 131 void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
 132 {
 133   // Main loop, called for each event
 134
 135   // remove the contents of output list set in the previous event
 136   if (fOutputAODBranch)
 137     fOutputAODBranch->Clear("C");
 138
 139   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 140   AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
 141
 142   if (!esdevent&&!aodevent) {
 143     Error("UserExec","Event not available");
 144     return;
 145   }
 146
 147   LoadBranches();
 148
 149   Init();
 150
 151   if (fJustUnfold) {
 152     AliWarning("Unfolding not implemented");
 153     return;
 154   }
 155
 156   FillDigitsArray();
 157
 158   if (fRecalibOnly) {
 159     UpdateCells();
 160     return; // not requested to run clusterizer
 161   }
 162
 163   Clusterize();
 164   UpdateCells();
 165   UpdateClusters();
 166
 167   if (fOutputAODBranch)
 168     RecPoints2Clusters(fOutputAODBranch);
 169 }
 170
 171 //________________________________________________________________________
 172 void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
 173 {
 174   // Clusterize
 175
 176   if (fSubBackground) {
 177     fClusterizer->SetInputCalibrated(kTRUE);
 178     fClusterizer->SetCalibrationParameters(0);
 179   }
 180   fClusterizer->Digits2Clusters("");
 181   if (fSubBackground) {
 182     fClusterizer->SetInputCalibrated(kFALSE);
 183     fClusterizer->SetCalibrationParameters(fCalibData);
 184   }
 185 }
 186
 187 //________________________________________________________________________
 188 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
 189 {
 190   // Fill digits from cells.
 191
 192   fDigitsArr->Clear("C");
 193   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 194   Double_t avgE = 0; // for background subtraction
 195   Int_t ncells = cells->GetNumberOfCells();
 196   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 197     Double_t cellAmplitude=0, cellTime=0;
 198     Short_t cellNumber=0;
 199     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 200       break;
 201     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 202     digit->SetId(cellNumber);
 203     digit->SetTime(cellTime);
 204     digit->SetTimeR(cellTime);
 205     digit->SetIndexInList(idigit);
 206     digit->SetType(AliEMCALDigit::kHG);
 207     if (fRecalibOnly||fSubBackground) {
 208       Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
 209       digit->SetAmplitude(energy);
 210       avgE += energy;
 211     } else {
 212       digit->SetAmplitude(cellAmplitude);
 213     }
 214     idigit++;
 215   }
 216
 217   if (fSubBackground) {
 218     avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
 219     Int_t ndigis = fDigitsArr->GetEntries();
 220     for (Int_t i = 0; i < ndigis; ++i) {
 221       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
 222       Double_t energy = digit->GetAmplitude() - avgE;
 223       if (energy<=0.001) {
 224         //fDigitsArr->RemoveAt(i);
 225         digit->SetAmplitude(0);
 226       } else {
 227         digit->SetAmplitude(energy);
 228       }
 229     }
 230   }
 231   fDigitsArr->Compress();
 232   fDigitsArr->Sort();
 233 }
 234
 235 //________________________________________________________________________________________
 236 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
 237 {
 238   // Cluster energy, global position, cells and their amplitude fractions are restored.
 239
 240   Bool_t esdobjects = 0;
 241   if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
 242     esdobjects = 1;
 243
 244   Int_t Ncls = fClusterArr->GetEntriesFast();
 245   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
 246     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 247     Int_t ncells_true = 0;
 248     const Int_t ncells = recpoint->GetMultiplicity();
 249     UShort_t   absIds[ncells];
 250     Double32_t ratios[ncells];
 251     Int_t *dlist = recpoint->GetDigitsList();
 252     Float_t *elist = recpoint->GetEnergiesList();
 253     for (Int_t c = 0; c < ncells; ++c) {
 254       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 255       absIds[ncells_true] = digit->GetId();
 256       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
 257       if (ratios[ncells_true] < 0.001)
 258         continue;
 259       ++ncells_true;
 260     }
 261
 262     if (ncells_true < 1) {
 263       AliWarning("Skipping cluster with no cells");
 264       continue;
 265     }
 266
 267     // calculate new cluster position
 268     TVector3 gpos;
 269     recpoint->GetGlobalPosition(gpos);
 270     Float_t g[3];
 271     gpos.GetXYZ(g);
 272
 273     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
 274     c->SetType(AliVCluster::kEMCALClusterv1);
 275     c->SetE(recpoint->GetEnergy());
 276     c->SetPosition(g);
 277     c->SetNCells(ncells_true);
 278     c->SetDispersion(recpoint->GetDispersion());
 279     c->SetEmcCpvDistance(-1);            //not yet implemented
 280     c->SetChi2(-1);                      //not yet implemented
 281     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 282     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 283     Float_t elipAxis[2];
 284     recpoint->GetElipsAxis(elipAxis);
 285     c->SetM02(elipAxis[0]*elipAxis[0]) ;
 286     c->SetM20(elipAxis[1]*elipAxis[1]) ;
 287     c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
 288     if (esdobjects) {
 289       AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
 290       cesd->SetCellsAbsId(absIds);
 291       cesd->SetCellsAmplitudeFraction(ratios);
 292     } else {
 293       AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
 294       caod->SetCellsAbsId(absIds);
 295       caod->SetCellsAmplitudeFraction(ratios);
 296     }
 297   }
 298
 299   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 300   if (!esdevent)
 301     return;
 302   if (!fRecoUtils)
 303     return;
 304
 305   AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
 306   fRecoUtils->FindMatches(esdevent,clus);
 307
 308   if (!esdobjects) {
 309     Int_t Nclus = clus->GetEntries();
 310     for(Int_t i=0; i < Nclus; ++i) {
 311       AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
 312       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 313       if(trackIndex >= 0) {
 314         Float_t dR, dZ;
 315         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
 316         c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
 317         c->SetTrackDistance(dR,dZ); // not implemented
 318         c->SetEmcCpvDistance(dR);
 319         c->SetChi2(dZ);
 320         if(DebugLevel() > 1)
 321           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
 322       }
 323     }
 324   } else {
 325     Int_t Nclus = clus->GetEntries();
 326     for(Int_t i=0; i < Nclus; ++i) {
 327       AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
 328       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 329       if(trackIndex >= 0) {
 330         Float_t dR, dZ;
 331         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
 332         c->SetTrackDistance(dR,dZ);
 333         c->SetEmcCpvDistance(dR); //to be consistent with AODs
 334         c->SetChi2(dZ);           //to be consistent with AODs
 335         TArrayI tm(1,&trackIndex);
 336         c->AddTracksMatched(tm);
 337         if(DebugLevel() > 1)
 338           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
 339       }
 340     }
 341   }
 342 }
 343
 344 //________________________________________________________________________
 345 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
 346 {
 347   // Update cells in case re-calibration was done.
 348
 349   if (!fCalibData&&!fSubBackground)
 350     return;
 351
 352   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 353   Int_t ncells = cells->GetNumberOfCells();
 354   Int_t ndigis = fDigitsArr->GetEntries();
 355
 356   if (ncells!=ndigis) {
 357     cells->DeleteContainer();
 358     cells->CreateContainer(ndigis);
 359   }
 360   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
 361     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
 362     Double_t cellAmplitude = digit->GetCalibAmp();
 363     Short_t cellNumber = digit->GetId();
 364     Double_t cellTime = digit->GetTime();
 365     cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
 366   }
 367 }
 368
 369 //________________________________________________________________________
 370 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
 371 {
 372   // Update cells in case re-calibration was done.
 373
 374   if (!fAttachClusters)
 375     return;
 376
 377   TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
 378   if (!clus)
 379     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
 380   if(!clus)
 381     return;
 382
 383   Int_t nents = clus->GetEntries();
 384   for (Int_t i=0;i<nents;++i) {
 385     AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
 386     if (!c)
 387       continue;
 388     if (c->IsEMCAL())
 389       clus->RemoveAt(i);
 390   }
 391   clus->Compress();
 392   RecPoints2Clusters(clus);
 393 }
 394
 395 //________________________________________________________________________________________
 396 void AliAnalysisTaskEMCALClusterizeFast::Init()
 397 {
 398   //Select clusterization/unfolding algorithm and set all the needed parameters
 399
 400   AliVEvent * event = InputEvent();
 401   if (!event) {
 402     AliWarning("Event not available!!!");
 403     return;
 404   }
 405
 406   if (event->GetRunNumber()==fRun)
 407     return;
 408   fRun = event->GetRunNumber();
 409
 410   if (fJustUnfold){
 411     // init the unfolding afterburner
 412     delete fUnfolder;
 413     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 414     return;
 415   }
 416
 417   AliCDBManager *cdb = AliCDBManager::Instance();
 418   if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
 419     cdb->SetDefaultStorage(fOCDBpath);
 420   if (fRun!=cdb->GetRun())
 421     cdb->SetRun(fRun);
 422
 423   AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
 424   if (!geometry) {
 425     AliFatal("Geometry not available!!!");
 426     return;
 427   }
 428
 429   if (!fGeomMatrixSet) {
 430     if (fLoadGeomMatrices) {
 431       for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
 432         if(fGeomMatrix[mod]){
 433           if(DebugLevel() > 2)
 434             fGeomMatrix[mod]->Print();
 435           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
 436         }
 437       }
 438     } else { // get matrix from file (work around bug in aliroot)
 439       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
 440         const TGeoHMatrix *gm = 0;
 441         AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event->GetHeader());
 442         if (esdevent) {
 443           gm = esdevent->GetEMCALMatrix(mod);
 444         } else {
 445           AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
 446           if (aodheader) {
 447             gm = aodheader->GetEMCALMatrix(mod);
 448           }
 449         }
 450         if (gm) {
 451           if(DebugLevel() > 2)
 452             gm->Print();
 453           geometry->SetMisalMatrix(gm,mod);
 454         }
 455       }
 456     }
 457     fGeomMatrixSet=kTRUE;
 458   }
 459
 460   // setup digit array if needed
 461   if (!fDigitsArr) {
 462     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
 463     fDigitsArr->SetOwner(1);
 464   }
 465
 466   // then setup clusterizer
 467   delete fClusterizer;
 468   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 469     fClusterizer = new AliEMCALClusterizerv1(geometry);
 470   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 471     fClusterizer = new AliEMCALClusterizerNxN(geometry);
 472   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
 473    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
 474    clusterizer->SetNRowDiff(2);
 475    clusterizer->SetNColDiff(2);
 476    fClusterizer = clusterizer;
 477   } else{
 478     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 479   }
 480   fClusterizer->InitParameters(fRecParam);
 481   if (!fCalibData&&fLoadCalib) {
 482     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
 483     if (entry)
 484       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
 485     if (!fCalibData)
 486       AliFatal("Calibration parameters not found in CDB!");
 487   }
 488   if (!fPedestalData&&fLoadPed) {
 489     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
 490     if (entry)
 491       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
 492   }
 493   if (fCalibData) {
 494     fClusterizer->SetInputCalibrated(kFALSE);
 495     fClusterizer->SetCalibrationParameters(fCalibData);
 496     fClusterizer->SetCaloCalibPedestal(fPedestalData);
 497   } else {
 498     fClusterizer->SetInputCalibrated(kTRUE);
 499   }
 500   fClusterizer->SetJustClusters(kTRUE);
 501   fClusterizer->SetDigitsArr(fDigitsArr);
 502   fClusterizer->SetOutput(0);
 503   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 504 }