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