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