JETAN/AliJetESDReader.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 // Jet ESD Reader
  18 // ESD reader for jet analysis
  19 // Authors: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
  20 //          Magali Estienne <magali.estienne@subatech.in2p3.fr>
  21 //-------------------------------------------------------------------------
  22
  23 // --- Standard library ---
  24 #include <Riostream.h>
  25
  26 // --- ROOT system ---
  27 #include <TSystem.h>
  28 #include <TStopwatch.h>
  29 #include <TLorentzVector.h>
  30 #include <TVector3.h>
  31 #include "TTask.h"
  32 #include "TTree.h"
  33 #include "TFile.h"
  34 #include <TGeoManager.h>
  35 #include <assert.h>
  36 #include <TRefArray.h>
  37 #include <TMath.h>
  38 #include <TProcessID.h>
  39 #include <TRandom3.h>
  40
  41 // --- AliRoot header files ---
  42 #include "AliJetESDReader.h"
  43 #include "AliJetESDReaderHeader.h"
  44 #include "AliESDEvent.h"
  45 #include "AliESD.h"
  46 #include "AliESDtrack.h"
  47 #include "AliJetDummyGeo.h"
  48 #include "AliJetESDFillUnitArrayTracks.h"
  49 #include "AliJetESDFillUnitArrayEMCalDigits.h"
  50 #include "AliJetHadronCorrectionv1.h"
  51 #include "AliJetUnitArray.h"
  52 #include "AliAnalysisTask.h"
  53
  54 ClassImp(AliJetESDReader)
  55
  56 AliJetESDReader::AliJetESDReader():
  57   AliJetReader(),
  58   fGeom(0),
  59   fHadCorr(0x0),
  60   fTpcGrid(0x0),
  61   fEmcalGrid(0x0),
  62   fGrid0(0),
  63   fGrid1(0),
  64   fGrid2(0),
  65   fGrid3(0),
  66   fGrid4(0),
  67   fApplyElectronCorrection(kFALSE),
  68   fApplyMIPCorrection(kTRUE),
  69   fApplyFractionHadronicCorrection(kFALSE),
  70   fFractionHadronicCorrection(0.3),
  71   fNumUnits(0),
  72   fDebug(0),
  73   fMass(0),
  74   fSign(0),
  75   fNIn(0),
  76   fOpt(0),
  77   fDZ(0),
  78   fNeta(0),
  79   fNphi(0),
  80   fRefArray(0x0),
  81   fProcId(kFALSE)
  82 {
  83   // Constructor
  84 }
  85
  86 //____________________________________________________________________________
  87 AliJetESDReader::~AliJetESDReader()
  88 {
  89   // Destructor
  90     delete fChain;
  91     delete fTree;
  92     delete fESD;
  93     delete fTpcGrid;
  94     delete fEmcalGrid;
  95     if(fDZ)
  96       {
  97         delete fGrid0;
  98         delete fGrid1;
  99         delete fGrid2;
 100         delete fGrid3;
 101         delete fGrid4;
 102       }
 103 }
 104
 105 //____________________________________________________________________________
 106 void AliJetESDReader::OpenInputFiles()
 107 {
 108   // chain for the ESDs
 109   fChain   = new TChain("esdTree");
 110
 111   // get directory and pattern name from the header
 112    const char* dirName=fReaderHeader->GetDirectory();
 113    const char* pattern=fReaderHeader->GetPattern();
 114
 115    // Add files matching patters to the chain
 116    void *dir  = gSystem->OpenDirectory(dirName);
 117    const char *name = 0x0;
 118    int nesd = ((AliJetESDReaderHeader*) fReaderHeader)->GetNesd();
 119    int a = 0;
 120    while ((name = gSystem->GetDirEntry(dir))){
 121        if (a>=nesd) continue;
 122
 123        if (strstr(name,pattern)){
 124          printf("Adding %s\n",name);
 125          fChain->AddFile(Form("%s/%s/AliESDs.root",dirName,name));
 126          a++;
 127        }
 128    }
 129
 130   gSystem->FreeDirectory(dir);
 131
 132   fESD = new AliESDEvent();
 133   fESD->ReadFromTree(fChain);
 134
 135   int nMax = fChain->GetEntries();
 136
 137   printf("\n AliJetESDReader: Total number of events in chain= %d \n",nMax);
 138
 139   // set number of events in header
 140   if (fReaderHeader->GetLastEvent() == -1)
 141     fReaderHeader->SetLastEvent(nMax);
 142   else {
 143     Int_t nUsr = fReaderHeader->GetLastEvent();
 144     fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
 145   }
 146
 147 }
 148
 149 //__________________________________________________________
 150 void AliJetESDReader::SetApplyMIPCorrection(Bool_t val)
 151 {
 152   //
 153   // Set flag to apply MIP correction fApplyMIPCorrection
 154   // - exclusive with fApplyFractionHadronicCorrection
 155   //
 156
 157   fApplyMIPCorrection = val;
 158   if(fApplyMIPCorrection == kTRUE)
 159     {
 160       SetApplyFractionHadronicCorrection(kFALSE);
 161       printf("Enabling MIP Correction \n");
 162     }
 163   else
 164     {
 165       printf("Disabling MIP Correction \n");
 166     }
 167 }
 168
 169 //__________________________________________________________
 170 void AliJetESDReader::SetApplyFractionHadronicCorrection(Bool_t val)
 171 {
 172   //
 173   // Set flag to apply EMC hadronic correction fApplyFractionHadronicCorrection
 174   // - exclusive with fApplyMIPCorrection
 175   //
 176
 177   fApplyFractionHadronicCorrection = val;
 178   if(fApplyFractionHadronicCorrection == kTRUE)
 179     {
 180       SetApplyMIPCorrection(kFALSE);
 181       printf("Enabling Fraction Hadronic Correction \n");
 182     }
 183   else
 184     {
 185       printf("Disabling Fraction Hadronic Correction \n");
 186     }
 187 }
 188
 189 //__________________________________________________________
 190 void AliJetESDReader::SetFractionHadronicCorrection(Double_t val)
 191 {
 192   //
 193   // Set value to fFractionHadronicCorrection (default is 0.3)
 194   // apply EMC hadronic correction fApplyFractionHadronicCorrection
 195   // - exclusive with fApplyMIPCorrection
 196   //
 197
 198   fFractionHadronicCorrection = val;
 199   if(fFractionHadronicCorrection > 0.0 && fFractionHadronicCorrection <= 1.0)
 200     {
 201       SetApplyFractionHadronicCorrection(kTRUE);
 202       printf("Fraction Hadronic Correction %1.3f \n",fFractionHadronicCorrection);
 203     }
 204   else
 205     {
 206       SetApplyFractionHadronicCorrection(kFALSE);
 207     }
 208 }
 209
 210
 211 //____________________________________________________________________________
 212 void AliJetESDReader::SetInputEvent(const TObject* esd, const TObject* /*aod*/, const TObject* /*mc*/) {
 213     // Connect the tree
 214      fESD   = (AliESDEvent*) esd;
 215 }
 216
 217 //____________________________________________________________________________
 218 Bool_t AliJetESDReader::FillMomentumArray()
 219 {
 220   // Fill momentum array
 221
 222   Int_t goodTrack = 0;
 223   Int_t nt = 0;
 224   Float_t pt, eta;
 225   TVector3 p3;
 226
 227   // clear momentum array
 228   ClearArray();
 229   fDebug = fReaderHeader->GetDebug();
 230   fOpt = fReaderHeader->GetDetector();
 231
 232   if (!fESD) {
 233       return kFALSE;
 234   }
 235
 236   // get number of tracks in event (for the loop)
 237   nt = fESD->GetNumberOfTracks();
 238   printf("Fill Momentum Array %5d  \n", nt);
 239
 240   // temporary storage of signal and pt cut flag
 241   Int_t* sflag  = new Int_t[nt];
 242   Int_t* cflag  = new Int_t[nt];
 243
 244   // get cuts set by user
 245   Float_t ptMin =  fReaderHeader->GetPtCut();
 246   Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
 247   Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
 248
 249   //loop over tracks in ESD
 250   for (Int_t it = 0; it < nt; it++) {
 251       AliESDtrack *track = fESD->GetTrack(it);
 252       UInt_t status = track->GetStatus();
 253
 254       Double_t mom[3];
 255       track->GetPxPyPz(mom);
 256       p3.SetXYZ(mom[0],mom[1],mom[2]);
 257       pt = p3.Pt();
 258       if ((status & AliESDtrack::kTPCrefit) == 0)    continue;      // quality check
 259       if ((status & AliESDtrack::kITSrefit) == 0)    continue;      // quality check
 260       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
 261           && TMath::Abs(track->GetLabel()) > 10000)  continue;      // quality check
 262       if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
 263           && TMath::Abs(track->GetLabel()) < 10000)  continue;      // quality check
 264       eta = p3.Eta();
 265       if ( (eta > etaMax) || (eta < etaMin))         continue;      // checking eta cut
 266
 267       new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
 268       sflag[goodTrack]=0;
 269       if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
 270       cflag[goodTrack]=0;
 271       if (pt > ptMin) cflag[goodTrack]=1;                           // pt cut
 272       goodTrack++;
 273   }
 274
 275   // set the signal flags
 276   fSignalFlag.Set(goodTrack,sflag);
 277   fCutFlag.Set(goodTrack,cflag);
 278
 279   delete[] sflag;
 280   delete[] cflag;
 281
 282   return kTRUE;
 283
 284 }
 285
 286 //____________________________________________________________________________
 287 void AliJetESDReader::CreateTasks(TChain* tree)
 288 {
 289   //
 290   // For reader task initialization
 291   //
 292
 293   fDebug = fReaderHeader->GetDebug();
 294   fDZ = fReaderHeader->GetDZ();
 295   fTree = tree;
 296
 297   // Init EMCAL geometry
 298   SetEMCALGeometry();
 299   // Init parameters
 300   InitParameters();
 301   // Create and init unit array
 302   InitUnitArray();
 303
 304   // Create global reader task for analysis
 305   fFillUnitArray = new TTask("fFillUnitArray","Fill unit array jet finder");
 306   // Create a task for to fill the charged particle information
 307   fFillUAFromTracks = new AliJetESDFillUnitArrayTracks();
 308   fFillUAFromTracks->SetReaderHeader(fReaderHeader);
 309   fFillUAFromTracks->SetGeom(fGeom);
 310   fFillUAFromTracks->SetTPCGrid(fTpcGrid);
 311   fFillUAFromTracks->SetEMCalGrid(fEmcalGrid);
 312   if(fDZ)
 313     { // Calo dead zones inclusion
 314       fFillUAFromTracks->SetGrid0(fGrid0);
 315       fFillUAFromTracks->SetGrid1(fGrid1);
 316       fFillUAFromTracks->SetGrid2(fGrid2);
 317       fFillUAFromTracks->SetGrid3(fGrid3);
 318       fFillUAFromTracks->SetGrid4(fGrid4);
 319     }
 320   fFillUAFromTracks->SetApplyMIPCorrection(fApplyMIPCorrection);
 321   fFillUAFromTracks->SetHadCorrector(fHadCorr);
 322   // Create a task for to fill the neutral particle information
 323   fFillUAFromEMCalDigits = new AliJetESDFillUnitArrayEMCalDigits();
 324   fFillUAFromEMCalDigits->SetReaderHeader(fReaderHeader);
 325   fFillUAFromEMCalDigits->SetGeom(fGeom);
 326   fFillUAFromEMCalDigits->SetTPCGrid(fTpcGrid);
 327   fFillUAFromEMCalDigits->SetEMCalGrid(fEmcalGrid);
 328   fFillUAFromEMCalDigits->SetApplyFractionHadronicCorrection(fApplyFractionHadronicCorrection);
 329   fFillUAFromEMCalDigits->SetFractionHadronicCorrection(fFractionHadronicCorrection);
 330    fFillUAFromEMCalDigits->SetApplyElectronCorrection(fApplyElectronCorrection);
 331   // Add the task to global task
 332   fFillUnitArray->Add(fFillUAFromTracks);
 333   fFillUnitArray->Add(fFillUAFromEMCalDigits);
 334   fFillUAFromTracks->SetActive(kFALSE);
 335   fFillUAFromEMCalDigits->SetActive(kFALSE);
 336
 337   cout << "Tasks instantiated at that stage ! " << endl;
 338   cout << "You can loop over events now ! " << endl;
 339
 340 }
 341
 342 //____________________________________________________________________________
 343 Bool_t AliJetESDReader::ExecTasks(Bool_t procid, TRefArray* refArray)
 344 {
 345   //
 346   // Reader task execussion
 347   //
 348
 349   fProcId = procid;
 350   fRefArray = refArray;
 351
 352   // clear momentum array
 353   ClearArray();
 354
 355   fDebug = fReaderHeader->GetDebug();
 356   fOpt = fReaderHeader->GetDetector();
 357
 358   if(!fESD) {
 359     return kFALSE;
 360   }
 361
 362   // TPC only or Digits+TPC or Clusters+TPC
 363   if(fOpt%2==!0 && fOpt!=0){
 364     fFillUAFromTracks->SetESD(fESD);
 365     fFillUAFromTracks->SetActive(kTRUE);
 366     fFillUAFromTracks->SetUnitArray(fUnitArray);
 367     fFillUAFromTracks->SetRefArray(fRefArray);
 368     fFillUAFromTracks->SetProcId(fProcId);
 369     //    fFillUAFromTracks->ExecuteTask("tpc"); // Temporarily changed
 370     fFillUAFromTracks->Exec("tpc");              // Temporarily added
 371     if(fOpt==1){
 372       fNumCandidate = fFillUAFromTracks->GetMult();
 373       fNumCandidateCut = fFillUAFromTracks->GetMultCut();
 374     }
 375   }
 376
 377   // Digits only or Digits+TPC
 378   if(fOpt>=2 && fOpt<=3){
 379     fFillUAFromEMCalDigits->SetESD(fESD);
 380     fFillUAFromEMCalDigits->SetActive(kTRUE);
 381     fFillUAFromEMCalDigits->SetUnitArray(fUnitArray);
 382     fFillUAFromEMCalDigits->SetRefArray(fRefArray);
 383     fFillUAFromEMCalDigits->SetProcId(fFillUAFromTracks->GetProcId());
 384     fFillUAFromEMCalDigits->SetInitMult(fFillUAFromTracks->GetMult());
 385     fFillUAFromEMCalDigits->SetInitMultCut(fFillUAFromTracks->GetMultCut());
 386     fFillUAFromEMCalDigits->Exec("digits");      // Temporarily added
 387     fNumCandidate = fFillUAFromEMCalDigits->GetMult();
 388     fNumCandidateCut = fFillUAFromEMCalDigits->GetMultCut();
 389   }
 390
 391   //  fFillUnitArray->ExecuteTask();             // Temporarily commented
 392
 393   return kTRUE;
 394 }
 395
 396 //____________________________________________________________________________
 397 Bool_t AliJetESDReader::SetEMCALGeometry()
 398 {
 399   //
 400   // Set the EMCal Geometry
 401   //
 402
 403   if (!fTree->GetFile())
 404     return kFALSE;
 405
 406   TString geomFile(fTree->GetFile()->GetName());
 407   geomFile.ReplaceAll("AliESDs", "geometry");
 408
 409   // temporary workaround for PROOF bug #18505
 410   geomFile.ReplaceAll("#geometry.root#geometry.root", "#geometry.root");
 411   if(fDebug>1) printf("Current geometry file %s \n", geomFile.Data());
 412
 413   // Define EMCAL geometry to be able to read ESDs
 414   fGeom = AliJetDummyGeo::GetInstance();
 415   if (fGeom == 0)
 416     fGeom = AliJetDummyGeo::GetInstance("EMCAL_COMPLETE","EMCAL");
 417
 418   // To be setted to run some AliEMCALGeometry functions
 419   TGeoManager::Import(geomFile);
 420   fGeom->GetTransformationForSM();
 421   printf("\n EMCal Geometry set ! \n");
 422
 423   return kTRUE;
 424 }
 425
 426 //____________________________________________________________________________
 427 void AliJetESDReader::InitParameters()
 428 {
 429   // Initialise parameters
 430   fOpt = fReaderHeader->GetDetector();
 431   //  fHadCorr        = 0;          // For hadron correction
 432   if(fEFlag==kFALSE){
 433     if(fOpt==0 || fOpt==1)
 434       fECorrection    = 0;        // For electron correction
 435     else fECorrection = 1;        // For electron correction
 436   }
 437   fNumUnits       = fGeom->GetNCells();      // Number of cells in EMCAL
 438   if(fDebug>1) printf("\n EMCal parameters initiated ! \n");
 439 }
 440
 441 //____________________________________________________________________________
 442 void AliJetESDReader::InitUnitArray()
 443 {
 444   //
 445   // Create and Initialise unit arrays
 446   //
 447
 448   Int_t nElements = fTpcGrid->GetNEntries();
 449   Float_t eta = 0., phi = 0., Deta = 0., Dphi = 0.;
 450   if(fArrayInitialised) fUnitArray->Delete();
 451
 452   if(fTpcGrid->GetGridType()==0)
 453     { // Fill the following quantities :
 454       // Good track ID, (Eta,Phi) position ID, eta, phi, energy, px, py, pz, Deta, Dphi,
 455       // detector flag, in/out jet, pt cut, mass, cluster ID)
 456       for(Int_t nBin = 1; nBin < nElements+1; nBin++)
 457         {
 458           fTpcGrid->GetEtaPhiFromIndex2(nBin,phi,eta);
 459           phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 460           Deta = fTpcGrid->GetDeta();
 461           Dphi = fTpcGrid->GetDphi();
 462           new ((*fUnitArray)[nBin-1]) AliJetUnitArray(nBin-1,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 463         }
 464     }
 465
 466   if(fTpcGrid->GetGridType()==1)
 467     {
 468       Int_t nGaps = 0;
 469       Int_t n0 = 0, n1 = 0, n2 = 0, n3 = 0, n4 = 0;
 470
 471       if(fDZ)
 472         {
 473           // Define a grid of cell for the gaps between SM
 474           Double_t phimin0 = 0., phimin1 = 0., phimin2 = 0., phimin3 = 0., phimin4 = 0.;
 475           Double_t phimax0 = 0., phimax1 = 0., phimax2 = 0., phimax3 = 0., phimax4 = 0.;
 476           fGeom->GetPhiBoundariesOfSMGap(0,phimin0,phimax0);
 477           fGrid0 = new AliJetGrid(0,95,phimin0,phimax0,-0.7,0.7); // 0.015 x 0.015
 478           fGrid0->SetGridType(0);
 479           fGrid0->SetMatrixIndexes();
 480           fGrid0->SetIndexIJ();
 481           n0 = fGrid0->GetNEntries();
 482           fGeom->GetPhiBoundariesOfSMGap(1,phimin1,phimax1);
 483           fGrid1 = new AliJetGrid(0,95,phimin1,phimax1,-0.7,0.7); // 0.015 x 0.015
 484           fGrid1->SetGridType(0);
 485           fGrid1->SetMatrixIndexes();
 486           fGrid1->SetIndexIJ();
 487           n1 = fGrid1->GetNEntries();
 488           fGeom->GetPhiBoundariesOfSMGap(2,phimin2,phimax2);
 489           fGrid2 = new AliJetGrid(0,95,phimin2,phimax2,-0.7,0.7); // 0.015 x 0.015
 490           fGrid2->SetGridType(0);
 491           fGrid2->SetMatrixIndexes();
 492           fGrid2->SetIndexIJ();
 493           n2 = fGrid2->GetNEntries();
 494           fGeom->GetPhiBoundariesOfSMGap(3,phimin3,phimax3);
 495           fGrid3 = new AliJetGrid(0,95,phimin3,phimax3,-0.7,0.7); // 0.015 x 0.015
 496           fGrid3->SetGridType(0);
 497           fGrid3->SetMatrixIndexes();
 498           fGrid3->SetIndexIJ();
 499           n3 = fGrid3->GetNEntries();
 500           fGeom->GetPhiBoundariesOfSMGap(4,phimin4,phimax4);
 501           fGrid4 = new AliJetGrid(0,95,phimin4,phimax4,-0.7,0.7); // 0.015 x 0.015
 502           fGrid4->SetGridType(0);
 503           fGrid4->SetMatrixIndexes();
 504           fGrid4->SetIndexIJ();
 505           n4 = fGrid4->GetNEntries();
 506
 507           nGaps = n0+n1+n2+n3+n4;
 508
 509         }
 510
 511       for(Int_t nBin = 0; nBin < fNumUnits+nElements+nGaps; nBin++)
 512         {
 513           if(nBin<fNumUnits)
 514             {
 515               fGeom->EtaPhiFromIndex(nBin, eta, phi); // From EMCal geometry
 516               // fEmcalGrid->GetEtaPhiFromIndex2(nBin,phi,eta); // My function from Grid
 517               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 518               Deta = fEmcalGrid->GetDeta(); // Modify with the exact detector values
 519               Dphi = fEmcalGrid->GetDphi(); // Modify with the exact detector values
 520               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 521             }
 522           else {
 523             if(nBin>=fNumUnits && nBin<fNumUnits+nElements){
 524               fTpcGrid->GetEtaPhiFromIndex2(nBin+1-fNumUnits,phi,eta);
 525               phi = ((phi < 0) ? phi + 2. * TMath::Pi() : phi);
 526               Deta = fTpcGrid->GetDeta();
 527               Dphi = fTpcGrid->GetDphi();
 528               new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 529             }
 530             else {
 531               if(fDZ) {
 532                 if(nBin>=fNumUnits+nElements && nBin<fNumUnits+nElements+nGaps){
 533                   if(nBin<fNumUnits+nElements+n0)
 534                     {
 535                       Float_t phi0 = eta = 0.;
 536                       fGrid0->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements),phi0,eta);
 537                       Deta = fGrid0->GetDeta();
 538                       Dphi = fGrid0->GetDphi();
 539                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 540                     }
 541                   else if(nBin>=fNumUnits+nElements+n0 && nBin<fNumUnits+nElements+n0+n1)
 542                     {
 543                       Float_t phi0 = eta = 0.;
 544                       fGrid1->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0),phi0,eta);
 545                       Deta = fGrid1->GetDeta();
 546                       Dphi = fGrid1->GetDphi();
 547                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 548                     }
 549                   else if(nBin>=fNumUnits+nElements+n0+n1 && nBin<fNumUnits+nElements+n0+n1+n2)
 550                     {
 551                       Float_t phi0 = eta = 0.;
 552                       fGrid2->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1),phi0,eta);
 553                       Deta = fGrid2->GetDeta();
 554                       Dphi = fGrid2->GetDphi();
 555                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 556                     }
 557                   else if(nBin>=fNumUnits+nElements+n0+n1+n2 && nBin<fNumUnits+nElements+n0+n1+n2+n3)
 558                     {
 559                       Float_t phi0 = eta = 0.;
 560                       fGrid3->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2),phi0,eta);
 561                       Deta = fGrid3->GetDeta();
 562                       Dphi = fGrid3->GetDphi();
 563                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 564                     }
 565                   else if(nBin>=fNumUnits+nElements+n0+n1+n2+n3 && nBin<fNumUnits+nElements+nGaps)
 566                     {
 567                       Float_t phi0 = eta = 0.;
 568                       fGrid4->GetEtaPhiFromIndex2(nBin+1-(fNumUnits+nElements+n0+n1+n2+n3),phi0,eta);
 569                       Deta = fGrid4->GetDeta();
 570                       Dphi = fGrid4->GetDphi();
 571                       new ((*fUnitArray)[nBin]) AliJetUnitArray(nBin,0,eta,phi0,0.,Deta,Dphi,kTpc,kOutJet,kPtSmaller,kPtSmaller,kBad,0.,-1);
 572                     }
 573                 }
 574               } // end if(fDZ)
 575             } // end else 2
 576           } // end else 1
 577         } // end loop on nBin
 578     } // end grid type == 1
 579   fArrayInitialised = 1;
 580 }
 581
 582
 583