- Clean-up of data members and methods in base classes and derived classes
[u/mrichter/AliRoot.git] / JETAN / AliJetESDmcReader.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 // 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 #include "AliHeader.h"
33
34 ClassImp(AliJetESDmcReader)
35
36 AliJetESDmcReader::AliJetESDmcReader():
37     AliJetESDReader(),
38     fChainMC(0x0),
39     fAliHeader(0)
40 {
41   // Constructor    
42 }
43
44 //____________________________________________________________________________
45
46 AliJetESDmcReader::~AliJetESDmcReader()
47 {
48   // Destructor
49     delete fChainMC;
50 }
51
52 //____________________________________________________________________________
53
54 void 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;
68   int nesd = ((AliJetESDReaderHeader*)fReaderHeader)->GetNesd();
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
102 Bool_t AliJetESDmcReader::FillMomentumArray(Int_t event)
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);
171   return kTRUE;
172   
173 }
174
175