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