JETAN/AliJetESDFillUnitArrayEMCalDigits.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 //---------------------------------------------------------------------
  18 // Fill Unit Array
  19 // Called by ESD reader for jet analysis
  20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  21 //---------------------------------------------------------------------
  22
  23 // --- ROOT system ---
  24 #include <TLorentzVector.h>
  25 #include <TMath.h>
  26
  27 // --- AliRoot header files ---
  28 #include "AliESDCaloCluster.h"
  29 #include "AliESDCaloCells.h"
  30 #include "AliJetUnitArray.h"
  31 #include "AliJetESDFillUnitArrayEMCalDigits.h"
  32 // #include "AliEMCALCalibData.h"
  33 // #include "AliCDBManager.h"
  34 // // class AliCDBStorage;
  35 // #include "AliCDBEntry.h"
  36
  37 // --- ROOT system ---
  38 class TSystem;
  39 class TGeoManager;
  40 class TArrayS;
  41
  42 // --- AliRoot header files ---
  43 class AliJetFinder;
  44 class AliJetReader;
  45 class AliJetESDReader;
  46 // class AliEMCALCalibData;
  47 // class AliCDBManager;
  48 // class AliCDBStorage;
  49 // class AliCDBEntry;
  50
  51 using std::cout;
  52 using std::endl;
  53 ClassImp(AliJetESDFillUnitArrayEMCalDigits)
  54
  55 //_____________________________________________________________________________
  56 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits()
  57   : AliJetFillUnitArray(),
  58     fESD(0),
  59     fNIn(0),
  60     fOpt(0),
  61     fCluster(0),
  62     fDebug(0),
  63     fNCEMCAL(0),
  64     fNCPHOS(0),
  65     fNCCalo(0),
  66     fApplyElectronCorrection(kFALSE),
  67     fApplyFractionHadronicCorrection(kFALSE),
  68     fFractionHadronicCorrection(0.3),
  69     fClus(0x0),
  70     fNDigitEmcal(0),
  71     fNDigitEmcalCut(0)
  72 {
  73   // constructor
  74 }
  75
  76 //_____________________________________________________________________________
  77 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(AliESDEvent *esd)
  78   : AliJetFillUnitArray(),
  79     fESD(esd),
  80     fNIn(0),
  81     fOpt(0),
  82     fCluster(0),
  83     fDebug(0),
  84     fNCEMCAL(0),
  85     fNCPHOS(0),
  86     fNCCalo(0),
  87     fApplyElectronCorrection(kFALSE),
  88     fApplyFractionHadronicCorrection(kFALSE),
  89     fFractionHadronicCorrection(0.3),
  90     fClus(0x0),
  91     fNDigitEmcal(0),
  92     fNDigitEmcalCut(0)
  93 {
  94   // constructor
  95 }
  96
  97 //_____________________________________________________________________________
  98 AliJetESDFillUnitArrayEMCalDigits::AliJetESDFillUnitArrayEMCalDigits(const AliJetESDFillUnitArrayEMCalDigits &det)
  99   : AliJetFillUnitArray(det),
 100     fESD(det.fESD),
 101     fNIn(det.fNIn),
 102     fOpt(det.fOpt),
 103     fCluster(det.fCluster),
 104     fDebug(det.fDebug),
 105     fNCEMCAL(det.fNCEMCAL),
 106     fNCPHOS(det.fNCPHOS),
 107     fNCCalo(det.fNCCalo),
 108     fApplyElectronCorrection(det.fApplyElectronCorrection),
 109     fApplyFractionHadronicCorrection(det.fApplyFractionHadronicCorrection),
 110     fFractionHadronicCorrection(det.fFractionHadronicCorrection),
 111     fClus(det.fClus),
 112     fNDigitEmcal(det.fNDigitEmcal),
 113     fNDigitEmcalCut(det.fNDigitEmcalCut)
 114 {
 115   // Copy constructor
 116 }
 117
 118 //_____________________________________________________________________________
 119 AliJetESDFillUnitArrayEMCalDigits& AliJetESDFillUnitArrayEMCalDigits::operator=(const AliJetESDFillUnitArrayEMCalDigits& other)
 120 {
 121   // Assignment
 122
 123   fESD = other.fESD;
 124   fNIn = other.fNIn;
 125   fOpt = other.fOpt;
 126   fCluster = other.fCluster;
 127   fDebug = other.fDebug;
 128   fNCEMCAL = other.fNCEMCAL;
 129   fNCPHOS = other.fNCPHOS;
 130   fNCCalo = other.fNCCalo;
 131   fApplyElectronCorrection = other.fApplyElectronCorrection;
 132   fApplyFractionHadronicCorrection = other.fApplyFractionHadronicCorrection;
 133   fFractionHadronicCorrection = other.fFractionHadronicCorrection;
 134   fClus = other.fClus;
 135   fNDigitEmcal = other.fNDigitEmcal;
 136   fNDigitEmcalCut = other.fNDigitEmcalCut;
 137
 138   return (*this);
 139
 140 }
 141
 142 //_____________________________________________________________________________
 143 AliJetESDFillUnitArrayEMCalDigits::~AliJetESDFillUnitArrayEMCalDigits()
 144 {
 145   // destructor
 146 }
 147
 148 //_____________________________________________________________________________
 149 void AliJetESDFillUnitArrayEMCalDigits::Exec(Option_t* const /*option*/)
 150 {
 151   //
 152   // Main method.
 153   // Fill the unit array with the neutral particle information from the EMCal cells in ESD
 154   //
 155
 156   fDebug = fReaderHeader->GetDebug();
 157   fOpt = fReaderHeader->GetDetector();
 158   fCluster = fReaderHeader->GetCluster();
 159
 160   // Init parameters
 161   InitParameters();
 162
 163 //(not used ?)  Int_t   nDigitTot      = 0;
 164   Int_t   goodDigit      = 0;
 165 //(not used ?)  Int_t   beg            = 0;
 166 //(not used ?)  Int_t   end            = 0;
 167   Int_t   index          = 0;
 168 //(not used ?)  Int_t   count          = 0;
 169
 170 //(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
 171
 172   if(!fCluster) { // Keep all digit information
 173     // Loop over all cell information
 174     //------------------------------------------------------------------
 175     AliESDCaloCells &cells= *(fESD->GetEMCALCells());
 176     Int_t ncell = cells.GetNumberOfCells() ;
 177 //(not used ?)    Int_t type = cells.GetType();
 178
 179     for (Int_t icell=  0; icell <  ncell; icell++) {
 180       Int_t   digitID   = cells.GetCellNumber(icell);
 181       Float_t digitAmp  = cells.GetAmplitude(icell);
 182 //(not used ?)      Float_t digitTime = cells.GetTime(icell);
 183       Float_t digitEn   = digitAmp*0.0153; // Last correct
 184       //  Float_t digitEn = Calibrate(digitAmp,digitID);
 185
 186       Float_t etaD=-10., phiD=-10.;
 187       fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 188       //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 189
 190       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 191
 192       Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 193
 194       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 195
 196       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 197       uArray->SetUnitTrackID(digitID);
 198
 199       Float_t unitEnergy = 0.;
 200       Bool_t ok = kFALSE;
 201       unitEnergy = uArray->GetUnitEnergy();
 202
 203       if(unitEnergy==0){
 204         if(!fProcId){
 205           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 206           fProcId = kTRUE;
 207         }
 208         fRefArray->Add(uArray);
 209         fNDigitEmcal++;
 210         ok = kTRUE;
 211       }
 212
 213       // Detector flag
 214       if(unitEnergy>0.)
 215         uArray->SetUnitDetectorFlag(kAll);
 216       else uArray->SetUnitDetectorFlag(kEmcal);
 217
 218       uArray->SetUnitEnergy(unitEnergy+digitEt);
 219
 220       uArray->SetUnitCutFlag(kPtHigher);
 221
 222       // To be modified !!!
 223       uArray->SetUnitSignalFlag(kGood);
 224
 225       // This is for jet multiplicity
 226       uArray->SetUnitClusterID(index);
 227
 228       if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
 229
 230     } // End loop over cells
 231   } // end if !fCluster
 232   else { // Keep digit information from clusterization
 233
 234     // Loop over calo clusters
 235     //------------------------------------------------------------------
 236
 237     //select EMCAL clusters only
 238     TRefArray * caloClusters  = new TRefArray();
 239     fESD->GetEMCALClusters(caloClusters);
 240
 241     // Total number of EMCAL cluster
 242     Int_t nclus = caloClusters->GetEntries() ;
 243     Int_t beg   = 0;
 244     Float_t pos[3] ;
 245
 246     // Get reconstructed vertex position
 247     Double_t vertexPosition[3] ;
 248     fESD->GetVertex()->GetXYZ(vertexPosition) ;
 249
 250     // Get CaloCells
 251     AliESDCaloCells &cells= *(fESD->GetEMCALCells());
 252 //(not used ?)    Int_t ncell = cells.GetNumberOfCells() ;
 253 //(not used ?)    Int_t type = cells.GetType();
 254
 255     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
 256       // Retrieve cluster from esd
 257       fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
 258
 259       // Get the cluster info
 260 //(not used ?)      Float_t energy         = fClus->E() ;
 261 //(not used ?)      Int_t   iprim          = fClus->GetLabel();
 262
 263       fClus->GetPosition(pos) ;
 264       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 265       TLorentzVector p ;
 266       fClus->GetMomentum(p,vertexPosition);
 267
 268       Int_t     digMult = fClus->GetNCells() ;
 269       UShort_t *digID   = fClus->GetCellsAbsId() ;
 270 //(not used ?)      Double_t *digAmpFrac = fClus->GetCellsAmplitudeFraction() ;
 271       Int_t     trackIndex = fClus->GetTrackMatchedIndex();
 272
 273       // Do double-counted electron correction
 274       if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
 275         {
 276           // The electron correction go there
 277           // Under construction !!!!
 278         }  // End of Electron correction
 279
 280       // Get CaloCells of cluster and fill the unitArray
 281       for(Int_t i = 0; i < digMult ; i++){
 282         Int_t    digitID   = digID[i]; // or clus->GetCellNumber(i) ;
 283 //(not used ?)        Double_t digitAmpFrac = digAmpFrac[i];
 284         Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
 285 //(not used ?)        Float_t  digitTime    = cells.GetCellTime(digitID);
 286
 287         // Calibration for an energy in GeV
 288         Float_t digitEn = digitAmp*0.0153;
 289
 290         Float_t etaD=-10., phiD=-10.;
 291         fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 292         //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 293
 294         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 295
 296         Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 297
 298         cout << "Digit " << i << ", eta: " << etaD << ", phi: " << phiD << endl;
 299
 300         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 301         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 302         uArray->SetUnitTrackID(digitID);
 303
 304         Float_t unitEnergy = 0.;
 305         Bool_t ok = kFALSE;
 306         unitEnergy = uArray->GetUnitEnergy();
 307
 308         if(unitEnergy==0){
 309           if(!fProcId){
 310             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 311             fProcId = kTRUE;
 312           }
 313           fRefArray->Add(uArray);
 314           fNDigitEmcal++;
 315           ok = kTRUE;
 316         }
 317
 318         // Detector flag
 319         if(unitEnergy>0.)
 320           uArray->SetUnitDetectorFlag(kAll);
 321         else uArray->SetUnitDetectorFlag(kEmcal);
 322
 323         uArray->SetUnitEnergy(unitEnergy+digitEt);
 324
 325         uArray->SetUnitCutFlag(kPtHigher);
 326
 327         // To be modified !!!
 328         uArray->SetUnitSignalFlag(kGood);
 329
 330         // This is for jet multiplicity
 331         uArray->SetUnitClusterID(index);
 332
 333         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
 334
 335       } // End loop over cells
 336     } // End loop over clusters
 337   } // end else
 338
 339   fNIn += goodDigit;
 340
 341   if(fDebug>1)
 342     {
 343       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
 344       printf("goodDigit : %d\n", goodDigit);
 345     }
 346 }
 347
 348 /*
 349 //____________________________________________________________________________
 350 void AliJetESDFillUnitArrayEMCalDigits::GetCalibrationParameters()
 351 {
 352   // Set calibration parameters:
 353   // if calibration database exists, they are read from database,
 354   // otherwise, they are taken from digitizer.
 355   //
 356   // It is a user responsilibity to open CDB before reconstruction,
 357   // for example:
 358   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 359
 360   //Check if calibration is stored in data base
 361
 362   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
 363     {
 364       AliCDBEntry *entry = (AliCDBEntry*)
 365         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 366       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 367     }
 368
 369   if(!fCalibData)
 370     printf("************* Calibration parameters not found in CDB! ****************");
 371 //    AliFatal("Calibration parameters not found in CDB!");
 372
 373
 374 }
 375
 376 //____________________________________________________________________________
 377 Float_t  AliJetESDFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
 378 {
 379
 380   // Convert digitized amplitude into energy.
 381   // Calibration parameters are taken from calibration data base for raw data,
 382   // or from digitizer parameters for simulated data.
 383
 384   if(fCalibData){
 385
 386     if (fGeom==0)
 387       printf("************* Did not get geometry from EMCALLoader ***************");
 388 //      AliFatal("Did not get geometry from EMCALLoader") ;
 389
 390     Int_t iSupMod = -1;
 391     Int_t nModule  = -1;
 392     Int_t nIphi   = -1;
 393     Int_t nIeta   = -1;
 394     Int_t iphi    = -1;
 395     Int_t ieta    = -1;
 396
 397     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 398     if(!bCell) {
 399       // fGeom->PrintGeometry();
 400       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
 401       assert(0);
 402     }
 403
 404     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 405
 406     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 407     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 408
 409    return -fADCpedestalECA + amp * fADCchannelECA ;
 410
 411   }
 412   else //Return energy with default parameters if calibration is not available
 413     return -fADCpedestalECA + amp * fADCchannelECA ;
 414
 415 }
 416 */
 417