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 mask1 = esdevent->GetESDRun()->GetDetectorsInDAQ();
 154     UInt_t mask2 = esdevent->GetESDRun()->GetDetectorsInReco();
 155     Bool_t desc1 = (mask1 >> 18) & 0x1;
 156     Bool_t desc2 = (mask2 >> 18) & 0x1;
 157     if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
 158       AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
 159                     mask1, esdevent->GetESDRun()->GetDetectorsInReco(),
 160                     mask2, esdevent->GetESDRun()->GetDetectorsInDAQ()));
 161       return;
 162     }
 163     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 164     offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
 165   } else if (aodevent) {
 166     offtrigger =  aodevent->GetHeader()->GetOfflineTrigger();
 167   }
 168   if (offtrigger & AliVEvent::kFastOnly) {
 169     AliWarning(Form("EMCAL not in fast only partition"));
 170     return;
 171   }
 172
 173   Init();
 174
 175   if (fJustUnfold) {
 176     AliWarning("Unfolding not implemented");
 177     return;
 178   }
 179
 180   FillDigitsArray();
 181
 182   if (fRecalibOnly) {
 183     UpdateCells();
 184     return; // not requested to run clusterizer
 185   }
 186
 187   Clusterize();
 188   UpdateCells();
 189   UpdateClusters();
 190
 191   if (fOutputAODBranch)
 192     RecPoints2Clusters(fOutputAODBranch);
 193 }
 194
 195 //________________________________________________________________________
 196 void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
 197 {
 198   // Clusterize
 199
 200   if (fSubBackground) {
 201     fClusterizer->SetInputCalibrated(kTRUE);
 202     fClusterizer->SetCalibrationParameters(0);
 203   }
 204   fClusterizer->Digits2Clusters("");
 205   if (fSubBackground) {
 206     if (fCalibData) {
 207       fClusterizer->SetInputCalibrated(kFALSE);
 208       fClusterizer->SetCalibrationParameters(fCalibData);
 209     }
 210   }
 211 }
 212
 213 //________________________________________________________________________
 214 void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
 215 {
 216   // Fill digits from cells.
 217
 218   fDigitsArr->Clear("C");
 219   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 220   Double_t avgE = 0; // for background subtraction
 221   Int_t ncells = cells->GetNumberOfCells();
 222   for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
 223     Double_t cellAmplitude=0, cellTime=0;
 224     Short_t cellNumber=0;
 225     if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
 226       break;
 227     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
 228     digit->SetId(cellNumber);
 229     digit->SetTime(cellTime);
 230     digit->SetTimeR(cellTime);
 231     digit->SetIndexInList(idigit);
 232     digit->SetType(AliEMCALDigit::kHG);
 233     if (fRecalibOnly||fSubBackground) {
 234       Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
 235       digit->SetAmplitude(energy);
 236       avgE += energy;
 237     } else {
 238       digit->SetAmplitude(cellAmplitude);
 239     }
 240     idigit++;
 241   }
 242
 243   if (fSubBackground) {
 244     avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
 245     Int_t ndigis = fDigitsArr->GetEntries();
 246     for (Int_t i = 0; i < ndigis; ++i) {
 247       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
 248       Double_t energy = digit->GetAmplitude() - avgE;
 249       if (energy<=0.001) {
 250         //fDigitsArr->RemoveAt(i);
 251         digit->SetAmplitude(0);
 252       } else {
 253         digit->SetAmplitude(energy);
 254       }
 255     }
 256   }
 257   fDigitsArr->Compress();
 258   fDigitsArr->Sort();
 259 }
 260
 261 //________________________________________________________________________________________
 262 void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
 263 {
 264   // Cluster energy, global position, cells and their amplitude fractions are restored.
 265
 266   Bool_t esdobjects = 0;
 267   if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
 268     esdobjects = 1;
 269
 270   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 271   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance(fGeomName);
 272
 273   Int_t Ncls = fClusterArr->GetEntriesFast();
 274   for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
 275     AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
 276     Int_t ncells_true = 0;
 277     const Int_t ncells = recpoint->GetMultiplicity();
 278     UShort_t   absIds[ncells];
 279     Double32_t ratios[ncells];
 280     Int_t *dlist = recpoint->GetDigitsList();
 281     Float_t *elist = recpoint->GetEnergiesList();
 282     for (Int_t c = 0; c < ncells; ++c) {
 283       AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
 284       absIds[ncells_true] = digit->GetId();
 285       ratios[ncells_true] = elist[c]/digit->GetAmplitude();
 286       if (ratios[ncells_true] < 0.001)
 287         continue;
 288       ++ncells_true;
 289     }
 290
 291     if (ncells_true < 1) {
 292       AliWarning("Skipping cluster with no cells");
 293       continue;
 294     }
 295
 296     // calculate new cluster position
 297     TVector3 gpos;
 298     recpoint->GetGlobalPosition(gpos);
 299     Float_t g[3];
 300     gpos.GetXYZ(g);
 301
 302     AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
 303     c->SetType(AliVCluster::kEMCALClusterv1);
 304     c->SetE(recpoint->GetEnergy());
 305     c->SetPosition(g);
 306     c->SetNCells(ncells_true);
 307     c->SetDispersion(recpoint->GetDispersion());
 308     c->SetEmcCpvDistance(-1);            //not yet implemented
 309     c->SetChi2(-1);                      //not yet implemented
 310     c->SetTOF(recpoint->GetTime()) ;     //time-of-flight
 311     c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
 312     Float_t elipAxis[2];
 313     recpoint->GetElipsAxis(elipAxis);
 314     c->SetM02(elipAxis[0]*elipAxis[0]) ;
 315     c->SetM20(elipAxis[1]*elipAxis[1]) ;
 316     if (fRecoUtils && fRecoUtils->IsBadChannelsRemovalSwitchedOn()) {
 317       fRecoUtils->RecalculateClusterDistanceToBadChannel(geom, cells, c);
 318     } else {
 319       if (fPedestalData)
 320         recpoint->EvalDistanceToBadChannels(fPedestalData);
 321       c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
 322     }
 323
 324     if (esdobjects) {
 325       AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
 326       cesd->SetCellsAbsId(absIds);
 327       cesd->SetCellsAmplitudeFraction(ratios);
 328     } else {
 329       AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
 330       caod->SetCellsAbsId(absIds);
 331       caod->SetCellsAmplitudeFraction(ratios);
 332     }
 333   }
 334
 335   AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
 336   if (!esdevent)
 337     return;
 338   if (!fRecoUtils)
 339     return;
 340
 341   AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
 342   fRecoUtils->FindMatches(esdevent,clus);
 343
 344   if (!esdobjects) {
 345     Int_t Nclus = clus->GetEntries();
 346     for(Int_t i=0; i < Nclus; ++i) {
 347       AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
 348       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 349       if(trackIndex >= 0) {
 350         Float_t dR, dZ;
 351         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
 352         c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
 353         c->SetTrackDistance(dR,dZ); // not implemented
 354         c->SetEmcCpvDistance(dR);
 355         c->SetChi2(dZ);
 356         if(DebugLevel() > 1)
 357           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
 358       }
 359     }
 360   } else {
 361     Int_t Nclus = clus->GetEntries();
 362     for(Int_t i=0; i < Nclus; ++i) {
 363       AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
 364       Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
 365       if(trackIndex >= 0) {
 366         Float_t dR, dZ;
 367         fRecoUtils->GetMatchedResiduals(i,dR, dZ);
 368         c->SetTrackDistance(dR,dZ);
 369         c->SetEmcCpvDistance(dR); //to be consistent with AODs
 370         c->SetChi2(dZ);           //to be consistent with AODs
 371         TArrayI tm(1,&trackIndex);
 372         c->AddTracksMatched(tm);
 373         if(DebugLevel() > 1)
 374           AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
 375       }
 376     }
 377   }
 378 }
 379
 380 //________________________________________________________________________
 381 void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
 382 {
 383   // Update cells in case re-calibration was done.
 384
 385   if (!fCalibData&&!fSubBackground)
 386     return;
 387
 388   AliVCaloCells *cells = InputEvent()->GetEMCALCells();
 389   Int_t ncells = cells->GetNumberOfCells();
 390   Int_t ndigis = fDigitsArr->GetEntries();
 391
 392   if (ncells!=ndigis) {
 393     cells->DeleteContainer();
 394     cells->CreateContainer(ndigis);
 395   }
 396   for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
 397     AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
 398     Double_t cellAmplitude = digit->GetCalibAmp();
 399     Short_t cellNumber = digit->GetId();
 400     Double_t cellTime = digit->GetTime();
 401     cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
 402   }
 403 }
 404
 405 //________________________________________________________________________
 406 void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
 407 {
 408   // Update cells in case re-calibration was done.
 409
 410   if (!fAttachClusters)
 411     return;
 412
 413   TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
 414   if (!clus)
 415     clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
 416   if(!clus)
 417     return;
 418
 419   Int_t nents = clus->GetEntries();
 420   for (Int_t i=0;i<nents;++i) {
 421     AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
 422     if (!c)
 423       continue;
 424     if (c->IsEMCAL())
 425       clus->RemoveAt(i);
 426   }
 427   clus->Compress();
 428   RecPoints2Clusters(clus);
 429 }
 430
 431 //________________________________________________________________________________________
 432 void AliAnalysisTaskEMCALClusterizeFast::Init()
 433 {
 434   //Select clusterization/unfolding algorithm and set all the needed parameters
 435
 436   AliVEvent * event = InputEvent();
 437   if (!event) {
 438     AliWarning("Event not available!!!");
 439     return;
 440   }
 441
 442   if (event->GetRunNumber()==fRun)
 443     return;
 444   fRun = event->GetRunNumber();
 445
 446   if (fJustUnfold){
 447     // init the unfolding afterburner
 448     delete fUnfolder;
 449     fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
 450     return;
 451   }
 452
 453   AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
 454   if (!geometry) {
 455     AliFatal("Geometry not available!!!");
 456     return;
 457   }
 458
 459   if (!fGeomMatrixSet) {
 460     if (fLoadGeomMatrices) {
 461       for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
 462         if(fGeomMatrix[mod]){
 463           if(DebugLevel() > 2)
 464             fGeomMatrix[mod]->Print();
 465           geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
 466         }
 467       }
 468     } else { // get matrix from file (work around bug in aliroot)
 469       for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
 470         const TGeoHMatrix *gm = 0;
 471         AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event);
 472         if (esdevent) {
 473           gm = esdevent->GetEMCALMatrix(mod);
 474         } else {
 475           AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
 476           if (aodheader) {
 477             gm = aodheader->GetEMCALMatrix(mod);
 478           }
 479         }
 480         if (gm) {
 481           if(DebugLevel() > 2)
 482             gm->Print();
 483           geometry->SetMisalMatrix(gm,mod);
 484         }
 485       }
 486     }
 487     fGeomMatrixSet=kTRUE;
 488   }
 489
 490   // setup digit array if needed
 491   if (!fDigitsArr) {
 492     fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
 493     fDigitsArr->SetOwner(1);
 494   }
 495
 496   // then setup clusterizer
 497   delete fClusterizer;
 498   if     (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
 499     fClusterizer = new AliEMCALClusterizerv1(geometry);
 500   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
 501     fClusterizer = new AliEMCALClusterizerNxN(geometry);
 502   else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
 503     fClusterizer = new AliEMCALClusterizerv2(geometry);
 504   else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
 505    AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
 506    clusterizer->SetNRowDiff(2);
 507    clusterizer->SetNColDiff(2);
 508    fClusterizer = clusterizer;
 509   } else{
 510     AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
 511   }
 512   fClusterizer->InitParameters(fRecParam);
 513
 514   if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
 515     AliCDBManager *cdb = AliCDBManager::Instance();
 516     if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
 517       cdb->SetDefaultStorage(fOCDBpath);
 518     if (fRun!=cdb->GetRun())
 519       cdb->SetRun(fRun);
 520   }
 521   if (!fCalibData&&fLoadCalib) {
 522     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
 523     if (entry)
 524       fCalibData =  static_cast<AliEMCALCalibData*>(entry->GetObject());
 525     if (!fCalibData)
 526       AliFatal("Calibration parameters not found in CDB!");
 527   }
 528   if (!fPedestalData&&fLoadPed) {
 529     AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
 530     if (entry)
 531       fPedestalData =  static_cast<AliCaloCalibPedestal*>(entry->GetObject());
 532   }
 533   if (fCalibData) {
 534     fClusterizer->SetInputCalibrated(kFALSE);
 535     fClusterizer->SetCalibrationParameters(fCalibData);
 536   } else {
 537     fClusterizer->SetInputCalibrated(kTRUE);
 538   }
 539   fClusterizer->SetCaloCalibPedestal(fPedestalData);
 540   fClusterizer->SetJustClusters(kTRUE);
 541   fClusterizer->SetDigitsArr(fDigitsArr);
 542   fClusterizer->SetOutput(0);
 543   fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
 544 }