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 // Fill Unit Array
  18 // Called by ESD reader for jet analysis
  19 // Author: Magali Estienne (magali.estienne@ires.in2p3.fr)
  20
  21 // --- Standard library ---
  22 #include <Riostream.h>
  23 #include <assert.h>
  24
  25 // --- ROOT system ---
  26 #include <TSystem.h>
  27 #include <TLorentzVector.h>
  28 #include <TRefArray.h>
  29 #include "TTask.h"
  30 #include <TGeoManager.h>
  31 #include <TMath.h>
  32 #include <TClonesArray.h>
  33 #include <TArrayS.h>
  34
  35 // --- AliRoot header files ---
  36 #include "AliJetFinder.h"
  37 #include "AliJetReaderHeader.h"
  38 #include "AliJetReader.h"
  39 #include "AliJetESDReader.h"
  40 #include "AliJetESDReaderHeader.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDVertex.h"
  43 #include "AliJetDummyGeo.h"
  44 #include "AliESDCaloCluster.h"
  45 #include "AliESDCaloCells.h"
  46 #include "AliJetUnitArray.h"
  47 #include "AliJetFillUnitArrayEMCalDigits.h"
  48 // Remove CDB dependence under construction
  49 //#include "AliEMCALCalibData.h"
  50 //#include "AliCDBManager.h"
  51
  52 //class AliCDBStorage;
  53 //#include "AliCDBEntry.h"
  54
  55
  56 ClassImp(AliJetFillUnitArrayEMCalDigits)
  57
  58 //_____________________________________________________________________________
  59 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits()
  60   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
  61     fESD(0),
  62     fNumUnits(0),
  63     fEtaMinCal(0.),
  64     fEtaMaxCal(0.),
  65     fPhiMinCal(0.),
  66     fPhiMaxCal(0.),
  67     fNIn(0),
  68     fOpt(0),
  69     fCluster(0),
  70     fDebug(0),
  71     fNCEMCAL(0),
  72     fNCPHOS(0),
  73     fNCCalo(0),
  74     fTPCGrid(0x0),
  75     fEMCalGrid(0x0),
  76     fECorrection(0),
  77     fReaderHeader(0x0),
  78     fMomentumArray(0x0),
  79     fUnitArray(0x0),
  80     fRefArray(0x0),
  81     fProcId(kFALSE),
  82     fGeom(0x0),
  83     fClus(0x0),
  84     fNDigitEmcal(0),
  85     fNDigitEmcalCut(0),
  86     fCalibData(0x0),
  87     fADCchannelECA(0),
  88     fADCpedestalECA(0)
  89 {
  90   // constructor
  91 }
  92
  93 //_____________________________________________________________________________
  94 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(AliESDEvent *esd)
  95   : TTask("AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
  96     fESD(esd),
  97     fNumUnits(0),
  98     fEtaMinCal(0.),
  99     fEtaMaxCal(0.),
 100     fPhiMinCal(0.),
 101     fPhiMaxCal(0.),
 102     fNIn(0),
 103     fOpt(0),
 104     fCluster(0),
 105     fDebug(0),
 106     fNCEMCAL(0),
 107     fNCPHOS(0),
 108     fNCCalo(0),
 109     fTPCGrid(0x0),
 110     fEMCalGrid(0x0),
 111     fECorrection(0),
 112     fReaderHeader(0x0),
 113     fMomentumArray(0x0),
 114     fUnitArray(0x0),
 115     fRefArray(0x0),
 116     fProcId(kFALSE),
 117     fGeom(0x0),
 118     fClus(0x0),
 119     fNDigitEmcal(0),
 120     fNDigitEmcalCut(0),
 121     fCalibData(0x0),
 122     fADCchannelECA(0),
 123     fADCpedestalECA(0)
 124 {
 125   // constructor 2
 126 }
 127
 128 //_____________________________________________________________________________
 129 AliJetFillUnitArrayEMCalDigits::AliJetFillUnitArrayEMCalDigits(const AliJetFillUnitArrayEMCalDigits &det)
 130   : TTask(det),//"AliJetFillUnitArrayEMCalDigits","Fill Unit Array with tpc/its and emcal information"),
 131     fESD(det.fESD),
 132     fNumUnits(det.fNumUnits),
 133     fEtaMinCal(det.fEtaMinCal),
 134     fEtaMaxCal(det.fEtaMaxCal),
 135     fPhiMinCal(det.fPhiMinCal),
 136     fPhiMaxCal(det.fPhiMaxCal),
 137     fNIn(det.fNIn),
 138     fOpt(det.fOpt),
 139     fCluster(det.fCluster),
 140     fDebug(det.fDebug),
 141     fNCEMCAL(det.fNCEMCAL),
 142     fNCPHOS(det.fNCPHOS),
 143     fNCCalo(det.fNCCalo),
 144     fTPCGrid(det.fTPCGrid),
 145     fEMCalGrid(det.fEMCalGrid),
 146     fECorrection(det.fECorrection),
 147     fReaderHeader(det.fReaderHeader),
 148     fMomentumArray(det.fMomentumArray),
 149     fUnitArray(det.fUnitArray),
 150     fRefArray(det.fRefArray),
 151     fProcId(det.fProcId),
 152     fGeom(det.fGeom),
 153     fClus(det.fClus),
 154     fNDigitEmcal(det.fNDigitEmcal),
 155     fNDigitEmcalCut(det.fNDigitEmcalCut),
 156     fCalibData(det.fCalibData),
 157     fADCchannelECA(det.fADCchannelECA),
 158     fADCpedestalECA(det.fADCpedestalECA)
 159 {
 160   // Copy constructor
 161 }
 162
 163 //_____________________________________________________________________________
 164 AliJetFillUnitArrayEMCalDigits& AliJetFillUnitArrayEMCalDigits::operator=(const AliJetFillUnitArrayEMCalDigits& other)
 165 {
 166   // Assignment
 167
 168   fESD = other.fESD;
 169   fNumUnits = other.fNumUnits;
 170   fEtaMinCal = other.fEtaMinCal;
 171   fEtaMaxCal = other.fEtaMaxCal;
 172   fPhiMinCal = other.fPhiMinCal;
 173   fPhiMaxCal = other.fPhiMaxCal;
 174   fNIn = other.fNIn;
 175   fOpt = other.fOpt;
 176   fCluster = other.fCluster;
 177   fDebug = other.fDebug;
 178   fNCEMCAL = other.fNCEMCAL;
 179   fNCPHOS = other.fNCPHOS;
 180   fNCCalo = other.fNCCalo;
 181   fTPCGrid = other.fTPCGrid;
 182   fEMCalGrid = other.fEMCalGrid;
 183   fECorrection = other.fECorrection;
 184   fReaderHeader = other.fReaderHeader;
 185   fMomentumArray = other.fMomentumArray;
 186   fUnitArray = other.fUnitArray;
 187   fRefArray = other.fRefArray;
 188   fProcId = other.fProcId;
 189   fGeom = other.fGeom;
 190   fClus = other.fClus;
 191   fNDigitEmcal = other.fNDigitEmcal;
 192   fNDigitEmcalCut = other.fNDigitEmcalCut;
 193   fCalibData = other.fCalibData;
 194   fADCchannelECA = other.fADCchannelECA;
 195   fADCpedestalECA = other.fADCpedestalECA;
 196
 197   return (*this);
 198
 199 }
 200
 201 //____________________________________________________________________________
 202 void AliJetFillUnitArrayEMCalDigits::InitParameters()
 203 {
 204   fNumUnits = fGeom->GetNCells();      // Number of towers in EMCAL
 205
 206   fEtaMinCal = fGeom->GetArm1EtaMin();
 207   fEtaMaxCal = fGeom->GetArm1EtaMax();
 208   fPhiMinCal = fGeom->GetArm1PhiMin();
 209   fPhiMaxCal = fGeom->GetArm1PhiMax();
 210   fClus      = 0;
 211
 212   // Get calibration parameters from file or digitizer default values.
 213   // Under construction
 214   //  GetCalibrationParameters() ;
 215
 216   if(fDebug>1) printf("\n EMCAL parameters initiated ! \n");
 217
 218 }
 219
 220 //_____________________________________________________________________________
 221 AliJetFillUnitArrayEMCalDigits::~AliJetFillUnitArrayEMCalDigits()
 222 {
 223   // destructor
 224 }
 225
 226 //_____________________________________________________________________________
 227 void AliJetFillUnitArrayEMCalDigits::Exec(Option_t* /*option*/)
 228 {
 229   //
 230   // Main method.
 231   //
 232
 233   fDebug   = fReaderHeader->GetDebug();
 234   fOpt     = fReaderHeader->GetDetector();
 235   fCluster = fReaderHeader->GetCluster();
 236
 237   // Init parameters
 238   InitParameters();
 239
 240   Int_t   goodDigit      = 0;
 241   Int_t   index          = 0;
 242
 243   if(!fCluster) { // Keep all digit information
 244     // Loop over all cell information
 245     //------------------------------------------------------------------
 246     AliESDCaloCells &cells = *(fESD->GetEMCALCells());
 247     Int_t            ncell = cells.GetNumberOfCells() ;
 248
 249     for (Int_t icell=  0; icell <  ncell; icell++) {
 250       Int_t     digitID   = cells.GetCellNumber(icell);
 251       Double_t  digitAmp  = cells.GetAmplitude(icell);
 252       Float_t digitEn   = digitAmp*0.0153; // Last correct
 253       //      Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
 254
 255       Float_t etaD=-10., phiD=-10.;
 256       fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 257       //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 258
 259       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 260
 261       Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 262
 263       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 264
 265       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 266       uArray->SetUnitTrackID(digitID);
 267
 268       Float_t unitEnergy = 0.;
 269       Bool_t ok = kFALSE;
 270       unitEnergy = uArray->GetUnitEnergy();
 271
 272       if(unitEnergy==0){
 273         if(!fProcId){
 274           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 275           fProcId = kTRUE;
 276         }
 277         fRefArray->Add(uArray);
 278         fNDigitEmcal++;
 279         ok = kTRUE;
 280       }
 281
 282       // Detector flag
 283       if(unitEnergy>0.)
 284         uArray->SetUnitDetectorFlag(kAll);
 285       else uArray->SetUnitDetectorFlag(kEmcal);
 286
 287       uArray->SetUnitEnergy(unitEnergy+digitEt);
 288
 289       uArray->SetUnitCutFlag(kPtHigher);
 290
 291       // To be modified !!!
 292       uArray->SetUnitSignalFlag(kGood);
 293
 294       // This is for jet multiplicity
 295       uArray->SetUnitClusterID(index);
 296
 297       if(fDebug > 1) printf("goodDigit : %d\n", goodDigit);
 298
 299     } // End loop over cells
 300
 301   } // end if !fCluster
 302   else { // Keep digit information from clusterization
 303
 304     // Loop over calo clusters
 305     //------------------------------------------------------------------
 306
 307     //select EMCAL clusters only
 308     TRefArray * caloClusters  = new TRefArray();
 309     fESD->GetEMCALClusters(caloClusters);
 310
 311     // Total number of EMCAL cluster
 312     Int_t nclus = caloClusters->GetEntries() ;
 313     Int_t beg   = 0;
 314     Float_t pos[3] ;
 315
 316     // Get reconstructed vertex position
 317     Double_t vertex_position[3] ;
 318     fESD->GetVertex()->GetXYZ(vertex_position) ;
 319
 320     // Get CaloCells
 321     AliESDCaloCells &cells= *(fESD->GetEMCALCells());
 322
 323     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
 324       // Retrieve cluster from esd
 325       fClus = (AliESDCaloCluster *) caloClusters->At(j) ;
 326
 327       // Get the cluster info
 328       Float_t energy       = fClus->E() ;
 329       Int_t iprim          = fClus->GetLabel();
 330       Int_t trackIndex     = fClus->GetTrackMatched();
 331
 332       fClus->GetPosition(pos) ;
 333       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 334       TLorentzVector p ;
 335       fClus->GetMomentum(p,vertex_position);
 336
 337       Int_t     digMult = fClus->GetNCells() ;
 338       UShort_t *digID   = fClus->GetCellsAbsId() ;
 339       //Print cluster info
 340       if(fDebug>2) cout<<"Cluster "<< j <<"; digits mult "<<digMult<<"; type "<<(Int_t )fClus->GetClusterType()
 341                        <<"; Energy "<<energy<< "; transverse energy:" <<  energy*TMath::Abs(TMath::Sin(EtaToTheta(vpos.Eta())))
 342                        <<"; Phi "<<vpos.Phi()<<"; Eta "<<vpos.Eta() <<"; label "<<iprim<<endl;
 343
 344       // Do double-counted electron correction
 345       if (fECorrection != 0 && trackIndex !=-1 )
 346         {
 347           // The electron correction go there
 348           // Under construction !!!!
 349
 350         } // End of Electron correction
 351
 352       // Get CaloCells of cluster and fill the unitArray
 353       for(Int_t i = 0; i < digMult ; i++)
 354         {
 355           Int_t     digitID     = digID[i]; // or clus->GetCellNumber(i) ;
 356           Double_t  digitAmp    = cells.GetCellAmplitude(digitID) ;
 357
 358           // Calibration for an energy in GeV
 359           Float_t digitEn = digitAmp*0.0153;
 360           //      Float_t digitEn = Calibrate((Int_t)digitAmp,digitID);
 361
 362           Float_t etaD=-10., phiD=-10.;
 363           fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 364           //          fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 365
 366           phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 367
 368           Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 369
 370           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 371           if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 372           uArray->SetUnitTrackID(digitID);
 373
 374           Float_t unitEnergy = 0.;
 375           Bool_t ok = kFALSE;
 376           unitEnergy = uArray->GetUnitEnergy();
 377
 378           if(unitEnergy==0){
 379             if(!fProcId){
 380               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 381               fProcId = kTRUE;
 382             }
 383             fRefArray->Add(uArray);
 384             fNDigitEmcal++;
 385             ok = kTRUE;
 386           }
 387
 388           // Detector flag
 389           if(unitEnergy>0.)
 390             uArray->SetUnitDetectorFlag(kAll);
 391           else uArray->SetUnitDetectorFlag(kEmcal);
 392
 393           uArray->SetUnitEnergy(unitEnergy+digitEt);
 394           uArray->SetUnitCutFlag(kPtHigher);
 395
 396           // To be modified !!!
 397           uArray->SetUnitSignalFlag(kGood);
 398
 399           // This is for jet multiplicity
 400           uArray->SetUnitClusterID(index);
 401
 402         } // End loop over cells
 403     } // End loop over clusters
 404   } // end else
 405
 406   fNIn += goodDigit;
 407
 408   if(fDebug>1)
 409     {
 410       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
 411       printf("goodDigit : %d\n", goodDigit);
 412     }
 413
 414 }
 415
 416 // //____________________________________________________________________________
 417 // void AliJetFillUnitArrayEMCalDigits::GetCalibrationParameters()
 418 // {
 419 //   // Set calibration parameters:
 420 //   // if calibration database exists, they are read from database,
 421 //   // otherwise, they are taken from digitizer.
 422 //   //
 423 //   // It is a user responsilibity to open CDB before reconstruction,
 424 //   // for example:
 425 //   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 426
 427 //   //Check if calibration is stored in data base
 428
 429 //   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
 430 //     {
 431 //       AliCDBEntry *entry = (AliCDBEntry*)
 432 //         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 433 //       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 434 //     }
 435
 436 //   if(!fCalibData)
 437 //     printf("************* Calibration parameters not found in CDB! ****************");
 438 // //    AliFatal("Calibration parameters not found in CDB!");
 439
 440
 441 // }
 442
 443 // //____________________________________________________________________________
 444 // Float_t  AliJetFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
 445 // {
 446
 447 //   // Convert digitized amplitude into energy.
 448 //   // Calibration parameters are taken from calibration data base for raw data,
 449 //   // or from digitizer parameters for simulated data.
 450
 451 //   if(fCalibData){
 452
 453 //     if (fGeom==0)
 454 //       printf("************* Did not get geometry from EMCALLoader ***************");
 455
 456 //     Int_t iSupMod  = -1;
 457 //     Int_t nModule  = -1;
 458 //     Int_t nIphi    = -1;
 459 //     Int_t nIeta    = -1;
 460 //     Int_t iphi     = -1;
 461 //     Int_t ieta     = -1;
 462
 463 //     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 464 //     if(!bCell) {
 465 //       // fGeom->PrintGeometry();
 466 //       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
 467 //       assert(0);
 468 //     }
 469
 470 //     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 471
 472 //     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 473 //     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 474
 475 //    return -fADCpedestalECA + amp * fADCchannelECA ;
 476
 477 //   }
 478 //   else //Return energy with default parameters if calibration is not available
 479 //     return -fADCpedestalECA + amp * fADCchannelECA ;
 480
 481 // }
 482
 483
 484 //_____________________________________________________________________________
 485 Float_t  AliJetFillUnitArrayEMCalDigits::EtaToTheta(Float_t arg)
 486 {
 487   //  return (180./TMath::Pi())*2.*atan(exp(-arg));
 488   return 2.*atan(exp(-arg));
 489
 490
 491 }