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   Int_t   goodDigit      = 0;
 144   Int_t   index          = 0;
 145
 146   if(!fCluster) { // Keep all digit information
 147
 148     // Loop over all cell information
 149     //------------------------------------------------------------------
 150     AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
 151     Int_t ncell = cells.GetNumberOfCells() ;
 152 //(not used ?)    Int_t type = cells.GetType();
 153
 154     for (Int_t icell=  0; icell <  ncell; icell++) {
 155       Int_t   digitID   = cells.GetCellNumber(icell);
 156       Float_t digitAmp  = cells.GetAmplitude(icell);
 157       Float_t digitEn   = digitAmp*0.0153; // Last correct
 158       //  Float_t digitEn = Calibrate(digitAmp,digitID);
 159
 160       Float_t etaD=-10., phiD=-10.;
 161       fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 162       //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 163
 164       phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 165
 166       //      Float_t digitEt = digitEn*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 167
 168       AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 169
 170       if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 171       uArray->SetUnitTrackID(digitID);
 172
 173       Float_t unitEnergy = 0.;
 174       Bool_t ok = kFALSE;
 175       unitEnergy = uArray->GetUnitEnergy();
 176
 177       if(unitEnergy==0){
 178         if(!fProcId){
 179           new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 180           fProcId = kTRUE;
 181         }
 182         fRefArray->Add(uArray);
 183         fNDigitEmcal++;
 184         ok = kTRUE;
 185       }
 186
 187       // Hadronic Correction
 188       double correction=0;
 189       if (fApplyFractionHadronicCorrection) {
 190
 191         TArrayS* matched = new TArrayS();
 192         GetTracksPointingToCell(matched, etaD, phiD,0.02);
 193
 194         double ptot = 0;
 195         for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
 196           {
 197             const short indexS = matched->At(itrack);
 198             if (indexS>-1)
 199               {
 200                 AliAODTrack *track = fAOD->GetTrack(indexS);
 201                 ptot += track->P();
 202               }
 203           }
 204
 205         correction = ptot * fFractionHadronicCorrection;
 206
 207         if (ptot>0 && fDebug>1 ) {
 208            printf("SUM of track momentum for this cell: p=%f \n", ptot);
 209            printf("fractional correction=%f \n", fFractionHadronicCorrection);
 210            printf("TOTAL correction=%f \n", correction );
 211         }
 212
 213         delete matched;
 214
 215       }//end hadronic correction
 216
 217       double enerCorr = digitEn;
 218       if (correction >= enerCorr) enerCorr = 0;
 219       else enerCorr -= correction;
 220       if (correction>0 && fDebug>1) {
 221         printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
 222       }
 223
 224       Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 225
 226       // Detector flag
 227       if(unitEnergy>0.)
 228         uArray->SetUnitDetectorFlag(kAll);
 229       else uArray->SetUnitDetectorFlag(kEmcal);
 230
 231       uArray->SetUnitEnergy(unitEnergy+digitEt);
 232
 233       uArray->SetUnitCutFlag(kPtHigher);
 234
 235       // To be modified !!!
 236       uArray->SetUnitSignalFlag(kGood);
 237
 238       // This is for jet multiplicity
 239       uArray->SetUnitClusterID(index);
 240
 241       if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
 242
 243     } // End loop over cells
 244   } // end if !fCluster
 245   else { // Keep digit information from clusterization
 246
 247     // Loop over calo clusters
 248     //------------------------------------------------------------------
 249
 250     // select EMCAL clusters only
 251     TRefArray * caloClusters  = new TRefArray();
 252     fAOD->GetEMCALClusters(caloClusters);
 253
 254     // Total number of EMCAL cluster
 255     Int_t nclus = caloClusters->GetEntries() ;
 256     Int_t beg   = 0;
 257     Float_t pos[3] ;
 258
 259     // Get CaloCells
 260     AliAODCaloCells &cells= *(fAOD->GetEMCALCells());
 261
 262     for(Int_t j = beg; j < nclus; j++) { // loop over clusters
 263       // Retrieve cluster from aod
 264       fClus = (AliAODCaloCluster *) caloClusters->At(j) ;
 265
 266       // Get the cluster info
 267
 268       fClus->GetPosition(pos) ;
 269       TVector3 vpos(pos[0],pos[1],pos[2]) ;
 270
 271       Int_t     digMult = fClus->GetNCells() ;
 272       UShort_t *digID   = fClus->GetCellsAbsId() ;
 273       Int_t     trackIndex = -1;
 274       if(fClus->GetNTracksMatched()!=0)
 275         trackIndex = ((AliAODTrack*)fClus->GetTrackMatched(0))->GetID();
 276
 277       // Do double-counted electron correction
 278       if (fApplyElectronCorrection != 0 && trackIndex !=-1 )
 279         {
 280           // The electron correction go there
 281           // Under construction !!!!
 282         }  // End of Electron correction
 283
 284       // Get CaloCells of cluster and fill the unitArray
 285       for(Int_t i = 0; i < digMult ; i++) {
 286         Int_t    digitID      = digID[i]; // or clus->GetCellNumber(i) ;
 287         Float_t  digitAmp     = cells.GetCellAmplitude(digitID) ;
 288
 289         // Calibration for an energy in GeV
 290         Float_t digitEn = digitAmp*0.0153;
 291
 292         Float_t etaD=-10., phiD=-10.;
 293         fGeom->EtaPhiFromIndex(digitID,etaD,phiD);
 294         //  fEMCalGrid->GetEtaPhiFromIndex2(digitID,phiD,etaD);
 295
 296         phiD = ((phiD < 0) ? phiD + 2.* TMath::Pi() : phiD);
 297
 298         // Hadronic Correction
 299         double correction=0;
 300         if (fApplyFractionHadronicCorrection) {
 301
 302           TArrayS* matched = new TArrayS();
 303           GetTracksPointingToCell(matched, etaD, phiD,0.02);
 304
 305           double ptot = 0;
 306           for(int itrack = 0; itrack <  matched->GetSize(); itrack++)
 307             {
 308               const short indexS = matched->At(itrack);
 309               if (indexS>-1)
 310                 {
 311                   AliAODTrack *track = fAOD->GetTrack(indexS);
 312                   ptot += track->P();
 313                 }
 314             }
 315
 316           correction = ptot * fFractionHadronicCorrection;
 317
 318           if (ptot>0 && fDebug>1 ) {
 319             printf("SUM of track momentum for this cell: p=%f \n", ptot);
 320             printf("fractional correction=%f \n", fFractionHadronicCorrection);
 321             printf("TOTAL correction=%f \n", correction );
 322           }
 323
 324           delete matched;
 325
 326         }//end hadronic correction
 327
 328         double enerCorr = digitEn;
 329         if (correction >= enerCorr) enerCorr = 0;
 330         else enerCorr -= correction;
 331         if (correction>0 && fDebug>1) {
 332           printf("AliJetAODFillUnitArrayEMCALDigits---Hadronic Correction: uncorrected E=%f, corrected E=%f \n", digitEn, enerCorr);
 333         }
 334
 335         Float_t digitEt = enerCorr*TMath::Abs(TMath::Sin(EtaToTheta(etaD)));
 336
 337         AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(digitID);
 338         if(uArray->GetUnitEnergy() == 0.) goodDigit++;
 339         uArray->SetUnitTrackID(digitID);
 340
 341         Float_t unitEnergy = 0.;
 342         Bool_t ok = kFALSE;
 343         unitEnergy = uArray->GetUnitEnergy();
 344
 345         if(unitEnergy==0){
 346           if(!fProcId){
 347             new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 348             fProcId = kTRUE;
 349           }
 350           fRefArray->Add(uArray);
 351           fNDigitEmcal++;
 352           ok = kTRUE;
 353         }
 354
 355         // Detector flag
 356         if(unitEnergy>0.)
 357           uArray->SetUnitDetectorFlag(kAll);
 358         else uArray->SetUnitDetectorFlag(kEmcal);
 359
 360         uArray->SetUnitEnergy(unitEnergy+digitEt);
 361
 362         uArray->SetUnitCutFlag(kPtHigher);
 363
 364         // Signal or background
 365         // To be modified !!!
 366         uArray->SetUnitSignalFlag(kGood);
 367
 368         // This is for jet multiplicity
 369         uArray->SetUnitClusterID(index);
 370
 371         if(fDebug > 12) printf("goodDigit : %d\n", goodDigit);
 372
 373       }// End loop over cells
 374     } // End loop over clusters
 375   } // end else
 376
 377   fNIn += goodDigit;
 378
 379   if(fDebug>1)
 380     {
 381       printf("End of digits %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%");
 382       printf("goodDigit : %d\n", goodDigit);
 383     }
 384 }
 385
 386 //____________________________________________________________________________
 387 void AliJetAODFillUnitArrayEMCalDigits::GetTracksPointingToCell(TArrayS* array,Double_t eta, Double_t phi, Double_t cut)
 388 {
 389 // Get all tracks pointing to cell
 390   int size=0;
 391
 392   for (Int_t itrk =  0; itrk <  fAOD->GetNumberOfTracks() ; itrk++) { //track loop
 393
 394     AliAODTrack * track = (AliAODTrack*) fAOD->GetTrack(itrk) ;
 395     AliAODPid*    pid   = (AliAODPid*) track->GetDetPid();
 396
 397     if(pid) {
 398       Double_t emcpos[3];
 399       pid->GetEMCALPosition(emcpos);
 400       TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
 401
 402       Double_t deta = tpos.Eta() - eta;
 403       Double_t dphi = tpos.Phi() - phi;
 404
 405       if(dphi > TMath::Pi()) dphi -= 2*TMath::Pi();
 406       if(dphi < -TMath::Pi()) dphi += 2*TMath::Pi();
 407
 408       Double_t res = sqrt(dphi*dphi + deta*deta);
 409
 410       if (res< cut) {
 411         //add this track-index
 412         size++;
 413         array->Set( size );
 414         array->AddAt( itrk, (size-1) );
 415         if(fDebug>1) printf("MTH:: track %d matched cell at eta=%f , phi=%f \n", itrk, eta, phi);
 416
 417       }
 418     }
 419   }
 420
 421 }
 422
 423 /*
 424 //____________________________________________________________________________
 425 void AliJetAODFillUnitArrayEMCalDigits::GetCalibrationParameters()
 426 {
 427   // Set calibration parameters:
 428   // if calibration database exists, they are read from database,
 429   // otherwise, they are taken from digitizer.
 430   //
 431   // It is a user responsilibity to open CDB before reconstruction,
 432   // for example:
 433   // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
 434
 435   //Check if calibration is stored in data base
 436
 437   if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet()))
 438     {
 439       AliCDBEntry *entry = (AliCDBEntry*)
 440         AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
 441       if (entry) fCalibData =  (AliEMCALCalibData*) entry->GetObject();
 442     }
 443
 444   if(!fCalibData)
 445     printf("************* Calibration parameters not found in CDB! ****************");
 446 //    AliFatal("Calibration parameters not found in CDB!");
 447
 448
 449 }
 450
 451 //____________________________________________________________________________
 452 Float_t  AliJetAODFillUnitArrayEMCalDigits::Calibrate(Int_t amp, Int_t AbsId)
 453 {
 454
 455   // Convert digitized amplitude into energy.
 456   // Calibration parameters are taken from calibration data base for raw data,
 457   // or from digitizer parameters for simulated data.
 458
 459   if(fCalibData){
 460
 461     if (fGeom==0)
 462       printf("************* Did not get geometry from EMCALLoader ***************");
 463 //      AliFatal("Did not get geometry from EMCALLoader") ;
 464
 465     Int_t iSupMod = -1;
 466     Int_t nModule  = -1;
 467     Int_t nIphi   = -1;
 468     Int_t nIeta   = -1;
 469     Int_t iphi    = -1;
 470     Int_t ieta    = -1;
 471
 472     Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
 473     if(!bCell) {
 474       // fGeom->PrintGeometry();
 475       Error("Calibrate()"," Wrong cell id number : %i", AbsId);
 476       assert(0);
 477     }
 478
 479     fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
 480
 481     fADCchannelECA  = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
 482     fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
 483
 484    return -fADCpedestalECA + amp * fADCchannelECA ;
 485
 486   }
 487   else //Return energy with default parameters if calibration is not available
 488     return -fADCpedestalECA + amp * fADCchannelECA ;
 489
 490 }
 491 */
 492
 493