JETAN/AliJetESDFillUnitArrayTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //======================================================================
  18 // ***July 2006
  19 // Fill Unit Array class
  20 // Class used by AliJetESDReader to fill a UnitArray from the information extracted
  21 // from the particle tracks
  22 // Author: magali.estienne@ires.in2p3.fr
  23 //======================================================================
  24
  25
  26 // --- ROOT system ---
  27 #include <TVector3.h>
  28 #include <TProcessID.h>
  29
  30 // --- AliRoot header files ---
  31 #include "AliJetUnitArray.h"
  32 #include "AliJetHadronCorrectionv1.h"
  33 #include "AliJetESDFillUnitArrayTracks.h"
  34
  35 // --- ROOT system ---
  36 class TSystem;
  37 class TLorentzVector;
  38 class TGeoManager;
  39
  40 // --- AliRoot header files ---
  41 class AliJetFinder;
  42 class AliJetReader;
  43 class AliJetESDReader;
  44
  45 ClassImp(AliJetESDFillUnitArrayTracks)
  46
  47 //_____________________________________________________________________________
  48 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks()
  49   : AliJetFillUnitArray(),
  50     fNumUnits(0),
  51     fHadCorr(0),
  52     fApplyMIPCorrection(kTRUE),
  53     fESD(0x0),
  54     fGrid0(0x0),
  55     fGrid1(0x0),
  56     fGrid2(0x0),
  57     fGrid3(0x0),
  58     fGrid4(0x0)
  59 {
  60   // constructor
  61 }
  62
  63 //_____________________________________________________________________________
  64 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(AliESDEvent* esd)
  65   : AliJetFillUnitArray(),
  66     fNumUnits(0),
  67     fHadCorr(0),
  68     fApplyMIPCorrection(kTRUE),
  69     fESD(esd),
  70     fGrid0(0x0),
  71     fGrid1(0x0),
  72     fGrid2(0x0),
  73     fGrid3(0x0),
  74     fGrid4(0x0)
  75 {
  76   // constructor
  77 }
  78
  79 //_____________________________________________________________________________
  80 AliJetESDFillUnitArrayTracks::AliJetESDFillUnitArrayTracks(const AliJetESDFillUnitArrayTracks &det)
  81   : AliJetFillUnitArray(det),
  82     fNumUnits(det.fNumUnits),
  83     fHadCorr(det.fHadCorr),
  84     fApplyMIPCorrection(det.fApplyMIPCorrection),
  85     fESD(det.fESD),
  86     fGrid0(det.fGrid0),
  87     fGrid1(det.fGrid1),
  88     fGrid2(det.fGrid2),
  89     fGrid3(det.fGrid3),
  90     fGrid4(det.fGrid4)
  91 {
  92   // Copy constructor
  93 }
  94
  95 //------------------------------------------------------------------
  96 AliJetESDFillUnitArrayTracks& AliJetESDFillUnitArrayTracks::operator=(const AliJetESDFillUnitArrayTracks& other)
  97 {
  98   // Assignment
  99
 100     fNumUnits = other.fNumUnits;
 101     fHadCorr = other.fHadCorr;
 102     fApplyMIPCorrection = other.fApplyMIPCorrection;
 103     fESD = other.fESD;
 104     fGrid0 = other.fGrid0;
 105     fGrid1 = other.fGrid1;
 106     fGrid2 = other.fGrid2;
 107     fGrid3 = other.fGrid3;
 108     fGrid4 = other.fGrid4;
 109
 110     return (*this);
 111 }
 112
 113 //____________________________________________________________________________
 114 void AliJetESDFillUnitArrayTracks::InitParameters()
 115 {
 116   //  fHadCorr        = 0;                 // For hadron correction
 117   fNumUnits = fGeom->GetEMCGeometry()->GetNCells();      // Number of towers in EMCAL
 118   cout << "In AliJetESDFillUnitArrayTracks:InitParameters(), Ncells : " << fNumUnits << endl;
 119
 120   fTPCGrid->GetAccParam(fNphi,fNeta,fPhiMin,
 121                         fPhiMax,fEtaMin,fEtaMax);
 122   fTPCGrid->GetBinParam(fPhiBinInTPCAcc,fEtaBinInTPCAcc,
 123                         fPhiBinInEMCalAcc,fEtaBinInEMCalAcc,fNbinPhi);
 124
 125   fIndex = fTPCGrid->GetIndexObject();
 126
 127   if(fDebug>20){
 128     for(Int_t i=0; i<fNphi+1; i++)
 129       for(Int_t j=0; j<fNeta+1; j++) {cout << "fIndex[" << i << "," << j << "] : " <<
 130           (*fIndex)(i,j) << endl; }
 131   }
 132   if(fDebug>1) printf("\n Parameters initiated ! \n");
 133 }
 134
 135 //_____________________________________________________________________________
 136 AliJetESDFillUnitArrayTracks::~AliJetESDFillUnitArrayTracks()
 137 {
 138   // destructor
 139 }
 140
 141 //_____________________________________________________________________________
 142 void AliJetESDFillUnitArrayTracks::Exec(Option_t* const /*option*/)
 143 {
 144   //
 145   // Main method.
 146   // Fill the unit array with the charged particle information in ESD
 147   //
 148
 149   fDebug = fReaderHeader->GetDebug();
 150
 151   // Set parameters
 152   InitParameters();
 153   //  fRef->Clear();
 154
 155   // get number of tracks in event (for the loop)
 156   Int_t goodTrack = 0;
 157   Int_t nt = 0;
 158 //(not used ?)  Int_t nt2 = 0;
 159   Int_t nmax = 0;
 160   Float_t pt,eta,phi;
 161   //  Int_t sflag = 0;
 162   TVector3 p3;
 163
 164   if(fDebug>1) cout << "fESD in Fill array : " << fESD << endl;
 165
 166   nt = fESD->GetNumberOfTracks();
 167   if(fDebug>1) cout << "Number of Tracks in ESD : " << nt << endl;
 168
 169   // temporary storage of signal and pt cut flag
 170   Int_t* sflag  = new Int_t[nt];
 171   Int_t* cflag  = new Int_t[nt];
 172
 173   // get cuts set by user
 174   Float_t ptMin  = fReaderHeader->GetPtCut();
 175   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 176   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 177   fOpt = fReaderHeader->GetDetector();
 178   fDZ  = fReaderHeader->GetDZ();
 179
 180   Int_t nTracksEmcal      = 0;
 181   Int_t nTracksEmcalDZ    = 0;
 182   Int_t nTracksTpc        = 0;
 183   Int_t nTracksTpcOnly    = 0;
 184   Int_t nTracksEmcalCut   = 0;
 185   Int_t nTracksEmcalDZCut = 0;
 186   Int_t nTracksTpcCut     = 0;
 187   Int_t nTracksTpcOnlyCut = 0;
 188
 189 //(not used ?)  Int_t nElements  = fTPCGrid->GetNEntries();
 190 //(not used ?)  Int_t nElements2 = fEMCalGrid->GetNEntries();
 191   fGrid = fTPCGrid->GetGridType();
 192
 193 //(not used ?)  Int_t nEntries = fUnitArray->GetEntries();
 194
 195 //(not used ?)  Int_t nRefEnt = fRefArray->GetEntries();
 196
 197   //loop over tracks
 198     nmax = nt;
 199     for (Int_t it = 0; it < nmax; it++) {
 200       AliESDtrack *track;
 201       track = fESD->GetTrack(it);
 202       UInt_t status = track->GetStatus();
 203
 204       Double_t mom[3];
 205       track->GetPxPyPz(mom);
 206
 207       p3.SetXYZ(mom[0],mom[1],mom[2]);
 208       pt = p3.Pt();
 209
 210       Float_t mass = 0.;
 211       mass = track->GetMass();
 212
 213       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
 214       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
 215       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 216           && TMath::Abs(track->GetLabel()) > 10000)  continue;   // quality check
 217       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 218           && TMath::Abs(track->GetLabel()) < 10000)  continue;   // quality check
 219        eta = p3.Eta();
 220        phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
 221
 222       if ( (eta > etaMax) || (eta < etaMin)) continue;           // checking eta cut
 223
 224       sflag[goodTrack]=0;
 225       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 226       cflag[goodTrack]=0;
 227       if (pt > ptMin) cflag[goodTrack]=1;                       // pt cut
 228       //      fRef->Add(track);
 229
 230       if(fGrid==0)
 231         {
 232           // Only TPC filled from its grid in its total acceptance
 233
 234           Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 235           Bool_t ok = kFALSE;
 236           Bool_t ref = kFALSE;
 237
 238           AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(idTPC-1);
 239           TRefArray *reference = uArray->GetUnitTrackRef();
 240           if (reference->GetEntries() == 0)  {
 241               new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
 242           }
 243
 244           reference->Add(track);
 245
 246           Float_t unitEnergy = 0.;
 247           unitEnergy = uArray->GetUnitEnergy();
 248           // nTracksTpcOnly is necessary to count the number of candidate cells
 249           // but it doesn't give the real multiplicity -> it will be extracted
 250           // in the jet finder through the reference to tracks
 251           if(unitEnergy==0.){
 252             nTracksTpcOnly++;
 253             ok = kTRUE;
 254             ref = kTRUE;
 255           }
 256
 257           // Fill energy in TPC acceptance
 258           uArray->SetUnitEnergy(unitEnergy + pt);
 259           uArray->SetUnitMass(mass);
 260
 261           // Pt cut flag
 262           if(uArray->GetUnitEnergy()<ptMin){
 263             uArray->SetUnitCutFlag(kPtSmaller);
 264           }
 265           else {
 266             uArray->SetUnitCutFlag(kPtHigher);
 267             if(ok) nTracksTpcOnlyCut++;
 268           }
 269
 270           // Signal flag
 271           if(sflag[goodTrack] == 1) {
 272             uArray->SetUnitSignalFlag(kGood);
 273             uArray->SetUnitSignalFlagC(kFALSE,kGood);
 274           } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 275
 276           if(uArray->GetUnitEnergy()>0 && ref){
 277             if(!fProcId) {
 278               fProcId = kTRUE;
 279               new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 280             }
 281             fRefArray->Add(uArray);
 282           }
 283         }
 284
 285       if(fGrid==1)
 286         {
 287           Int_t nElements = fTPCGrid->GetNEntries();
 288           // Fill track information in EMCAL acceptance
 289             if((eta >= fEtaMin && eta <= fEtaMax) &&
 290              (phi >= fPhiMin && phi <= fPhiMax))// &&
 291             {
 292
 293               // Include dead-zones
 294               if(fDZ)
 295                 {
 296                   Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 297                   Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 298                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 299                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 300                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 301                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 302                   fGeom->GetEMCGeometry()->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 303                   Int_t n0 = fGrid0->GetNEntries();
 304                   Int_t n1 = fGrid1->GetNEntries();
 305                   Int_t n2 = fGrid2->GetNEntries();
 306                   Int_t n3 = fGrid3->GetNEntries();
 307 //(not used ?)            Int_t n4 = fGrid4->GetNEntries();
 308
 309                   if(phi >= phimin0 && phi <= phimax0){
 310                     Int_t id0 = fGrid0->GetIndex(phi,eta)-1;
 311                     AliJetUnitArray *uArray0 = (AliJetUnitArray*)fUnitArray->At(id0+fNumUnits+nElements);
 312                     TRefArray *reference0 = uArray0->GetUnitTrackRef();
 313                     if (reference0->GetEntries() == 0)  {
 314                         new(reference0) TRefArray(TProcessID::GetProcessWithUID(track));
 315                     }
 316
 317                     reference0->Add(track);
 318
 319                     Float_t uEnergy0 = uArray0->GetUnitEnergy();
 320                     Bool_t ok0 = kFALSE;
 321                     Bool_t ref0 = kFALSE;
 322                     if(uEnergy0==0.){
 323                       nTracksEmcalDZ++;
 324                       ok0 = kTRUE;
 325                       ref0 = kTRUE;
 326                     }
 327                     uArray0->SetUnitEnergy(uEnergy0+pt);
 328
 329                     if(uArray0->GetUnitEnergy()<ptMin)
 330                       uArray0->SetUnitCutFlag(kPtSmaller);
 331                     else {
 332                       uArray0->SetUnitCutFlag(kPtHigher);
 333                       if(ok0) nTracksEmcalDZCut++;
 334                     }
 335                     if(sflag[goodTrack] == 1) {
 336                       uArray0->SetUnitSignalFlag(kGood);
 337                       uArray0->SetUnitSignalFlagC(kFALSE,kGood);
 338                     } else uArray0->SetUnitSignalFlagC(kFALSE,kBad);
 339
 340                     if(uArray0->GetUnitEnergy()>0 && ref0){
 341                       if(!fProcId){
 342                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray0));
 343                         fProcId = kTRUE;
 344                       }
 345                       fRefArray->Add(uArray0);
 346                     }
 347                   }
 348
 349                   if(phi >= phimin1 && phi <= phimax1){
 350                     Int_t id1 = fGrid1->GetIndex(phi,eta)-1+n0;
 351                     AliJetUnitArray *uArray1 = (AliJetUnitArray*)fUnitArray->At(id1+fNumUnits+nElements);
 352                     TRefArray *reference1 = uArray1->GetUnitTrackRef();
 353                     if (reference1->GetEntries() == 0)  {
 354                         new(reference1) TRefArray(TProcessID::GetProcessWithUID(track));
 355                     }
 356
 357                     reference1->Add(track);
 358
 359                     Float_t uEnergy1 = uArray1->GetUnitEnergy();
 360                     Bool_t ok1 = kFALSE;
 361                     Bool_t ref1 = kFALSE;
 362                     if(uEnergy1==0.){
 363                       nTracksEmcalDZ++;
 364                       ok1 = kTRUE;
 365                       ref1 = kTRUE;
 366                     }
 367                     uArray1->SetUnitEnergy(uEnergy1+pt);
 368
 369                     if(uArray1->GetUnitEnergy()<ptMin)
 370                       uArray1->SetUnitCutFlag(kPtSmaller);
 371                     else {
 372                       uArray1->SetUnitCutFlag(kPtHigher);
 373                       if(ok1) nTracksEmcalDZCut++;
 374                     }
 375                     if(sflag[goodTrack] == 1) {
 376                       uArray1->SetUnitSignalFlag(kGood);
 377                       uArray1->SetUnitSignalFlagC(kFALSE,kGood);
 378                     } else uArray1->SetUnitSignalFlagC(kFALSE,kBad);
 379
 380                     if(uArray1->GetUnitEnergy()>0 && ref1){
 381                       if(!fProcId){
 382                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray1));
 383                         fProcId = kTRUE;
 384                       }
 385                       fRefArray->Add(uArray1);
 386                     }
 387                   }
 388
 389                   if(phi >= phimin2 && phi <= phimax2){
 390                     Int_t id2 = fGrid2->GetIndex(phi,eta)-1+n0+n1;
 391                     AliJetUnitArray *uArray2 = (AliJetUnitArray*)fUnitArray->At(id2+fNumUnits+nElements);
 392                     TRefArray *reference2 = uArray2->GetUnitTrackRef();
 393                     if (reference2->GetEntries() == 0)  {
 394                         new(reference2) TRefArray(TProcessID::GetProcessWithUID(track));
 395                     }
 396
 397                     reference2->Add(track);
 398
 399                     Float_t uEnergy2 = uArray2->GetUnitEnergy();
 400                     Bool_t ok2 = kFALSE;
 401                     Bool_t ref2 = kFALSE;
 402                     if(uEnergy2==0.){
 403                       nTracksEmcalDZ++;
 404                       ok2 = kTRUE;
 405                       ref2 = kTRUE;
 406                     }
 407                     uArray2->SetUnitEnergy(uEnergy2+pt);
 408
 409                     if(uArray2->GetUnitEnergy()<ptMin)
 410                       uArray2->SetUnitCutFlag(kPtSmaller);
 411                     else {
 412                       uArray2->SetUnitCutFlag(kPtHigher);
 413                       if(ok2) nTracksEmcalDZCut++;
 414                     }
 415                     if(sflag[goodTrack] == 1) {
 416                       uArray2->SetUnitSignalFlag(kGood);
 417                       uArray2->SetUnitSignalFlagC(kFALSE,kGood);
 418                     } else uArray2->SetUnitSignalFlagC(kFALSE,kBad);
 419
 420                     if(uArray2->GetUnitEnergy()>0 && ref2){
 421                       if(!fProcId){
 422                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray2));
 423                         fProcId = kTRUE;
 424                       }
 425                       fRefArray->Add(uArray2);
 426                     }
 427                   }
 428
 429                   if(phi >= phimin3 && phi <= phimax3){
 430                     Int_t id3 = fGrid3->GetIndex(phi,eta)-1+n0+n1+n2;
 431                     AliJetUnitArray *uArray3 = (AliJetUnitArray*)fUnitArray->At(id3+fNumUnits+nElements);
 432                     TRefArray *reference3 = uArray3->GetUnitTrackRef();
 433                     if (reference3->GetEntries() == 0)  {
 434                         new(reference3) TRefArray(TProcessID::GetProcessWithUID(track));
 435                     }
 436
 437                     reference3->Add(track);
 438
 439                     Float_t uEnergy3 = uArray3->GetUnitEnergy();
 440                     Bool_t ok3 = kFALSE;
 441                     Bool_t ref3 = kFALSE;
 442                     if(uEnergy3==0.){
 443                       nTracksEmcalDZ++;
 444                       ok3 = kTRUE;
 445                       ref3 = kTRUE;
 446                     }
 447                     uArray3->SetUnitEnergy(uEnergy3+pt);
 448
 449                     if(uArray3->GetUnitEnergy()<ptMin)
 450                       uArray3->SetUnitCutFlag(kPtSmaller);
 451                     else {
 452                       uArray3->SetUnitCutFlag(kPtHigher);
 453                       if(ok3) nTracksEmcalDZCut++;
 454                     }
 455                     if(sflag[goodTrack] == 1) {
 456                       uArray3->SetUnitSignalFlag(kGood);
 457                       uArray3->SetUnitSignalFlagC(kFALSE,kGood);
 458                     } else uArray3->SetUnitSignalFlagC(kFALSE,kBad);
 459
 460                     if(uArray3->GetUnitEnergy()>0 && ref3){
 461                       if(!fProcId){
 462                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray3));
 463                         fProcId = kTRUE;
 464                       }
 465                       fRefArray->Add(uArray3);
 466                     }
 467                   }
 468
 469                   if(phi >= phimin4 && phi <= phimax4){
 470                     Int_t id4 = fGrid4->GetIndex(phi,eta)-1+n0+n1+n2+n3;
 471                     AliJetUnitArray *uArray4 = (AliJetUnitArray*)fUnitArray->At(id4+fNumUnits+nElements);
 472                     TRefArray *reference4 = uArray4->GetUnitTrackRef();
 473                     if (reference4->GetEntries() == 0)  {
 474                         new(reference4) TRefArray(TProcessID::GetProcessWithUID(track));
 475                     }
 476
 477                     reference4->Add(track);
 478
 479                     Float_t uEnergy4 = uArray4->GetUnitEnergy();
 480                     Bool_t ok4 = kFALSE;
 481                     Bool_t ref4 = kFALSE;
 482                     if(uEnergy4==0.){
 483                       nTracksEmcalDZ++;
 484                       ok4 = kTRUE;
 485                       ref4 = kTRUE;
 486                     }
 487                     uArray4->SetUnitEnergy(uEnergy4+pt);
 488
 489                     if(uArray4->GetUnitEnergy()<ptMin)
 490                       uArray4->SetUnitCutFlag(kPtSmaller);
 491                     else {
 492                       uArray4->SetUnitCutFlag(kPtHigher);
 493                       if(ok4) nTracksEmcalDZCut++;
 494                     }
 495                     if(sflag[goodTrack] == 1) {
 496                       uArray4->SetUnitSignalFlag(kGood);
 497                       uArray4->SetUnitSignalFlagC(kFALSE,kGood);
 498                     } else uArray4->SetUnitSignalFlagC(kFALSE,kBad);
 499
 500                     if(uArray4->GetUnitEnergy()>0 && ref4){
 501                       if(!fProcId){
 502                         new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray4));
 503                         fProcId = kTRUE;
 504                       }
 505                       fRefArray->Add(uArray4);
 506                     }
 507                   }
 508                 } // end fDZ
 509
 510               Int_t towerID = 0;
 511               fGeom->GetAbsCellIdFromEtaPhi(eta,phi,towerID);
 512               if(towerID==-1) continue;
 513
 514               AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(towerID);
 515               TRefArray *reference = uArray->GetUnitTrackRef();
 516               if (reference->GetEntries() == 0)  {
 517                   new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
 518               }
 519
 520               reference->Add(track);
 521
 522               Float_t unitEnergy = uArray->GetUnitEnergy();
 523               Bool_t ok = kFALSE;
 524               Bool_t ref = kFALSE;
 525               if(unitEnergy==0.){
 526                 nTracksEmcal++;
 527                 ok=kTRUE;
 528                 ref=kTRUE;
 529               }
 530
 531               // Do Hadron Correction
 532               // This is under construction !!!!!!!!!!!!!!!!!!!!!!!
 533               // For the moment I apply MIP correction if p >= 0.5 GeV/c
 534               if (fApplyMIPCorrection != 0 && p3.Mag() >= 0.5)
 535                 {
 536                   ((AliJetHadronCorrectionv1*)fHadCorr)->SetGeometry("EMCAL_COMPLETE",1.);
 537
 538                   // Get track position at the outer part of the reconstruction ~ TRD
 539                   Double_t phiOut = track->GetOuterParam()->Phi();
 540                   Double_t etaOut = track->GetOuterParam()->Eta();
 541
 542                   // If the track in the outer part of the TPC/TDR ? is inside
 543                   // the calorimeter, it can deposit part of its energy
 544                   // We can then correct on average for these particles
 545                   if((etaOut >= fEtaMin && etaOut <= fEtaMax) &&
 546                      (phiOut >= fPhiMin && phiOut <= fPhiMax))// &&
 547                     {
 548                       Double_t   hCEnergy = (Double_t)fHadCorr->GetEnergy(p3.Mag(), (Double_t)eta,0);
 549                       unitEnergy -= hCEnergy*TMath::Sin(2.0*TMath::ATan(TMath::Exp(-eta)));
 550                     }
 551                 } //end Hadron Correction loop
 552
 553               cout << "unitEnergy + pt = " << unitEnergy << " + " << pt << " = " << unitEnergy + pt << endl;
 554
 555               if((unitEnergy + pt) > 0.) uArray->SetUnitEnergy(unitEnergy + pt);
 556               else uArray->SetUnitEnergy(0.);
 557
 558               // Put a pt cut flag
 559               if(uArray->GetUnitEnergy()<ptMin){
 560                 uArray->SetUnitCutFlag(kPtSmaller);
 561               }
 562               else {
 563                 uArray->SetUnitCutFlag(kPtHigher);
 564                 if(ok) nTracksEmcalCut++;
 565               }
 566
 567               // Signal flag
 568               if(sflag[goodTrack] == 1) {
 569                 uArray->SetUnitSignalFlag(kGood);
 570                 uArray->SetUnitSignalFlagC(kFALSE,kGood);
 571               } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 572
 573
 574               if(uArray->GetUnitEnergy()>0 && ref){
 575                 if(!fProcId){
 576                   new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 577                   fProcId = kTRUE;
 578                 }
 579                 fRefArray->Add(uArray);
 580               }
 581             } // end loop on EMCal acceptance cut + tracks quality
 582           else{
 583             // Outside EMCal acceptance
 584
 585             Int_t idTPC = fTPCGrid->GetIndex(phi,eta);
 586
 587             AliJetUnitArray *uArray = (AliJetUnitArray*)fUnitArray->At(fNumUnits-1+idTPC);
 588             TRefArray *reference = uArray->GetUnitTrackRef();
 589             if (reference->GetEntries() == 0)  {
 590                 new(reference) TRefArray(TProcessID::GetProcessWithUID(track));
 591             }
 592
 593             reference->Add(track);
 594
 595             Float_t unitEnergy2 = uArray->GetUnitEnergy();
 596             Bool_t okout = kFALSE;
 597             Bool_t refout = kFALSE;
 598             if(unitEnergy2==0.){
 599               nTracksTpc++;
 600               okout=kTRUE;
 601               refout=kTRUE;
 602             }
 603             // Fill energy outside emcal acceptance
 604             uArray->SetUnitEnergy(unitEnergy2 + pt);
 605
 606             // Pt cut flag
 607             if(uArray->GetUnitEnergy()<ptMin){
 608               uArray->SetUnitCutFlag(kPtSmaller);
 609             }
 610             else {
 611               uArray->SetUnitCutFlag(kPtHigher);
 612               if(okout) nTracksTpcCut++;
 613             }
 614             // Signal flag
 615             if(sflag[goodTrack] == 1) {
 616               uArray->SetUnitSignalFlag(kGood);
 617               uArray->SetUnitSignalFlagC(kFALSE,kGood);
 618             } else uArray->SetUnitSignalFlagC(kFALSE,kBad);
 619
 620
 621             if(uArray->GetUnitEnergy()>0 && refout){
 622               if(!fProcId){
 623                 new(fRefArray) TRefArray(TProcessID::GetProcessWithUID(uArray));
 624                 fProcId = kTRUE;
 625               }
 626               fRefArray->Add(uArray);
 627             }
 628           }
 629         } // end fGrid==1
 630
 631       goodTrack++;
 632
 633     } // end loop on entries (tpc tracks)
 634
 635   if(fDebug>0)
 636     {
 637       cout << "End of Tracks %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
 638       cout << "goodTracks: " << goodTrack << endl;
 639     }
 640
 641   //  fSignalFlag.Set(goodTrack,sflag);
 642   //  fCutFlag.Set(goodTrack,cflag);
 643
 644   delete[] sflag;
 645   delete[] cflag;
 646
 647   if(fGrid==0) {
 648     fNTracks = nTracksTpcOnly;
 649     fNTracksCut = nTracksTpcOnlyCut;
 650     if(fDebug>10){
 651       cout << "fNTracks : " << fNTracks << endl;
 652       cout << "fNTracksCut : " << fNTracksCut << endl;
 653     }
 654   }
 655   if(fGrid==1) {
 656     fNTracks = nTracksEmcal+nTracksEmcalDZ+nTracksTpc;
 657     fNTracksCut = nTracksEmcalCut+nTracksEmcalDZCut+nTracksTpcCut;
 658     if(fDebug>10){
 659       cout << "fNTracks : " << fNTracks << endl;
 660       cout << "fNTracksCut : " << fNTracksCut << endl;
 661     }
 662   }
 663
 664 }
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678