Classed moved form libJETANMC to libJETAN.
[u/mrichter/AliRoot.git] / JETAN / AliJetESDmcReader.cxx
CommitLineData
1994c1b0 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// Jet ESD Reader
17// ESD reader for jet analysis (it reads the esd and the MC trees)
18// Author: Mercedes Lopez Noriega (mercedes.lopez.noriega@cern.ch)
19
20#include <Riostream.h>
21#include <TSystem.h>
22#include <TLorentzVector.h>
23#include <TVector3.h>
24#include <TPDGCode.h>
25#include <TParticle.h>
26#include <TParticlePDG.h>
27#include <TMath.h>
28#include "AliJetESDmcReader.h"
29#include "AliJetESDReaderHeader.h"
1b307662 30#include "AliESDEvent.h"
1994c1b0 31#include "AliESDtrack.h"
0ffa8579 32#include "AliHeader.h"
1994c1b0 33
34ClassImp(AliJetESDmcReader)
35
1b7d5d7e 36AliJetESDmcReader::AliJetESDmcReader():
b45b0c92 37 AliJetESDReader(),
38 fChainMC(0x0),
39 fAliHeader(0)
1994c1b0 40{
41 // Constructor
1994c1b0 42}
43
44//____________________________________________________________________________
45
46AliJetESDmcReader::~AliJetESDmcReader()
47{
48 // Destructor
b45b0c92 49 delete fChainMC;
1994c1b0 50}
51
52//____________________________________________________________________________
53
54void AliJetESDmcReader::OpenInputFiles()
55
56{
57 // chain for the ESDs
58 fChain = new TChain("esdTree");
59 fChainMC = new TChain("mcStackTree");
60
61 // get directory and pattern name from the header
62 const char* dirName=fReaderHeader->GetDirectory();
63 const char* pattern=fReaderHeader->GetPattern();
64
65 // Add files matching patters to the chain
66 void *dir = gSystem->OpenDirectory(dirName);
67 const char *name = 0x0;
b45b0c92 68 int nesd = ((AliJetESDReaderHeader*)fReaderHeader)->GetNesd();
1994c1b0 69 int a = 0;
70 while ((name = gSystem->GetDirEntry(dir))){
71 if (a>=nesd) continue;
72 if (strstr(name,pattern)){
73 printf("Adding %s\n",name);
74 char path[256];
75 sprintf(path,"%s/%s",dirName,name);
76 fChain->AddFile(path,-1);
77 fChainMC->AddFile(path,-1);
78 a++;
79 }
80 }
81 printf("%d ESDs added\n",a);
82
83 gSystem->FreeDirectory(dir);
84 fChain->SetBranchAddress("ESD", &fESD);
85 fChainMC->SetBranchAddress("Header", &fAliHeader);
86 fChainMC->SetBranchAddress("Stack", &fArrayMC);
87
88 int nMax = fChain->GetEntries();
89 printf("\nTotal number of events in chain= %d\n",nMax);
90
91 // set number of events in header
92 if (fReaderHeader->GetLastEvent() == -1)
93 fReaderHeader->SetLastEvent(nMax);
94 else {
95 Int_t nUsr = fReaderHeader->GetLastEvent();
96 fReaderHeader->SetLastEvent(TMath::Min(nMax,nUsr));
97 }
98}
99
100//____________________________________________________________________________
101
0ffa8579 102Bool_t AliJetESDmcReader::FillMomentumArray(Int_t event)
1994c1b0 103{
104 // Fill momentum array
105
106 Int_t goodTrack = 0;
107 Int_t nt = 0;
108 Int_t pdgCode = 0;
109 Float_t pt, eta;
110 Float_t energyMC, pxMC, pyMC, pzMC, ptMC; // Monte Carlo
111 TVector3 p3;
112
113 // clear momentum array
114 ClearArray();
115
116 // get event from chain
117 fChain->GetEntry(event);
118 fChainMC->GetEntry(event);
119
120 // get number of tracks in event (for the loop)
121 nt = fESD->GetNumberOfTracks();
122
123 // temporary storage of signal and pt cut flag
124 Int_t* sflag = new Int_t[nt];
125 Int_t* cflag = new Int_t[nt];
126
127 // get cuts set by user
128 Float_t ptMin = fReaderHeader->GetPtCut();
129 Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
130 Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
131
132 //loop over tracks
133 for (Int_t it = 0; it < nt; it++) {
134 AliESDtrack *track = fESD->GetTrack(it);
135 UInt_t status = track->GetStatus();
136
137 Double_t mom[3];
138 track->GetPxPyPz(mom);
139 p3.SetXYZ(mom[0],mom[1],mom[2]);
140 pt = p3.Pt();
141 if (((status & AliESDtrack::kITSrefit) == 0) ||
142 ((status & AliESDtrack::kTPCrefit) == 0)) continue; // quality check
143 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadSignalOnly()
144 && TMath::Abs(track->GetLabel()) > 10000) continue; // quality check
145 if (((AliJetESDReaderHeader*) fReaderHeader)->ReadBkgdOnly()
146 && TMath::Abs(track->GetLabel()) < 10000) continue; // quality check
147 eta = p3.Eta();
148 if ( (eta > etaMax) || (eta < etaMin)) continue; // checking eta cut
149
150 new ((*fMomentumArray)[goodTrack]) TLorentzVector(p3,p3.Mag());
151
152 sflag[goodTrack]=0;
153 if (TMath::Abs(track->GetLabel()) < 10000) sflag[goodTrack]=1;
154 cflag[goodTrack]=0;
155 if (pt > ptMin) cflag[goodTrack]=1; // pt cut
156
157 goodTrack++;
158
159 // Monte Carlo information
160 Int_t label = TMath::Abs(track->GetLabel());
161 TClonesArray &arrayMC = *fArrayMC;
162 TParticle *part = (TParticle*)arrayMC[label]; //particle
163 ptMC = part->Pt();
164 pdgCode = part->GetPdgCode();
165 energyMC = part->Energy();
166 pxMC = part->Px(); pyMC = part->Py(); pzMC = part->Pz();
167 }
168 // set the signal flags
169 fSignalFlag.Set(goodTrack,sflag);
170 fCutFlag.Set(goodTrack,cflag);
0ffa8579 171 return kTRUE;
172
1994c1b0 173}
174
175