StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StGeant2LcpTreeMaker.cxx
1 // *-- Author : Jan Balewski
2 //
3 // $Id: StGeant2LcpTreeMaker.cxx,v 1.4 2007/07/12 19:24:55 fisyak Exp $
4 
5 #include <math.h>
6 #include <TFile.h>
7 #include <TH2.h>
8 
9 #include "StGeant2LcpTreeMaker.h"
10 
11 #include "TChain.h"
12 #include "TClonesArray.h"
13 
14 #include "StEventInfo.h"
15 #include "StMuDSTMaker/COMMON/StMuEvent.h"
16 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
17 
18 #include "St_DataSet.h"
19 #include "St_DataSetIter.h"
20 #include "tables/St_g2t_event_Table.h"
21 #include "tables/St_particle_Table.h"
22 
23 #include "StMuDSTMaker/COMMON/StMuTrack.h"
24 #include "TMath.h"
25 ClassImp(StGeant2LcpTreeMaker)
26 
27 
28 //________________________________________________
29 //________________________________________________
30 StGeant2LcpTreeMaker::StGeant2LcpTreeMaker(const char* self ,const char* muDstMakerName) : StMaker(self){
31  mMuDstMaker = (StMuDstMaker*)GetMaker(muDstMakerName);
32  assert(mMuDstMaker);
33 
34  runID=101;// default not initialized value
35  treeName="./";
36  SetMaxEta(1.0);
37  SetMinPt(0.4); // GeV/c
38 }
39 
40 
41 //___________________ _____________________________
42 //________________________________________________
43 StGeant2LcpTreeMaker::~StGeant2LcpTreeMaker(){
44  if(tree) tree->Print();
45 
46  // Save all objects in this file
47  hfile->Write();
48 
49 }
50 
51 
52 //________________________________________________
53 //________________________________________________
54 Int_t StGeant2LcpTreeMaker::Init(){
55 
56  return StMaker::Init();
57 }
58 
59 
60 //________________________________________________
61 //________________________________________________
62 void StGeant2LcpTreeMaker::clearLCP(){
63  lcp_eta=lcp_phi=lcp_pt=lcp_e=0;
64  lcp_idhep=0;
65 }
66 
67 //________________________________________________
68 //________________________________________________
69 void StGeant2LcpTreeMaker::InitRunFromMake (int runNumber){
70  if(runID==runNumber) return ;
71  if(runID==888999) return ; // special run number
72 
73 
74  if(runID !=101) {
75  printf("\n\n FATAL %s::InitRunFromMake() runID=%d changed to %d, (not implemented),STOP \n",GetName(),runID,runNumber);
76  assert(runID==runNumber);
77  }
78 
79  runID=runNumber;
80  float Pi=3.14159;
81  treeName+="/G";
82  treeName+=runID;
83  treeName+=".tree.root";
84 
85  printf("%s::Init() \n runID=%d Setup output tree & histo to '%s' , CUTS: maxEta=%.2f minPt/GeV/c=%.2f\n",GetName(),runID,treeName.Data(),C_maxEta,C_minPt);
86 
87  hfile=new TFile(treeName,"RECREATE"," histograms & trees with LCP ");
88  assert(hfile->IsOpen());
89 
90  // Create some histograms and a profile histogram
91  h[0]=0;
92  h[1]=0;
93  h[2]=0;
94  h[3]=0;
95  h[4]=0;
96  h[5] = new TH1F("nPKpi","No. of stable p,pbar,K+/-,pi+/- with |eta|<1 and pT>0.4",50,-0.5,49.5);
97  h[6] =new TH1F("eta","eta of LCP",100,-2.,2.);
98 
99  h[7] = new TH1F("phi"," phi/rad of LCP", 240,-Pi,Pi);
100  h[8] = new TH1F("pT","pT/GeV of LCP",100,0,10.);
101  h[9] = 0;
102  h[10] = 0;
103 
104 
105  if (runID>100) {
106  // Create tree
107  tree = new TTree("G-LCP","LCP selected from Geant data");
108  tree->Branch("id", &eve_id,"id/I");
109  tree->Branch("sub", &eve_sub,"sub/I");
110  tree->Branch("nPKPi",&eve_nPKPi,"nPKPi/I");
111  tree->Branch("pt",&lcp_pt ,"pt/F");
112  tree->Branch("phi",&lcp_phi ,"phi/F");
113  tree->Branch("eta",&lcp_eta ,"eta/F");
114  tree->Branch("e", &lcp_e ,"e/F");
115  } else {
116  printf(" %s::Init(RunNumber=%d) WARN : TTRee NOT created\n",GetName(),runID);
117  tree=0;
118  }
119 return;
120 }
121 
122 //________________________________________________
123 //________________________________________________
124 Int_t StGeant2LcpTreeMaker::Make(){
125 
126  static int kEve=0;
127  kEve++;
128  clearLCP();
129 
130  printf("%s::Make() is called ..........\n",GetName());
131 
132  StMuDst* dst = mMuDstMaker->muDst();
133 
134  StMuEvent* muEve = dst->event();
135  StEventInfo &info=muEve->eventInfo();
136 
137  if(runID>100) InitRunFromMake(info.runId());
138 
139  printf(" eve ID=%d\n",info.id());
140 
141 
142  // Access to geant-tables .......................
143 
144  St_DataSet* Event = GetDataSet("geant");
145  //Event->ls(3);
146  St_DataSetIter geantDstI(Event);
147 
148  // Event generator info ........................
149  St_g2t_event *Pg2t_event=(St_g2t_event *) geantDstI("g2t_event");
150 
151  g2t_event_st *g2t_event1=Pg2t_event->GetTable();
152  printf("nr=%d \n",(int)Pg2t_event->GetNRows());
153  int k1= g2t_event1->eg_label;
154  int k2= g2t_event1->n_event;
155  int k3= g2t_event1->subprocess_id;
156 
157  printf("eg_label=%d n_event=%d subprocess_id=%d\n", k1,k2,k3);
158 
159  assert(info.id()==g2t_event1->n_event);
160  eve_id=info.id();
161  eve_sub=g2t_event1->subprocess_id;
162 
163  // This is an example to access the particle table directly.
164  St_particle *particleTabPtr = (St_particle *) geantDstI("particle");
165 
166  findGeantLcp( particleTabPtr);
167 
168  if(tree)tree->Fill();
169 
170  //if(kEve%500==0)
171  printf("#%s::Make(%d) nPKPi=%d lcp: pt=%f idHep=%d \n",GetName(),kEve,eve_nPKPi, lcp_pt,lcp_idhep);
172 
173 
174  return kStOK;
175 }
176 
177 
178 //================================================
179 //================================================
180 void StGeant2LcpTreeMaker::printTrack( particle_st* part) {
181  cout << " ist = " << part->isthep;
182  cout << " id = " << part->idhep << endl;
183  cout << " px = " << part->phep[0];
184  cout << " py = " << part->phep[1];
185  cout << " pz = " << part->phep[2] << endl;
186  cout << " e = " << part->phep[3];
187  cout << " m = " << part->phep[4] << endl;
188  //cout << " moth1 = " << part->jmohep[0] << endl;
189  //cout << " moth2 = " << part->jmohep[1] << endl;
190 
191 }
192 
193 //================================================
194 //================================================
195 particle_st* StGeant2LcpTreeMaker::findGeantLcp( St_particle *tab) {
196  assert(tab);
197  particle_st* part = tab->GetTable();
198  printf("tab dim=%d\n", (int)tab->GetNRows());
199  int i;
200  float maxPt=0.4;
201  float maxEta=1.;
202  float etaI=99;
203  int iMax=-1;
204  int ns=0, np=0,n2=0,n3=0,npkpi=0;
205  for(i=0;i<tab->GetNRows();i++) {
206  if(part[i].isthep!=1) continue;
207  ns++;
208  int id=part[i].idhep;
209  int idM=id>0 ? id : -id ;
210  if(idM!=211 && idM!=321 && idM!=2212) continue;
211  // accept only p,pbar, pi+/-, K+/-
212  np++;
213 
214  float px=part[i].phep[0];
215  float py=part[i].phep[1];
216  float pz=part[i].phep[2];
217 
218  float pt=sqrt(px*px+py*py);
219  if(pt<0.001) continue; // to allow for math
220  float eta= TMath::ASinH(pz/pt);
221 
222 #if 0
223  {
224  float phi=atan2(part[i].phep[1],part[i].phep[0]);
225  printf("i=%d pT=%f phi/deg=%f eta=%f pz=%f idHep=%d\n",i,pt,phi/3.1416*180,eta,pz,id);
226  }
227 #endif
228  //printf("px=%f py=%f pt=%f\n",px,py,pt);
229  if(pt>0.4 && fabs(eta)<1) npkpi++;
230  if(pt<maxPt) continue;
231  n2++;
232 
233  if(fabs(eta)>maxEta) continue;
234  // track is in eta range
235 
236  n3++;
237  maxPt=pt;
238  etaI=eta;
239  iMax=i;
240  }
241 
242  printf("#total %3d stable particles, np=%3d n2=%3d n3=%3d maxPt=%f iMax=%d eta=%f npkpi=%d\n",ns,np,n2,n3,maxPt,iMax,etaI,npkpi);
243 
244  eve_nPKPi=npkpi;
245  h[5]->Fill( eve_nPKPi);
246 
247  if(iMax>=0) {
248  printTrack(part+iMax);
249  lcp_eta=etaI;
250  lcp_phi=atan2(part[iMax].phep[1],part[iMax].phep[0]);
251  lcp_pt=maxPt;
252  lcp_idhep=part[iMax].idhep;
253  lcp_e=part[iMax].phep[3];
254  h[6]->Fill(lcp_eta);
255  h[7]->Fill(lcp_phi);
256  h[8]->Fill(lcp_pt);
257 
258  return part+iMax;
259  }
260 
261  return 0;
262 
263 }
264 
265 
266 
267 // $Log: StGeant2LcpTreeMaker.cxx,v $
268 // Revision 1.4 2007/07/12 19:24:55 fisyak
269 // Add includes for TMath for ROOT 5.16
270 //
271 // Revision 1.3 2004/02/24 21:33:23 balewski
272 // *** empty log message ***
273 //
274 // Revision 1.2 2004/01/26 22:49:37 perev
275 // WarnOff
276 //
277 // Revision 1.1 2004/01/06 17:25:26 balewski
278 // get LCP from Geant info
279 //
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StGeant2LcpTreeMaker
This commented block at the top ...
Definition: StGeant2LcpTreeMaker.h:29
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMaker::GetName
virtual const char * GetName() const
special overload
Definition: StMaker.cxx:237
kStOK
Definition: Stypes.h:40
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5