StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFwdQAMaker.cxx
1 #include "StFwdTrackMaker/StFwdQAMaker.h"
2 #include "StFwdQAMaker.h"
3 #include "St_base/StMessMgr.h"
4 #include "StBFChain/StBFChain.h"
5 #include "StFwdTrackMaker/StFwdTrackMaker.h"
6 
7 #include "StFwdTrackMaker/include/Tracker/FwdTracker.h"
8 #include "StFwdTrackMaker/include/Tracker/ObjExporter.h"
9 // StEvent includes
10 #include "StEvent/StBTofCollection.h"
11 #include "StEvent/StBTofHeader.h"
12 #include "StEvent/StEvent.h"
13 #include "StEvent/StFttCluster.h"
14 #include "StEvent/StFttCollection.h"
15 #include "StEvent/StFcsCluster.h"
16 #include "StEvent/StFcsCollection.h"
17 #include "StFcsDbMaker/StFcsDb.h"
18 #include "StRoot/StEpdUtil/StEpdGeom.h"
19 #include "StEvent/StFwdTrackCollection.h"
20 #include "StEvent/StFwdTrack.h"
21 
22 
23 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
24 #include "StMuDSTMaker/COMMON/StMuDst.h"
25 #include "StMuDSTMaker/COMMON/StMuEvent.h"
26 #include "StMuDSTMaker/COMMON/StMuFstCollection.h"
27 #include "StMuDSTMaker/COMMON/StMuFstHit.h"
28 #include "StMuDSTMaker/COMMON/StMuPrimaryVertex.h"
29 #include "StMuDSTMaker/COMMON/StMuFwdTrack.h"
30 #include "StMuDSTMaker/COMMON/StMuFwdTrackCollection.h"
31 #include "StMuDSTMaker/COMMON/StMuFcsCollection.h"
32 #include "StMuDSTMaker/COMMON/StMuFcsCluster.h"
33 #include "StMuDSTMaker/COMMON/StMuFcsHit.h"
34 #include "StMuDSTMaker/COMMON/StMuFttCluster.h"
35 #include "StMuDSTMaker/COMMON/StMuFttPoint.h"
36 #include "StMuDSTMaker/COMMON/StMuMcTrack.h"
37 #include "StMuDSTMaker/COMMON/StMuFstHit.h"
38 
40 void FwdQATreeData::clear(){
41  header.clear();
42  mcTracks.reset();
43  fttPoints.reset();
44  fttClusters.reset();
45  fstPoints.reset();
46  reco.reset();
47  seeds.reset();
48  wcal.reset();
49  hcal.reset();
50  wcalHits.reset();
51  hcalHits.reset();
52  epdHits.reset();
53 }
54 
55 FcsClusterWithStarXYZ::FcsClusterWithStarXYZ( StMuFcsCluster *clu, StFcsDb *fcsDb ) {
56  if ( nullptr == clu ) return;
57  StThreeVectorD xyz = fcsDb->getStarXYZfromColumnRow(clu->detectorId(),clu->x(),clu->y());
58  float detOffset = 0.0;
59  if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){
60  detOffset = 715.0; // cm from IP
61  } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){
62  detOffset = 807.0; // cm from IP
63  }
64  mXYZ.SetXYZ( xyz.x(), xyz.y(), xyz.z() + detOffset );
65  mClu = clu;
66 }
67 
68 FcsHitWithStarXYZ::FcsHitWithStarXYZ( StMuFcsHit *hit, StFcsDb *fcsDb ) {
69  if ( nullptr == hit ) return;
70  StThreeVectorD xyz = fcsDb->getStarXYZ(hit->detectorId(),hit->id());
71  float detOffset = 0.0;
72  if ( hit->detectorId() == kFcsEcalNorthDetId || hit->detectorId() == kFcsEcalSouthDetId ){
73  detOffset = 715.0; // cm from IP
74  } else if ( hit->detectorId() == kFcsHcalNorthDetId || hit->detectorId() == kFcsHcalSouthDetId ){
75  detOffset = 807.0; // cm from IP
76  } else if ( hit->detectorId() == kFcsPresNorthDetId || hit->detectorId() == kFcsPresSouthDetId ){
77  StEpdGeom epdgeo;
78  double zepd=375.0;
79  int pp,tt,n;
80  double x[5],y[5];
81  fcsDb->getEPDfromId(hit->detectorId(),hit->id(),pp,tt);
82  epdgeo.GetCorners(100*pp+tt,&n,x,y);
83  double x0 = (x[0] + x[1] + x[2] + x[3]) / 4.0;
84  double y0 = (y[0] + y[1] + y[2] + y[3]) / 4.0;
85  xyz.setX(x0);
86  xyz.setY(y0);
87  xyz.setZ(zepd);
88  }
89  mXYZ.SetXYZ( xyz.x(), xyz.y(), xyz.z() + detOffset );
90  mHit = hit;
91 }
92 
93 StFwdQAMaker::StFwdQAMaker() : StMaker("fwdQAMaker"), mTreeFile(nullptr), mTree(nullptr) {
94  setLocalOutputFile( "./fwdHists.root" ); // default off
95 }
96 
97 int StFwdQAMaker::Init() {
98 
99  mTreeFile = new TFile( mTreeFilename.Data(), "RECREATE");
100  mTree = new TTree("fwd", "fwd tracking tree");
101 
102  mTree->Branch("header", &mTreeData. header, 3200, 99 );
103  mTreeData.mcTracks.createBranch(mTree, "mcTracks");
104  mTree->Branch("nSeedTracks", &mTreeData.nSeedTracks, "nSeedTracks/I");
105  mTreeData.fstPoints.createBranch(mTree, "fstHits");
106  mTreeData.fttPoints.createBranch(mTree, "fttPoints");
107  mTreeData.fttClusters.createBranch(mTree, "fttClusters");
108  mTreeData.wcal.createBranch(mTree, "wcalClusters");
109  mTreeData.hcal.createBranch(mTree, "hcalClusters");
110 
111  mTreeData.wcalHits.createBranch(mTree, "wcalHits");
112  mTreeData.hcalHits.createBranch(mTree, "hcalHits");
113  mTreeData.epdHits.createBranch(mTree, "epdHits");
114 
115  mTreeData.reco.createBranch(mTree, "reco");
116  mTreeData.seeds.createBranch(mTree, "seeds");
117 
118 
119 
120  //========================================================================================================= adding histograms (new)
121  AddHist( mHists["fwdMultFailed"] = new TH1F("fwdMultFailed", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
122  AddHist( mHists["fwdMultAll"] = new TH1F("fwdMultAll", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
123  AddHist( mHists["fwdMultGood"] = new TH1F("fwdMultGood", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
124  AddHist( mHists["fwdMultFST"] = new TH1F("fwdMultFST", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
125  AddHist( mHists["nHitsFit"] = new TH1F("nHitsFit", ";nHitsFit; counts", 10, 0, 10) );
126  AddHist( mHists["fwdMultEcalMatch"] = new TH1F("fwdMultEcalMatch", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
127  AddHist( mHists["fwdMultHcalMatch"] = new TH1F("fwdMultHcalMatch", ";N_{ch}^{FWD}; counts", 100, 0, 100) );
128  AddHist( mHists["fwdMultEcalClusters"] = new TH1F("fwdMultEcalClusters", ";N_{Clu}^{ECAL}; counts", 100, 0, 100) );
129  AddHist( mHists["fwdMultHcalClusters"] = new TH1F("fwdMultHcalClusters", ";N_{Clu}^{HCAL}; counts", 100, 0, 100) );
130  AddHist( mHists["eta"] = new TH1F("eta", ";#eta; counts", 100, 0, 5) );
131  AddHist( mHists["phi"] = new TH1F("phi", ";#phi; counts", 100, -3.1415926, 3.1415926) );
132  AddHist( mHists["pt"] = new TH1F("pt", "; pT; counts", 500, 0, 10) );
133  AddHist( mHists["charge"] = new TH1F("charge", "; charge; counts", 4, -2, 2) );
134  AddHist( mHists["ecalMatchPerTrack"] = new TH1F("ecalMatchPerTrack", ";N_{match} / track; counts", 5, 0, 5) );
135  AddHist( mHists["hcalMatchPerTrack"] = new TH1F("hcalMatchPerTrack", ";N_{match} / track; counts", 5, 0, 5) );
136  AddHist( mHists["matchedEcalEnergy"] = new TH1F("matchedEcalEnergy", ";Energy; counts", 100, 0, 15) );
137  AddHist( mHists["matchedHcalEnergy"] = new TH1F("matchedHcalEnergy", ";Energy; counts", 100, 0, 15) );
138  AddHist( mHists["ecalEnergy"] = new TH1F("ecalEnergy", ";Energy; counts", 100, 0, 15) );
139  AddHist( mHists["hcalEnergy"] = new TH1F("hcalEnergy", ";Energy; counts", 100, 0, 15) );
140  AddHist( mHists["ecalXY"] = new TH2F( "ecalXY", ";ecalX;ecalY", 200, -200, 200, 200, -200, 200 ) );
141  AddHist( mHists["hcalXY"] = new TH2F( "hcalXY", ";hcalX;hcalY", 200, 0, 50, 200, 0, 50 ) );
142  AddHist( mHists["ecaldX"] = new TH1F( "ecaldX", ";dx (trk - ecal); counts", 400, -200, 200 ) );
143  AddHist( mHists["matchedEcaldX"] = new TH1F( "matchedEcaldX", ";dx (trk - ecal); counts", 400, -200, 200 ) );
144  AddHist( mHists["ecaldY"] = new TH1F( "ecaldY", ";dy (trk - ecal); counts", 400, -200, 200 ) );
145  AddHist( mHists["matchedEcaldY"] = new TH1F( "matchedEcaldY", ";dy (trk - ecal); counts", 400, -200, 200 ) );
146  AddHist( mHists["ecaldR"] = new TH1F( "ecaldR", ";dr (trk - ecal); counts", 400, 0, 400 ) );
147  AddHist( mHists["ecalMindR"] = new TH1F( "ecalMindR", ";dr (trk - ecal); counts", 400, 0, 400 ) );
148  AddHist( mHists["matchedEcaldR"] = new TH1F( "matchedEcaldR", ";dr (trk - ecal); counts", 400, 0, 400 ) );
149  AddHist( mHists["hcaldX"] = new TH1F( "hcaldX", ";dx (trk - hcal); counts", 400, -200, 200 ) );
150  AddHist( mHists["hcaldXdNFit"] = new TH2F( "hcaldXdNFit", ";dx (trk - hcal); nFit", 400, -200, 200, 10, 0, 10 ) );
151  AddHist( mHists["matchedHcaldX"] = new TH1F( "matchedHcaldX", ";dx (trk - hcal); counts", 400, -200, 200 ) );
152  AddHist( mHists["hcaldY"] = new TH1F( "hcaldY", ";dy (trk - hcal); counts", 400, -200, 200 ) );
153  AddHist( mHists["hcaldYdNFit"] = new TH2F( "hcaldYdNFit", ";dy (trk - hcal); nFit", 400, -200, 200, 10, 0, 10 ) );
154  AddHist( mHists["matchedHcaldY"] = new TH1F( "matchedHcaldY", ";dy (trk - hcal); counts", 400, -200, 200 ) );
155  AddHist( mHists["hcaldR"] = new TH1F( "hcaldR", ";dr (trk - hcal); counts", 400, 0, 400 ) );
156  AddHist( mHists["hcalMindR"] = new TH1F( "hcalMindR", ";dr (trk - hcal); counts", 400, 0, 400 ) );
157  AddHist( mHists["matchedHcaldR"] = new TH1F( "matchedHcaldR", ";dr (trk - hcal); counts", 400, 0, 400 ) );
158  AddHist( mHists["trkEcalX"] = new TH2F( "trkEcalX", ";trkX;ecalX", 300, -150, 150, 300, -150, 150 ) );
159  AddHist( mHists["trkEcalY"] = new TH2F( "trkEcalY", ";trkY;ecalY", 300, -150, 150, 300, -150, 150 ) );
160  AddHist( mHists["trkEcalMinX"] = new TH2F( "trkEcalMinX", ";trkX;ecalX", 300, -150, 150, 300, -150, 150 ) );
161  AddHist( mHists["trkEcalMinY"] = new TH2F( "trkEcalMinY", ";trkY;ecalY", 300, -150, 150, 300, -150, 150 ) );
162  AddHist( mHists["trkHcalX"] = new TH2F( "trkHcalX", ";trkX;hcalX", 300, -150, 150, 300, -150, 150 ) );
163  AddHist( mHists["trkHcalY"] = new TH2F( "trkHcalY", ";trkY;hcalY", 300, -150, 150, 300, -150, 150 ) );
164  AddHist( mHists["trkHcalMinX"] = new TH2F( "trkHcalMinX", ";trkX;hcalX", 300, -150, 150, 300, -150, 150 ) );
165  AddHist( mHists["trkHcalMinY"] = new TH2F( "trkHcalMinY", ";trkY;hcalY", 300, -150, 150, 300, -150, 150 ) );
166 
167 //===================================================================================================================================
168 
169  return kStOk;
170 }
171 
172 
173 
174 int StFwdQAMaker::Finish() {
175 
176  if ( mTreeFile && mTree ){
177  mTreeFile->cd();
178  mTree->Write();
179  mTreeFile->Write();
180  LOG_DEBUG << "StFwdQA File written" << endm;
181  }
182 
183  //beginning new
184  if ( mLocalOutputFile != "" ){
185  auto prevDir = gDirectory;
186 
187  // output file name
188  TFile *fOutput = new TFile(mLocalOutputFile, "RECREATE");
189  fOutput->cd();
190  for (auto nh : mHists) {
191  nh.second->SetDirectory(gDirectory);
192  nh.second->Write();
193  }
194 
195  // restore previous directory
196  gDirectory = prevDir;
197 
198  LOG_INFO << "Done writing StFwdQAMaker output to local file : " << mLocalOutputFile << endm;
199  }//end new
200 
201  return kStOk;
202 }
203 
204 
205 int StFwdQAMaker::Make() {
206  LOG_INFO << "FWD Report:" << endm;
207  StEvent *mStEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
208  if ( mStEvent ){
209  // report number of fwd tracks
210  auto fwdTracks = mStEvent->fwdTrackCollection();
211  LOG_INFO << "Number of FwdTracks (StFwdTrackCollection): " << fwdTracks->tracks().size() << endm;
212  if ( mStEvent->fttCollection() ){
213  LOG_INFO << "Number of Ftt Points (StEvent)" << mStEvent->fttCollection()->points().size() << endm;
214  }
215  }
216  LOG_INFO << "SETUP START" << endm;
217  // setup the datasets / makers
218 
219 
220  mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
221  if(mMuDstMaker) {
222  mMuDst = mMuDstMaker->muDst();
223  mMuForwardTrackCollection = mMuDst->muFwdTrackCollection();
224  mMuFcsCollection = mMuDst->muFcsCollection();
225  if (mMuDst->event())
226  mTreeData.header.run = mMuDst->event()->runNumber();
227  if (mMuForwardTrackCollection){
228  LOG_DEBUG << "Number of StMuFwdTracks: " << mMuForwardTrackCollection->numberOfFwdTracks() << endm;
229  }
230  else{
231  LOG_DEBUG << "No muFwdTrackCollection " << endm;
232  }
233  } else {
234  LOG_DEBUG << "No StMuDstMaker found: " << mMuDstMaker << endm;
235  }
236 
237  mFcsDb = static_cast<StFcsDb *>(GetDataSet("fcsDb"));
238  mFwdTrackMaker = (StFwdTrackMaker*) GetMaker( "fwdTrack" );
239  if (!mFwdTrackMaker) {
240  LOG_WARN << "No StFwdTrackMaker found, skipping StFwdQAMaker" << endm;
241  // return kStOk;
242  }
243 
244 
245  LOG_DEBUG << "SETUP COMPLETE" << endm;
246  ProcessFwdTracks();
247  ProcessFwdMuTracks();
248 
249  FillMcTracks();
250  FillTracks(); //maybe not running
251  FillFstPoints(); //no fst
252  FillFttClusters();
253  FillFcsStMuDst();
254  mTree->Fill();
255 
256  return kStOk;
257 }
258 void StFwdQAMaker::Clear(const Option_t *opts) {
259  mTreeData.clear();
260  return;
261 }
262 
263 void StFwdQAMaker::FillFstPoints(){
264  StMuFstCollection * fst = mMuDst->muFstCollection();
265  if (!fst) {
266  LOG_WARN << "No StMuFstCollection ... bye-bye" << endm;
267  return;
268  }
269 
270  LOG_INFO << "Loading " << fst->numberOfHits() << " StMuFstHits" << endm;
271  for ( unsigned int index = 0; index < fst->numberOfHits(); index++){
272  StMuFstHit * muFstHit = fst->getHit( index );
273  mTreeData.fstPoints.add( muFstHit );
274  } // index
275 }
276 
277 void StFwdQAMaker::FillTracks() {
278  mTreeData.nSeedTracks = 0;
279  if ( mMuForwardTrackCollection ){
280  LOG_DEBUG << "Adding " << mMuForwardTrackCollection->numberOfFwdTracks() << " FwdTracks (MuDst)" << endm;
281  for ( size_t iTrack = 0; iTrack < mMuForwardTrackCollection->numberOfFwdTracks(); iTrack++ ){
282  auto muTrack = mMuForwardTrackCollection->getFwdTrack(iTrack);
283  mTreeData.reco.add( muTrack );
284 
285  for (auto fsth : muTrack->mFSTPoints){
286  mTreeData.seeds.add( fsth );
287  mTreeData.nSeedTracks++;
288  }
289  for (auto ftth : muTrack->mFTTPoints){
290  mTreeData.seeds.add( ftth );
291  mTreeData.nSeedTracks++;
292  }
293  if ( iTrack > 5000 ) {
294  LOG_WARN << "Truncating to 5000 tracks" << endm;
295  break;
296  }
297  }
298  } else {
299  LOG_WARN << "No StMuFwdTrackCollection found" << endm;
300  }
301  LOG_DEBUG << "TRACKS COMPLETE" << endm;
302 }
303 
304 void StFwdQAMaker::FillFcsStMuDst( ) {
305 
306  if ( !mMuDst ){
307  LOG_DEBUG << "No mMuDst found, skipping StFwdQAMaker::FillFcsStEvent" << endm;
308  return;
309  }
310  StMuFcsCollection* fcs = mMuDst->muFcsCollection();
311  if ( !fcs ){
312  LOG_DEBUG << "No muFcsCollection found, skipping StFwdQAMaker::FillFcsStEvent" << endm;
313  return;
314  }
315 
316  StFcsDb* fcsDb=static_cast<StFcsDb*>(GetDataSet("fcsDb"));
317 
318  // LOAD ECAL / HCAL CLUSTERS
319  LOG_INFO << "MuDst has #fcs clusters: " << fcs->numberOfClusters() << endm;
320  for( size_t i = 0; i < fcs->numberOfClusters(); i++){
321  StMuFcsCluster * clu = fcs->getCluster(i);
322  FcsClusterWithStarXYZ *cluSTAR = new FcsClusterWithStarXYZ(clu, fcsDb);
323  if ( clu->detectorId() == kFcsEcalNorthDetId || clu->detectorId() == kFcsEcalSouthDetId ){
324  LOG_INFO << "Adding WCAL Cluster to FwdTree" << endm;
325  mTreeData.wcal.add( cluSTAR );
326  } else if ( clu->detectorId() == kFcsHcalNorthDetId || clu->detectorId() == kFcsHcalSouthDetId ){
327  LOG_INFO << "Adding HCAL Cluster to FwdTree" << endm;
328  mTreeData.hcal.add( cluSTAR );
329  }
330 
331  delete cluSTAR;
332  }
333 
334  // LOAD ECAL / HCAL CLUSTERS
335  LOG_INFO << "MuDst has #fcs hits: " << fcs->numberOfHits() << endm;
336  for( size_t i = 0; i < fcs->numberOfHits(); i++){
337  StMuFcsHit * hit = fcs->getHit(i);
338  FcsHitWithStarXYZ *hitSTAR = new FcsHitWithStarXYZ(hit, fcsDb);
339  if ( hit->detectorId() == kFcsEcalNorthDetId || hit->detectorId() == kFcsEcalSouthDetId ){
340  LOG_DEBUG << "Adding WCAL Cluster to FwdTree" << endm;
341  mTreeData.wcalHits.add( hitSTAR );
342  } else if ( hit->detectorId() == kFcsHcalNorthDetId || hit->detectorId() == kFcsHcalSouthDetId ){
343  LOG_DEBUG << "Adding HCAL Cluster to FwdTree" << endm;
344  mTreeData.hcalHits.add( hitSTAR );
345  } else if ( hit->detectorId() == kFcsPresNorthDetId || hit->detectorId() == kFcsPresSouthDetId ){
346  LOG_DEBUG << "Adding PRES hit to FwdTree" << endm;
347  mTreeData.epdHits.add( hitSTAR );
348  }
349  delete hitSTAR;
350  }
351 
352  // TODO, cleanup?
353 }
354 
355 void StFwdQAMaker::FillMcTracks(){
356  // Retrieve pointer to MC tracks
357  TClonesArray *mcTracks = mMuDst->mcArray(1);
358  LOG_INFO << "MuDst has #mc tracks: " << mcTracks->GetEntriesFast() << endm;
359  // Loop over MC vertices
360  for (Int_t iVtx=0; iVtx<mcTracks->GetEntriesFast(); iVtx++) {
361  // Retrieve i-th MC vertex from MuDst
362  StMuMcTrack *mcTrack = (StMuMcTrack*)mcTracks->UncheckedAt(iVtx);
363  if ( !mcTrack ) continue;
364 
365  // Add MC track to the tree
366  mTreeData.mcTracks.add( mcTrack );
367  }
368 }
369 
370 
371 void StFwdQAMaker::FillFttClusters(){
372 
373  auto muFttCollection = mMuDst->muFttCollection();
374  if ( muFttCollection ){
375  LOG_DEBUG << "MuDst has #ftt clusters: " << muFttCollection->numberOfClusters() << endm;
376  for ( size_t i = 0; i < muFttCollection->numberOfClusters(); i++ ){
377  StMuFttCluster * c = muFttCollection->getCluster(i);
378  mTreeData.fttClusters.add( c );
379  }
380 
381  for ( size_t i = 0; i < muFttCollection->numberOfPoints(); i++ ){
382  StMuFttPoint * c = muFttCollection->getPoint(i);
383  mTreeData.fttPoints.add( c );
384  }
385  }
386  else{
387  LOG_INFO << "no muFttCollection " << endm;
388  }
389 }
390 
391 //==================================================adding new function
392 
393 void StFwdQAMaker::ProcessFwdTracks( ){
394  // This is an example of how to process fwd track collection
395  LOG_DEBUG << "StFwdAnalysisMaker::ProcessFwdTracks" << endm;
396  StEvent *stEvent = static_cast<StEvent *>(GetInputDS("StEvent"));
397  if (!stEvent)
398  return;
399 
400  if (stEvent){
401  StFttCollection *fttCol = stEvent->fttCollection();
402  if (fttCol){
403  LOG_DEBUG << "The Ftt Collection has " << fttCol->numberOfPoints() << " points" << endm;
404  }
405  }
406  StFwdTrackCollection * ftc = stEvent->fwdTrackCollection();
407  if (!ftc) {
408  LOG_DEBUG << "Forward Track Collection is not present" << endm;
409  return;
410  }
411 
412  LOG_DEBUG << "Checking FcsCollection" << endm;
413  StFcsCollection *fcs = stEvent->fcsCollection();
414  if (!fcs) return;
415 
416  StFcsDb *mFcsDb = static_cast<StFcsDb *>(GetDataSet("fcsDb"));
417 
418  size_t fwdMultEcalMatch = 0;
419  size_t fwdMultHcalMatch = 0;
420  size_t fwdMultFST = 0;
421 
422  LOG_INFO << "FwdTrackCollection has: " << ftc->tracks().size() << " tracks" << endm;
423 
424  getHist( "fwdMultAll" )->Fill( ftc->tracks().size() );
425 
426  // Cluster info (independen t of tracks)
427  size_t fwdMultEcalClusters = 0;
428  size_t fwdMultHcalClusters = 0;
429  for ( int iDet = 0; iDet < 4; iDet++ ){
430  for( size_t i = 0; i < fcs->clusters(iDet).size(); i++){
431  StFcsCluster * clu = fcs->clusters(iDet)[i];
432 
433  if ( iDet < 2 ){
434  fwdMultEcalClusters++;
435  getHist( "ecalEnergy" )->Fill( clu->energy() );
436  } else if ( iDet < 4 ){
437  fwdMultHcalClusters++;
438  getHist( "hcalEnergy" )->Fill( clu->energy() );
439  }
440  }
441  }
442 
443  getHist( "fwdMultEcalClusters" )->Fill( fwdMultEcalClusters );
444  getHist( "fwdMultHcalClusters" )->Fill( fwdMultHcalClusters );
445 
446 
447  size_t nGood = 0;
448  size_t nFailed = 0;
449  for ( auto fwdTrack : ftc->tracks() ){
450  if ( !fwdTrack->didFitConvergeFully() ) {
451  nFailed++;
452  continue;
453  }
454  nGood++;
455  LOG_DEBUG << TString::Format("StFwdTrack[ nProjections=%lu, nFTTSeeds=%lu, nFSTSeeds=%lu, mPt=%f ]", fwdTrack->mProjections.size(), fwdTrack->mFTTPoints.size(), fwdTrack->mFSTPoints.size(), fwdTrack->momentum().perp()) << endm;
456  LOG_DEBUG << "track fit momentum " << TString::Format( "(pt=%f, eta=%f, phi=%f)", fwdTrack->momentum().perp(), fwdTrack->momentum().pseudoRapidity(), fwdTrack->momentum().phi() ) << endm;
457  LOG_DEBUG << "StFwdTrack has " << fwdTrack->ecalClusters().size() << " ecal matches" << endm;
458  LOG_DEBUG << "StFwdTrack has " << fwdTrack->hcalClusters().size() << " hcal matches" << endm;
459 
460  getHist("ecalMatchPerTrack")->Fill( fwdTrack->ecalClusters().size() );
461  getHist("hcalMatchPerTrack")->Fill( fwdTrack->hcalClusters().size() );
462 
463  getHist( "nHitsFit" )->Fill( fwdTrack->numberOfFitPoints() );
464 
465  if (fwdTrack->mFSTPoints.size() > 0){
466  fwdMultFST ++;
467  }
468 
469  getHist("eta")->Fill( fwdTrack->momentum().pseudoRapidity() );
470  getHist("phi")->Fill( fwdTrack->momentum().phi() );
471  getHist("pt")->Fill( fwdTrack->momentum().perp() );
472 
473  getHist("charge")->Fill( fwdTrack->charge() );
474 
475  // ecal proj
476  int detId = kFcsWcalId;
477  TVector3 ecalXYZ;
478  TVector3 ecapP;
479 
480  StFwdTrackProjection ecalProj = fwdTrack->getProjectionFor( detId, 0 );
481  StFwdTrackProjection hcalProj = fwdTrack->getProjectionFor( kFcsHcalId, 0 );
482  LOG_DEBUG << "EcalProj z= " << ecalProj.mXYZ.z() << endm;
483  LOG_DEBUG << "HcalProj z= " << hcalProj.mXYZ.z() << endm;
484  LOG_DEBUG << "EcalProj Mom" << TString::Format( "(pt=%f, eta=%f, phi=%f)", ecalProj.mMom.perp(), ecalProj.mMom.pseudoRapidity(), ecalProj.mMom.phi() ) << endm;
485 
486  for ( size_t iEcal = 0; iEcal < fwdTrack->ecalClusters().size(); iEcal++ ){
487  StFcsCluster *clu = fwdTrack->ecalClusters()[iEcal];
488  LOG_DEBUG << "Ecal clu detId = " << clu->detectorId() << endm;
489  getHist("matchedEcalEnergy")->Fill( clu->energy() );
490 
491  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(clu->detectorId(), clu->x(), clu->y());
492  float dx = ecalProj.mXYZ.x() - xyz.x();
493  float dy = ecalProj.mXYZ.y() - xyz.y();
494  float dr = sqrt(dx*dx + dy*dy);
495  getHist("matchedEcaldX")->Fill( dx );
496  getHist("matchedEcaldY")->Fill( dy );
497  getHist("matchedEcaldR")->Fill( dr );
498  }
499 
500  if (ecalProj.mXYZ.z() > 500){
501  double mindR = 999;
502  StFcsCluster * cclu = nullptr; // closet cluster
503  for ( int iDet = 0; iDet < 2; iDet++ ){
504  for( size_t i = 0; i < fcs->clusters(iDet).size(); i++){
505  StFcsCluster * clu = fcs->clusters(iDet)[i];
506 
507  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(clu->detectorId(), clu->x(), clu->y());
508  getHist("ecalXY")->Fill( xyz.x(), xyz.y() );
509 
510  float dx = ecalProj.mXYZ.x() - xyz.x();
511  float dy = ecalProj.mXYZ.y() - xyz.y();
512  float dr = sqrt(dx*dx + dy*dy);
513 
514  if ( fabs(dy) < 25 )
515  getHist( "ecaldX" )->Fill( dx );
516  if ( fabs(dx) < 25 )
517  getHist( "ecaldY" )->Fill( dy );
518  getHist( "ecaldR" )->Fill( dr );
519  if ( dr < mindR ){
520  mindR = dr;
521  cclu = clu;
522  }
523 
524  getHist( "trkEcalX" ) -> Fill( ecalProj.mXYZ.x(), xyz.x() );
525  getHist( "trkEcalY" ) -> Fill( ecalProj.mXYZ.y(), xyz.y() );
526 
527  }
528  }
529  getHist( "ecalMindR" )->Fill( mindR );
530  if (cclu){
531  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(cclu->detectorId(), cclu->x(), cclu->y());
532  getHist( "trkEcalMinX" ) -> Fill( ecalProj.mXYZ.x(), xyz.x() );
533  getHist( "trkEcalMinY" ) -> Fill( ecalProj.mXYZ.y(), xyz.y() );
534  }
535  }
536 
537  if (hcalProj.mXYZ.z() > 500){
538 
539  double mindR = 999;
540  StFcsCluster * cclu = nullptr;
541  for ( int iDet = 2; iDet < 4; iDet++ ){
542  for( size_t i = 0; i < fcs->clusters(iDet).size(); i++){
543  StFcsCluster * clu = fcs->clusters(iDet)[i];
544  if (!clu) continue;
545  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(clu->detectorId(), clu->x(), clu->y());
546  getHist("hcalXY")->Fill( xyz.x(), xyz.y() );
547 
548  float dx = hcalProj.mXYZ.x() - xyz.x();
549  float dy = hcalProj.mXYZ.y() - xyz.y();
550  float dr = sqrt(dx*dx + dy*dy);
551 
552  if ( fabs(dy) < 25 ){
553  getHist( "hcaldX" )->Fill( dx );
554  getHist( "hcaldXdNFit" )->Fill( dx, fwdTrack->numberOfFitPoints() );
555 
556  }
557  if ( fabs(dx) < 25 ){
558  getHist( "hcaldY" )->Fill( dy );
559  getHist( "hcaldYdNFit" )->Fill( dy, fwdTrack->numberOfFitPoints() );
560  }
561  getHist( "hcaldR" )->Fill( dr );
562 
563  if ( dr < mindR ){
564  mindR = dr;
565  cclu = clu;
566  }
567 
568  getHist( "trkHcalX" ) -> Fill( hcalProj.mXYZ.x(), xyz.x() );
569  getHist( "trkHcalY" ) -> Fill( hcalProj.mXYZ.y(), xyz.y() );
570  }
571  }
572  getHist( "hcalMindR" )->Fill( mindR );
573  if (cclu){
574  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(cclu->detectorId(), cclu->x(), cclu->y());
575  getHist( "trkHcalMinX" ) -> Fill( hcalProj.mXYZ.x(), xyz.x() );
576  getHist( "trkHcalMinY" ) -> Fill( hcalProj.mXYZ.y(), xyz.y() );
577  }
578  }
579 
580  if (fwdTrack->ecalClusters().size() > 0)
581  fwdMultEcalMatch++;
582  if (fwdTrack->hcalClusters().size() > 0)
583  fwdMultHcalMatch++;
584 
585  } // Loop ftc->tracks()
586 
587  getHist( "fwdMultGood" )->Fill( nGood );
588  getHist( "fwdMultFailed" )->Fill( nFailed );
589  getHist("fwdMultFST")->Fill( fwdMultFST );
590  getHist("fwdMultHcalMatch")->Fill( fwdMultHcalMatch );
591  getHist("fwdMultEcalMatch")->Fill( fwdMultEcalMatch );
592 
593  LOG_INFO << "Found " << nFailed << " failed track fits out of " << ftc->tracks().size() << endm;
594 } // ProcessFwdTracks
595 //========================================================end of added function
596 
597 
598 
599 void StFwdQAMaker::ProcessFwdMuTracks( ){
600  // This is an example of how to process fwd track collection
601  LOG_DEBUG << "StFwdAnalysisMaker::ProcessFwdMuTracks" << endm;
602  StMuDstMaker *mMuDstMaker = (StMuDstMaker *)GetMaker("MuDst");
603  if(!mMuDstMaker) {
604  LOG_WARN << " No MuDstMaker ... bye-bye" << endm;
605  return;
606  }
607  StMuDst *mMuDst = mMuDstMaker->muDst();
608  if(!mMuDst) {
609  LOG_WARN << " No MuDst ... bye-bye" << endm;
610  return;
611  }
612  StMuFwdTrackCollection * ftc = mMuDst->muFwdTrackCollection();
613  if (!ftc) return;
614 
615  StMuFcsCollection *fcs = mMuDst->muFcsCollection();
616  if (!fcs) return;
617 
618  LOG_INFO << "Number of StMuFwdTracks: " << ftc->numberOfFwdTracks() << endl;
619 
620  StFcsDb *mFcsDb = static_cast<StFcsDb *>(GetDataSet("fcsDb"));
621 
622  size_t fwdMultFST = 0;
623  size_t fwdMultEcalMatch = 0;
624  size_t fwdMultHcalMatch = 0;
625 
626  for ( size_t iTrack = 0; iTrack < ftc->numberOfFwdTracks(); iTrack++ ){
627  StMuFwdTrack * muFwdTrack = ftc->getFwdTrack( iTrack );
628  // LOG_DEBUG << TString::Format("StMuFwdTrack[ nProjections=%lu, nFTTSeeds=%lu, nFSTSeeds=%lu, mPt=%f ]", muFwdTrack->mProjections.size(), muFwdTrack->mFTTPoints.size(), muFwdTrack->mFSTPoints.size(), muFwdTrack->momentum().Pt()) << endm;
629 
630  LOG_DEBUG << "StMuFwdTrack has " << muFwdTrack->mEcalClusters.GetEntries() << " Ecal matched" << endm;
631  LOG_DEBUG << "StMuFwdTrack has " << muFwdTrack->mHcalClusters.GetEntries() << " Hcal matched" << endm;
632 
633  getHist("eta")->Fill( muFwdTrack->momentum().Eta() );
634  getHist("phi")->Fill( muFwdTrack->momentum().Phi() );
635 
636  if (muFwdTrack->mFSTPoints.size() > 0){
637  fwdMultFST ++;
638  }
639 
640  if (muFwdTrack->mEcalClusters.GetEntries() > 0)
641  fwdMultEcalMatch++;
642  if (muFwdTrack->mHcalClusters.GetEntries() > 0)
643  fwdMultHcalMatch++;
644 
645 
646  // ecal proj
647  int detId = kFcsWcalId;
648  TVector3 ecalXYZ;
649  TVector3 ecapP;
650 
651  StMuFwdTrackProjection ecalProj;
652  bool foundEcalProj = muFwdTrack->getProjectionFor( detId, ecalProj, 0 );
653 
654  if (foundEcalProj){
655  for( size_t i = 0; i < fcs->numberOfClusters(); i++){
656  StMuFcsCluster * clu = fcs->getCluster(i);
657 
658  if ( clu->detectorId() > 1 ) continue;
659 
660  if ( clu->energy() < 1 ) continue;
661  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(clu->detectorId(), clu->x(), clu->y());
662 
663  float dx = ecalProj.mXYZ.X() - xyz.x();
664  float dy = ecalProj.mXYZ.Y() - xyz.y();
665  float dr = sqrt(dx*dx + dy*dy);
666 
667  getHist( "ecaldX" )->Fill( dx );
668  getHist( "ecaldY" )->Fill( dy );
669  getHist( "ecaldR" )->Fill( dr );
670 
671  getHist( "trkEcalX" ) -> Fill( ecalProj.mXYZ.X(), xyz.x() );
672 
673  } // i
674  } // foundEcalProj
675 
676 
677  for ( int i = 0; i < muFwdTrack->mEcalClusters.GetEntries(); i++ ){
678  auto c = (StMuFcsCluster*) muFwdTrack->mEcalClusters.At(i);
679  if (!c) continue;
680  getHist("ecalEnergy")->Fill( c->energy() );
681 
682  LOG_DEBUG << "eCal Cluster detId = " << c->detectorId() << endm;
683  StThreeVectorD xyz = mFcsDb->getStarXYZfromColumnRow(c->detectorId(), c->x(), c->y());
684  getHist("ecalXY")->Fill( xyz.x(), xyz.y() );
685 
686  if (foundEcalProj){
687  getHist("matchedEcaldX")->Fill( ecalProj.mXYZ.X() - xyz.x() );
688  }
689  } // i
690 
691  getHist("ecalMatchPerTrack")->Fill( muFwdTrack->mEcalClusters.GetEntries() );
692  getHist("hcalMatchPerTrack")->Fill( muFwdTrack->mHcalClusters.GetEntries() );
693 
694  for ( int i = 0; i < muFwdTrack->mHcalClusters.GetEntries(); i++ ){
695  auto c = (StMuFcsCluster*) muFwdTrack->mHcalClusters.At(i);
696  if (!c) continue;
697  getHist("hcalEnergy")->Fill( c->energy() );
698 
699  getHist("hcalXY")->Fill( c->x(), c->y() );
700  } // i
701  } // iTrack
702 
703  getHist("fwdMult")->Fill( ftc->numberOfFwdTracks() );
704  getHist("fwdMultFST")->Fill( fwdMultFST );
705  getHist("fwdMultHcalMatch")->Fill( fwdMultHcalMatch );
706  getHist("fwdMultEcalMatch")->Fill( fwdMultEcalMatch );
707 }
StMuDst::muFcsCollection
static StMuFcsCollection * muFcsCollection()
returns pointer to current StMuFcsCollection
Definition: StMuDst.h:395
StMuDst::muFttCollection
static StMuFttCollection * muFttCollection()
returns pointer to current StMuFttCollection
Definition: StMuDst.h:397
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
FwdQATreeData::header
FwdTreeHeader header
Primary event vertex.
Definition: StFwdQAMaker.h:154
FwdQATreeData::mcTracks
TClonesArrayWriter< StMuMcTrack > mcTracks
MC tracks.
Definition: StFwdQAMaker.h:156
StFcsDb::getStarXYZfromColumnRow
StThreeVectorD getStarXYZfromColumnRow(int det, float col, float row, float FcsZ=-1.0) const
get the STAR frame cooridnates from other way
Definition: StFcsDb.cxx:863
StMuDst::muFwdTrackCollection
static StMuFwdTrackCollection * muFwdTrackCollection()
returns pointer to current StMuFwdTrackCollection
Definition: StMuDst.h:401
StFcsDb::getStarXYZ
StThreeVectorD getStarXYZ(int det, float FcsX, float FcsY, float FcsZ=-1.0, float zVertex=0.0) const
get the STAR frame cooridnates from local XYZ [cm]
Definition: StFcsDb.cxx:868
StFcsDb::getEPDfromId
void getEPDfromId(int det, int id, int &pp, int &tt)
Get FCS&#39;s EPD map foom EPD mapping.
Definition: StFcsDb.cxx:1484
StMuDst::muFstCollection
static StMuFstCollection * muFstCollection()
returns pointer to current StMuFstCollection
Definition: StMuDst.h:399
StFwdQAMaker::getHist
TH1 * getHist(TString n)
Get the Hist object from the map.
Definition: StFwdQAMaker.h:233
FcsClusterWithStarXYZ
Store Cluster with STAR XYZ position.
Definition: StFwdQAMaker.h:119
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
FcsHitWithStarXYZ
Store Hit with STAR XYZ position.
Definition: StFwdQAMaker.h:135
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5