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