StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StJetMaker2012.cxx
1 //******************************
2 //2012 version of Jet Maker
3 //with off-axis cone underlying
4 //-event study included the unde-
5 //rlying event will be saved in
6 //the tree called "ue". The Mak-
7 //-er is derivatived from the
8 //older 2009 version.
9 //
10 //Author: Zilong Chang
11 //user name: zchang
12 //Cyclotron Institute
13 //Texas A&M University
14 //******************************
15 #include "StThreeVector.hh"
16 #include "TFile.h"
17 #include "TTree.h"
18 #include "StAnaPars.h"
19 #include "StjTPCMuDst.h"
20 //#include "StjTPCRandomMuDst.h"
21 #include "StjBEMCMuDst.h"
22 #include "StjEEMCMuDst.h"
23 #include "StjFMSMuDst.h"
24 #include "StjFCSMuDst.h"
25 #include "StjMCMuDst.h"
26 #include "StjTPCNull.h"
27 #include "StjBEMCNull.h"
28 #include "StjEEMCNull.h"
29 #include "StjFMSNull.h"
30 #include "StjAbstractTowerEnergyCorrectionForTracks.h"
31 #include "StjeTrackListToStMuTrackFourVecList.h"
32 #include "StjeTowerEnergyListToStMuTrackFourVecList.h"
33 #include "StjMCParticleToStMuTrackFourVec.h"
34 #include "StJetFinder/StJetFinder.h"
35 #include "StSpinPool/StJetEvent/StJetEventTypes.h"
36 #include "StMuTrackFourVec.h"
37 #include "StMuTrackEmu.h"
38 #include "StMuTowerEmu.h"
39 #include "StMcTrackEmu.h"
40 
41 #include "StarGenerator/StarGenEventReader/StarGenEventReader.h"
42 
43 #include "StSpinPool/StUeEvent/StUeOffAxisConesEvent.h"
44 #include "StSpinPool/StUeEvent/StUeVertex.h"
45 #include "StSpinPool/StUeEvent/StUeOffAxisConesJet.h"
46 #include "StSpinPool/StUeEvent/StUeOffAxisCones.h"
47 
48 #include "StJetMaker2012.h"
49 
50 ClassImp(StJetMaker2012);
51 
52 void StJetMaker2012::Clear(Option_t* option)
53 {
54  StJetMaker2009::Clear(option);
55  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ++ueBranch){
56  StJetUeBranch* jetuebranch = mJetUeBranches[ueBranch];
57  for (size_t iBranch = 0; iBranch < mJetBranches.size(); ++iBranch) {
58  (jetuebranch->eventUe).at(iBranch)->Clear(option);
59  // jetbranch->event->Clear(option);
60  }
61  }
62  StJetMaker2009::Clear(option);
63 }
64 
65 int StJetMaker2012::Init()
66 {
67  assert(!mFileNameUe.IsNull());
68  mFileUe = TFile::Open(mFileNameUe,"recreate");
69  assert(mFileUe);
70 
71  mTreeUe = new TTree("ue","ueTree");
72 
73  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ++ueBranch){
74  StJetUeBranch* jetuebranch = mJetUeBranches[ueBranch];
75  for (size_t iBranch = 0; iBranch < mJetBranches.size(); ++iBranch) {
76  StJetBranch* jetbranch = mJetBranches[iBranch];
77  TString branchname(jetbranch->name);
78  branchname.Append(jetuebranch->name);
79  mTreeUe->Branch(branchname, "StUeOffAxisConesEvent",&((jetuebranch->eventUe)[iBranch]));
80  LOG_INFO<<jetbranch->name<<" and "<<jetuebranch->name<<" initialized!"<<endm;
81  }
82  }
83 
84  mTreeUe->BranchRef();
85 
86  return StJetMaker2009::Init();
87  //return kStOk;
88 }
89 
90 int StJetMaker2012::Make()
91 {
92  // Loop over jet branches
93  for (size_t iBranch = 0; iBranch < mJetBranches.size(); ++iBranch) {
94  StJetBranch* jetbranch = mJetBranches[iBranch];
95 
96  // Fill header
97  jetbranch->event->mRunId = GetRunNumber();
98  jetbranch->event->mEventId = GetEventNumber();
99  jetbranch->event->mDatime = GetDateTime();
100 
101  if (jetbranch->anapars->useTpc) {
102 // StjTPCRandomMuDst tpc((StMuDstMaker*)0,
103 // jetbranch->anapars->randomSelectorProb,
104 // jetbranch->anapars->randomSelectorAt,
105 // jetbranch->anapars->randomSelectorSeed);
106  StjTPCMuDst tpc;
107 
108  // Save vertex index
109  int savedVertexIndex = tpc.currentVertexIndex();
110 
111  // Keep track of number of good vertices
112  int nvertices = 0;
113 
114  int totNumVertices = tpc.numberOfVertices();
115  // store the total number of available vertices if set by anapars
116  if(jetbranch->anapars->storeOnlyDefaultVertex && totNumVertices > 1) {
117  totNumVertices = 1;
118  }
119 
120  // Vertex loop
121  for (int iVertex = 0; iVertex < totNumVertices; ++iVertex) {
122  tpc.setVertexIndex(iVertex);
123 
124  // Get TPC vertex and tracks
125  StjPrimaryVertex vertex = tpc.getVertex();
126 
127  // Skip pile-up vertex
128  if (vertex.ranking() <= 0) continue;
129 
130  // Found good vertex
131  ++nvertices;
132  StjTrackList trackList = tpc.getTrackList();
133  if (jetbranch->anapars->changeTracks) trackList = (*jetbranch->anapars->changeTracks)(trackList);
134  trackList = jetbranch->anapars->tpcCuts()(trackList);
135 
136  // Get BEMC towers
137  StjTowerEnergyList bemcEnergyList;
138 
139  if (jetbranch->anapars->useBemc) {
140  StjBEMCMuDst bemc;
141  bemcEnergyList = bemc.getEnergyList();
142  if (jetbranch->anapars->changeTowers) (*jetbranch->anapars->changeTowers)(bemcEnergyList);
143  bemcEnergyList = jetbranch->anapars->bemcCuts()(bemcEnergyList);
144  }
145 
146  // Get EEMC towers
147  StjTowerEnergyList eemcEnergyList;
148 
149  if (jetbranch->anapars->useEemc) {
150  StjEEMCMuDst eemc;
151  eemcEnergyList = eemc.getEnergyList();
152  if (jetbranch->anapars->changeTowers) (*jetbranch->anapars->changeTowers)(eemcEnergyList);
153  eemcEnergyList = jetbranch->anapars->eemcCuts()(eemcEnergyList);
154  }
155 
156  // Get FMS towers
157  StjTowerEnergyList fmsEnergyList;
158  if (jetbranch->anapars->useFms) {
159  StjFMSMuDst fms;
160  fmsEnergyList = fms.getEnergyList();
161  }
162 
163  // Get FCS Towers
164 
165  StjTowerEnergyList fcsECalEnergyList; //FCS ECAL
166  if(jetbranch->anapars->useFcsECal){
167  StjFCSMuDst fcsECal;
168  fcsECal.useECAL();
169  fcsECalEnergyList = fcsECal.getEnergyList();
170  fcsECalEnergyList = jetbranch->anapars->FCSEcalemcCuts()(fcsECalEnergyList);
171 
172  }
173  StjTowerEnergyList fcsHCalEnergyList; //FCS HCAL
174  if(jetbranch->anapars->useFcsHCal){
175  StjFCSMuDst fcsHCal;
176  fcsHCal.useHCAL();
177  fcsHCalEnergyList = fcsHCal.getEnergyList();
178  fcsHCalEnergyList = jetbranch->anapars->FCSHcalhcCuts()(fcsHCalEnergyList);
179 
180  }
181 
182  // Merge BEMC and EEMC towers
183  StjTowerEnergyList energyList;
184 
185  copy(bemcEnergyList.begin(),bemcEnergyList.end(),back_inserter(energyList));
186  copy(eemcEnergyList.begin(),eemcEnergyList.end(),back_inserter(energyList));
187  copy(fmsEnergyList.begin(),fmsEnergyList.end(),back_inserter(energyList));
188 
189  // Apply hadronic correction to towers
190  energyList = jetbranch->anapars->correctTowerEnergyForTracks()(energyList,trackList);
191 
192  //Merge FCS towers with general energyList
193  copy(fcsECalEnergyList.begin(),fcsECalEnergyList.end(),back_inserter(energyList));
194  copy(fcsHCalEnergyList.begin(),fcsHCalEnergyList.end(),back_inserter(energyList));
195 
196  // Currently no hadronic corrections to towers for FCS
197 
198  // Convert tracks and towers to Lorentz vectors
199  FourList tpc4pList = StjeTrackListToStMuTrackFourVecList()(trackList);
200  FourList energy4pList = StjeTowerEnergyListToStMuTrackFourVecList()(energyList);
201 
202  // Merge tracks and towers
203  StProtoJet::FourVecList particles; // vector<const AbstractFourVec*>
204 
205  copy(tpc4pList.begin(),tpc4pList.end(),back_inserter(particles));
206  copy(energy4pList.begin(),energy4pList.end(), back_inserter(particles));
207 
208  // Run jet finder
209  StJetFinder::JetList protojets; // list<StProtoJet>
210  jetbranch->jetfinder->findJets(protojets,particles);
211 
212  // Filter jets
213  protojets = jetbranch->anapars->jetCuts()(protojets);
214 
215  // Add vertex
216  StJetVertex* jetvertex = jetbranch->event->newVertex();
217  copyVertex(vertex,jetvertex);
218  //Add UE vertex
219  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ueBranch++){
220  ((mJetUeBranches[ueBranch]->eventUe).at(iBranch))->setEventId(GetEventNumber());
221  StUeVertex *uevertex = (mJetUeBranches[ueBranch]->eventUe).at(iBranch)->newVertex();
222  uevertex->setPosition(vertex.position());
223  uevertex->setRanking(vertex.ranking());
224  }
225 
226  // Add jets
227  for (StJetFinder::JetList::const_iterator iProtoJet = protojets.begin(); iProtoJet != protojets.end(); ++iProtoJet){
228  addJet(*iProtoJet,jetbranch->event,jetvertex);
229  //UE jets
230  StJetCandidate* currentjet = jetbranch->event->jet(jetbranch->event->numberOfJets()-1);
231  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ueBranch++){
232  StJetUeBranch *jetuebranch = mJetUeBranches.at(ueBranch);
233  double radius = jetuebranch->uePars->coneRadius();
234  double uedensity = addJetUe(particles, (jetuebranch->eventUe).at(iBranch), *iProtoJet, radius);
235  currentjet->setUeDensity(jetuebranch->name, uedensity);
236  }
237  }
238  // Clean up particles
239  for (StProtoJet::FourVecList::const_iterator i = particles.begin(); i != particles.end(); ++i)
240  delete *i;
241  } // End vertex loop
242 
243  // No good vertex was found. Use (0,0,0) for vertex and EMC-only jets.
244  // Useful for beam-gas background studies.
245  if (!nvertices) {
246  // Get BEMC towers
247  StjTowerEnergyList bemcEnergyList;
248 
249  if (jetbranch->anapars->useBemc) {
250  StjBEMCMuDst bemc;
251  bemc.setVertex(0,0,0);
252  bemcEnergyList = bemc.getEnergyList();
253  if (jetbranch->anapars->changeTowers) (*jetbranch->anapars->changeTowers)(bemcEnergyList);
254  bemcEnergyList = jetbranch->anapars->bemcCuts()(bemcEnergyList);
255  }
256 
257  // Get EEMC towers
258  StjTowerEnergyList eemcEnergyList_temp, eemcEnergyList;
259 
260  if (jetbranch->anapars->useEemc) {
261  StjEEMCMuDst eemc;
262  eemcEnergyList = eemc.getEnergyList();
263  if (jetbranch->anapars->changeTowers) (*jetbranch->anapars->changeTowers)(eemcEnergyList);
264  eemcEnergyList = jetbranch->anapars->eemcCuts()(eemcEnergyList);
265  }
266 
267  // Get FMS towers
268  StjTowerEnergyList fmsEnergyList;
269  if (jetbranch->anapars->useFms) {
270  StjFMSMuDst fms;
271  fmsEnergyList = fms.getEnergyList();
272  }
273 
274  // Get FCS Towers
275 
276  StjTowerEnergyList fcsECalEnergyList; //FCS ECAL
277  if(jetbranch->anapars->useFcsECal){
278  StjFCSMuDst fcsECal;
279  fcsECal.useECAL();
280  fcsECalEnergyList = fcsECal.getEnergyList();
281  fcsECalEnergyList = jetbranch->anapars->FCSEcalemcCuts()(fcsECalEnergyList);
282 
283  }
284  StjTowerEnergyList fcsHCalEnergyList; //FCS HCAL
285  if(jetbranch->anapars->useFcsHCal){
286  StjFCSMuDst fcsHCal;
287  fcsHCal.useHCAL();
288  fcsHCalEnergyList = fcsHCal.getEnergyList();
289  fcsHCalEnergyList = jetbranch->anapars->FCSHcalhcCuts()(fcsHCalEnergyList);
290 
291  }
292 
293  // Merge BEMC and EEMC towers
294  StjTowerEnergyList energyList;
295 
296  copy(bemcEnergyList.begin(),bemcEnergyList.end(),back_inserter(energyList));
297  copy(eemcEnergyList.begin(),eemcEnergyList.end(),back_inserter(energyList));
298  copy(fmsEnergyList.begin(),fmsEnergyList.end(),back_inserter(energyList));
299  copy(fcsECalEnergyList.begin(),fcsECalEnergyList.end(),back_inserter(energyList));
300  copy(fcsHCalEnergyList.begin(),fcsHCalEnergyList.end(),back_inserter(energyList));
301 
302  // Convert towers to Lorentz vectors
303  FourList energy4pList = StjeTowerEnergyListToStMuTrackFourVecList()(energyList);
304 
305  // Merge tracks and towers
306  StProtoJet::FourVecList particles; // vector<const AbstractFourVec*>
307 
308  copy(energy4pList.begin(),energy4pList.end(),back_inserter(particles));
309 
310  // Run jet finder
311  StJetFinder::JetList protojets; // list<StProtoJet>
312  jetbranch->jetfinder->findJets(protojets,particles);
313 
314  // Filter jets
315  protojets = jetbranch->anapars->jetCuts()(protojets);
316 
317  // Add dummy vertex (0,0,0)
318  StJetVertex* jetvertex = jetbranch->event->newVertex();
319  jetvertex->mPosition.SetXYZ(0,0,0);
320 
321  // Add jets
322  for (StJetFinder::JetList::const_iterator iProtoJet = protojets.begin(); iProtoJet != protojets.end(); ++iProtoJet)
323  addJet(*iProtoJet,jetbranch->event,jetvertex);
324 
325  // Clean up particles
326  for (StProtoJet::FourVecList::const_iterator i = particles.begin(); i != particles.end(); ++i)
327  delete *i;
328  }
329 
330  // Restore vertex index
331  tpc.setVertexIndex(savedVertexIndex);
332 
333  } // End useTpc
334 
335 
336  if (jetbranch->anapars->useMonteCarlo) {
337 
338  StjMCMuDst mc(this);
339 
340  if(GetMaker("genEvent")){ // To be used with .genevent.root for MC particle from StarGenEventReader
341  eventreader = (StarGenEventReader*)GetMaker("genEvent");
342  sEvent = eventreader->Event();
343  mc.setGenEvent(sEvent); // In case no geant on chain
344  }
345 
346  StjPrimaryVertex mcvertex = mc.getMCVertex();
347  StjMCParticleList mcparticles = jetbranch->anapars->mcCuts()(mc.getMCParticleList());
348  StProtoJet::FourVecList particles; // vector<const AbstractFourVec*>
349  transform(mcparticles.begin(),mcparticles.end(),back_inserter(particles),StjMCParticleToStMuTrackFourVec());
350 
351  // Run jet finder
352  StJetFinder::JetList protojets; // list<StProtoJet>
353  jetbranch->jetfinder->findJets(protojets,particles);
354 
355  // Filter jets
356  protojets = jetbranch->anapars->jetCuts()(protojets);
357 
358  // Add vertex
359  StJetVertex* jetvertex = jetbranch->event->newVertex();
360  copyVertex(mcvertex,jetvertex);
361  //Add UE vertex
362  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ueBranch++){
363  ((mJetUeBranches[ueBranch]->eventUe).at(iBranch))->setEventId(GetEventNumber());
364  StUeVertex *uevertex = (mJetUeBranches[ueBranch]->eventUe).at(iBranch)->newVertex();
365  uevertex->setPosition(mcvertex.position());
366  uevertex->setRanking(mcvertex.ranking());
367  }
368 
369  // Add jets
370  for (StJetFinder::JetList::const_iterator iProtoJet = protojets.begin(); iProtoJet != protojets.end(); ++iProtoJet){
371  addJet(*iProtoJet,jetbranch->event,jetvertex);
372  //add UE jets
373 
374  StJetCandidate* currentjet = jetbranch->event->jet(jetbranch->event->numberOfJets()-1);//pointer to the current jet
375  for(size_t ueBranch = 0; ueBranch < mJetUeBranches.size(); ueBranch++){
376  StJetUeBranch *jetuebranch = mJetUeBranches.at(ueBranch);
377  double radius = jetuebranch->uePars->coneRadius();
378  double uedensity = addJetUe(particles, (jetuebranch->eventUe).at(iBranch), *iProtoJet, radius);
379  currentjet->setUeDensity(jetuebranch->name, uedensity);
380  }
381  }
382 
383  // Clean up particles
384  for (StProtoJet::FourVecList::const_iterator i = particles.begin(); i != particles.end(); ++i)
385  delete *i;
386  } // End useMonteCarlo
387 
388  } // End loop over jet branches
389  mTree->Fill();
390 
391  mTreeUe->Fill();
392 
393  return kStOk;
394 }
395 
396 int StJetMaker2012::Finish()
397 {
398  mFileUe->Write();
399  mFileUe->Close();
400 
401  return StJetMaker2009::Finish();
402 }
403 
404 // Setters
405 
406 void StJetMaker2012::addUeBranch(const char* name, StOffAxisConesPars* pars)
407 {
408  StJetUeBranch *ue = new StJetUeBranch(name, pars);
409  for(size_t iBranch = 0; iBranch < mJetBranches.size(); iBranch++)
410  (ue->eventUe).push_back(new StUeOffAxisConesEvent);
411  mJetUeBranches.push_back(ue);
412 
413 }
414 
415 void StJetMaker2012::setJetFileUe(const char* filename)
416 {
417  mFileNameUe = filename;
418 }
419 TTree* StJetMaker2012::treeUe()
420 {
421  return mTreeUe;
422 }
423 
424 StUeOffAxisConesEvent* StJetMaker2012::eventUe(const char* branchname)
425 {
426  TBranch* branch = mTreeUe->GetBranch(branchname);
427  if (branch) return *(StUeOffAxisConesEvent**)branch->GetAddress();
428  return 0;
429 }
430 
431 double StJetMaker2012::addJetUe(StProtoJet::FourVecList particles, StUeOffAxisConesEvent *ueEvent, const StProtoJet &jet, double radius)
432 {
433 
434  double jpt = jet.pt();
435  double jeta = jet.eta();
436  double jphi = jet.phi();
437  double je = jet.e();
438  double jarea = jet.area();
439 
440  const int Ncone = 2;
441  StUeOffAxisCones *cones[Ncone];
442  const double PI = 3.14159;
443 
444  double pt[Ncone];
445  int number[Ncone];
446  double cone_eta[Ncone];
447  double cone_phi[Ncone];
448  for(int ii = 0; ii < Ncone; ii++){
449  cone_eta[ii] = jeta;
450  cone_phi[ii] = jphi + (PI/2.)*(2*ii-1);
451  cones[ii] = ueEvent->newUeOffAxisCones();
452  pt[ii] = 0.;
453  number[ii] = 0;
454  }
455 
456  for(size_t ii = 0; ii < particles.size(); ii++){
457  StMuTrackFourVec * vect = (StMuTrackFourVec*)particles.at(ii);
458  StThreeVector<double> mom(vect->px(), vect->py(), vect->pz());
459  double p_pt = mom.perp();
460  double p_eta = vect->eta();
461  double p_phi = vect->phi();
462 
463  int index = -1;
464  double dR = DeltaR(p_eta, p_phi, cone_eta[0], cone_phi[0]);
465  if(dR < radius){
466  pt[0] += p_pt;
467  number[0]++;
468  index = 0;
469  }else{
470  dR = DeltaR(p_eta, p_phi, cone_eta[1], cone_phi[1]);
471  if(dR < radius){
472  pt[1] += p_pt;
473  number[1]++;
474  index = 1;
475  }
476  }
477  if(index == 0 || index == 1){
478  //track
479  if (const StMuTrackEmu* t = vect->track()) {
480  StJetTrack* track = ueEvent->newTrack();
481  copyTrack(t,track);
482  cones[index]->addTrack(track);
483  }
484  //tower
485  if (const StMuTowerEmu* t = vect->tower()) {
486  StJetTower* tower = ueEvent->newTower();
487  //copyTower(t, vertex, tower);
488  copyTower(t, tower);
489  cones[index]->addTower(tower);
490  }
491  //particle
492  if (const StMcTrackEmu* t = vect->mctrack()) {
493  StJetParticle* part = ueEvent->newParticle();
494  copyParticle(t,part);
495  cones[index]->addParticle(part);
496  }
497  }
498  }
499 
500  double density = 0.;
501  for(int ii = 0; ii < Ncone; ii++){
502  cones[ii]->setPt(pt[ii]);
503  cones[ii]->setEtaPhi(cone_eta[ii], cone_phi[ii]);
504  cones[ii]->setRadius(radius);
505  cones[ii]->setMult(number[ii]);
506  density += pt[ii]/(PI*radius*radius);
507  }
508  density = density/(Ncone+0.);
509 
510  StUeOffAxisConesJet *jcones = ueEvent->newUeOffAxisConesJet();
511  jcones->setPt(jpt);
512  jcones->setEta(jeta);
513  jcones->setPhi(jphi);
514  jcones->setE(je);
515  jcones->setArea(jarea);
516  jcones->setDensity(density);
517  jcones->addCone(cones[0]);
518  jcones->addCone(cones[1]);
519  StUeVertex *ueVertex = ueEvent->lastVertex();
520  ueVertex->addUeJet(jcones);
521 
522  return density;
523 }
524 double StJetMaker2012::DeltaR(double etaA, double phiA, double etaB, double phiB)
525 {
526  const double PI = 3.14159;
527 
528  double deta = etaA - etaB;
529  double dphi = phiA - phiB;
530 
531  if(dphi > PI) dphi -= 2.*PI;
532  if(dphi < -1.*PI) dphi += 2.*PI;
533 
534  double dR = TMath::Sqrt(deta*deta+dphi*dphi);
535  return dR;
536 }
StarGenerator::Event
StarGenEvent * Event()
Retrieves the event record.
Definition: StarGenerator.h:157
StMuTrackFourVec::phi
double phi() const
angles
Definition: StMuTrackFourVec.h:50
StJetMaker2012::Clear
void Clear(Option_t *option="")
User defined functions.
Definition: StJetMaker2012.cxx:52
StMaker::GetRunNumber
virtual Int_t GetRunNumber() const
Returns the current RunNumber.
Definition: StMaker.cxx:1054
StJetMaker2009::Clear
void Clear(Option_t *option="")
User defined functions.
Definition: StJetMaker2009.cxx:37
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5