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