JETAN/AliJetFillUnitArrayEMCalDigits.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 with EMCal information
  19 // Called by ESD reader for jet analysis
  20 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  21 //-------------------------------------------------------------
  22
  23 // --- Standard library ---
  24 #include <Riostream.h>
  25
  26 // --- ROOT system ---
  27 #include <TSystem.h>
  28 #include <TLorentzVector.h>
  29 #include <TRefArray.h>
  30 #include "TTask.h"
  31 #include <TGeoManager.h>
  32 #include <TMath.h>
  33 #include <TClonesArray.h>
  34 #include <TArrayS.h>
  35
  36 // --- AliRoot header files ---
  37 #include "AliJetFinder.h"
  38 #include "AliJetReaderHeader.h"
  39 #include "AliJetReader.h"
  40 #include "AliJetESDReader.h"
  41 #include "AliJetESDReaderHeader.h"
  42 //#include "AliESD.h"
  43 #include "AliESDEvent.h"
  44 #include "AliJetDummyGeo.h"
  45 #include "AliESDCaloCluster.h"
  46 #include "AliJetUnitArray.h"
  47 #include "AliJetFillUnitArrayEMCalDigits.h"
  48 #include "AliCDBManager.h"
  49 #include "AliCDBStorage.h"
  50 #include "AliCDBEntry.h"
  51
  52 ClassImp(AliJetFillUnitArrayEMCalDigits)
  53
  54 //_____________________________________________________________________________
  55 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
  56   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
  57     fESD(0),
  58     fNumUnits(0),
  59     fEtaMinCal(0.),
  60     fEtaMaxCal(0.),
  61     fPhiMinCal(0.),
  62     fPhiMaxCal(0.),
  63     fNIn(0),
  64     fOpt(0),
  65     fDebug(0),
  66     fNCEMCAL(0),
  67     fNCPHOS(0),
  68     fNCCalo(0),
  69     fTPCGrid(0x0),
  70     fEMCalGrid(0x0),
  71     fReaderHeader(0x0),
  72     fMomentumArray(0x0),
  73     fUnitArray(0x0),
  74     fRefArray(0x0),
  75     fGeom(0x0),
  76     fClus(0x0),
  77     fNDigitEmcal(0),
  78     fNDigitEmcalCut(0)
  79 {
  80   // constructor
  81 }
  82
  83 //_____________________________________________________________________________
  84 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent */*esd*/)
  85   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
  86     fESD(0),
  87     fNumUnits(0),
  88     fEtaMinCal(0.),
  89     fEtaMaxCal(0.),
  90     fPhiMinCal(0.),
  91     fPhiMaxCal(0.),
  92     fNIn(0),
  93     fOpt(0),
  94     fDebug(0),
  95     fNCEMCAL(0),
  96     fNCPHOS(0),
  97     fNCCalo(0),
  98     fTPCGrid(0x0),
  99     fEMCalGrid(0x0),
 100     fReaderHeader(0x0),
 101     fMomentumArray(0x0),
 102     fUnitArray(0x0),
 103     fRefArray(0x0),
 104     fGeom(0x0),
 105     fClus(0x0),
 106     fNDigitEmcal(0),
 107     fNDigitEmcalCut(0)
 108 {
 109   // constructor
 110 }
 111
 112
 113 //____________________________________________________________________________
 114 void AliJetFillUnitArrayEMCalDigits::InitParameters()
 115 {
 116   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
 117
 118   fEtaMinCal = fGeom->GetArm1EtaMin();
 119   fEtaMaxCal = fGeom->GetArm1EtaMax();
 120   fPhiMinCal = fGeom->GetArm1PhiMin();
 121   fPhiMaxCal = fGeom->GetArm1PhiMax();
 122   fClus      = 0;
 123
 124   if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
 125
 126 }
 127
 128 //_____________________________________________________________________________
 129 AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
 130 {
 131   // destructor
 132 }
 133
 134 //_____________________________________________________________________________
 135 void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
 136 {
 137   //
 138   // Main method.
 139   // Explain
 140
 141   fDebug = fReaderHeader->GetDebug();
 142   fOpt = fReaderHeader->GetDetector();
 143
 144   // Init parameters
 145   InitParameters();
 146
 147   // Get number of clusters from EMCAL
 148
 149   Int_t   nDigitTot      = 0;
 150   Int_t   goodDigit      = 0;
 151   Int_t   beg            = 0;
 152   Int_t   end            = 0;
 153   Float_t ptMin          = fReaderHeader->GetPtCut();
 154
 155   // Loop over calo clusters
 156   //------------------------------------------------------------------
 157   Int_t type = 0;
 158   Int_t index = 0;
 159
 160   // Total number of EMCAL cluster
 161   end =  fESD->GetNumberOfCaloClusters();
 162
 163   for(Int_t j = beg; j < end; j++) {
 164       fClus = fESD->GetCaloCluster(j);
 165       if(!fClus->IsEMCAL()) continue;
 166
 167       type = fClus->GetClusterType();
 168       index = fClus->GetID();
 169       nDigitTot = fClus->GetNumberOfDigits();
 170
 171       // Keep clusters or pseudo clusters
 172       if (type != AliESDCaloCluster::kClusterv1) continue;
 173       //      if (type != AliESDCaloCluster::kPseudoCluster) continue;
 174
 175       // Get the digit index and the digit information
 176       //============================================================
 177
 178       // Get number of digits in a cluster
 179       Int_t nD = fClus->GetNumberOfDigits();
 180
 181       TArrayS *digID = fClus->GetDigitIndex();
 182       TArrayS *digEnergy = fClus->GetDigitAmplitude();
 183       Float_t *digitEnergy = new Float_t[nD];
 184       //      Float_t digitEn = 0.;
 185
 186       // Loop over digits
 187       for(Int_t k=0; k<nD; k++) {
 188
 189         // Convert energy in GeV
 190         Int_t idF = (Int_t)digID->At(k);
 191         // Calibration for an energy in GeV
 192         digitEnergy[k] = (Float_t)digEnergy->At(k)/500.;
 193
 194         // Second method to extract eta, phi positions of a digit
 195         //=================================================================
 196
 197         Float_t etaD=-10., phiD=-10.;
 198         fGeom->EtaPhiFromIndex(idF,etaD,phiD);
 199         //          fEMCalGrid->GetEtaPhiFromIndex2(idF,phiD,etaD);
 200         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 201
 202         Float_t etDigit = digitEnergy[k]*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 203
 204         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idF);
 205         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 206
 207         Float_t unitEnergy = 0.;
 208         Bool_t ok = kFALSE;
 209         unitEnergy = uArray->GetUnitEnergy();
 210         if(unitEnergy==0){
 211           fRefArray->Add(uArray);
 212           fNDigitEmcal++;
 213           ok = kTRUE;
 214         }
 215         uArray->SetUnitEnergy(unitEnergy+etDigit);
 216         // Put a cut flag
 217         if(uArray->GetUnitEnergy()<ptMin)
 218           uArray->SetUnitCutFlag(kPtSmaller);
 219         else {
 220           uArray->SetUnitCutFlag(kPtHigher);
 221           if(ok) fNDigitEmcalCut++;
 222         }
 223         // Detector flag
 224         if(unitEnergy>0.)
 225           uArray->SetUnitDetectorFlag(kAll);
 226         else uArray->SetUnitDetectorFlag(kEmcal);
 227
 228         // This is for jet multiplicity
 229         uArray->SetUnitClusterID(index);
 230
 231         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
 232
 233       } // End loop over digits
 234
 235   } // End loop over clusters
 236
 237   fNIn += goodDigit;
 238
 239 }
 240
 241 //_____________________________________________________________________________
 242 Float_t  AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
 243 {
 244   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 245   return 2.*atan(exp(-arg));
 246
 247
 248 }