StRoot: StMcAnalysisMaker/StMcAnalysisMaker.cxx Source File

00001 /*************************************************
00002  *
00003  * $Id: StMcAnalysisMaker.cxx,v 1.38 2011/04/03 16:02:40 fisyak Exp $
00004  * $Log: StMcAnalysisMaker.cxx,v $
00005  * Revision 1.38  2011/04/03 16:02:40  fisyak
00006  * Fix effect of constness in StAssociationMaker
00007  *
00008  * Revision 1.37  2010/06/22 16:23:25  fine
00009  * introduce the correct const StMcTarck * pointer cast
00010  *
00011  * Revision 1.36  2010/05/10 17:15:34  fine
00012  * RT # 1932 Remove the redundant correction. Restore 1.34 rev
00013  *
00014  * Revision 1.35  2010/05/07 20:17:18  fine
00015  * Add CPP macro to separate McTracks
00016  *
00017  * Revision 1.34  2010/01/28 18:11:46  perev
00018  * WarnOff
00019  *
00020  * Revision 1.33  2007/04/17 05:08:36  perev
00021  * GetTFile()==>StMaker. Jerome request
00022  *
00023  * Revision 1.32  2007/03/21 17:14:36  fisyak
00024  * Add keys for switch off Tpc, Svt and Ssd NTuples
00025  *
00026  * Revision 1.31  2006/06/20 02:39:38  calderon
00027  * Fixed several changes from the last version which made this package not
00028  * executable using the default StAssociator.C macro.
00029  * Removed dependence on StBFChain, so the default StAssociator.C macro can work.
00030  * Removed #if 0  commands which had taken out the track ntuple.
00031  * Added back the missing code to produce the first 8 entries of the ntuple.
00032  * Reverted to the file handling as was done before so that StAssociator.C can work.
00033  * Kept the ntuples added by Yuri (but did not check them).
00034  *
00035  * Revision 1.30  2005/11/22 21:51:53  fisyak
00036  * Add NTuple for svt and ssd hit
00037  *
00038  * Revision 1.29  2005/09/28 22:52:52  calderon
00039  * Changed access to StMcEvent to use GetDataSet to be consistent with persistent StMcEvent.
00040  *
00041  * Revision 1.28  2005/07/19 22:04:49  perev
00042  * MultiVertex
00043  *
00044  * Revision 1.27  2004/03/30 03:11:18  calderon
00045  * Added information about the matching into the track Ntuple:
00046  *  - Dominatrack... (track with most common IdTruth, a.k.a. dominant contributor
00047  *  - nHitsIdTruth (should be very similar to commonTpcHits, but the distance cut
00048  *    option might make them different)
00049  *  - n MC hits for the dominatrack
00050  *  - n Fit points for the reconstructed track
00051  *  - n Points (from StDetectorInfo).
00052  *  - "quality" = average hit quality of the reco hits belonging to the dominatrack.
00053  *
00054  * Revision 1.26  2004/02/08 00:13:05  calderon
00055  * Check that the size of the map is non-zero to get the first track.
00056  *
00057  * Revision 1.25  2004/02/08 00:10:58  calderon
00058  * Histogram of rc hits is now not restricted to TPC only.   Histogram of mc hits is
00059  * done looping over tpc and svt hits.
00060  * If the first entry of the track nodes is not associated, use the first entry of
00061  * the map.
00062  *
00063  * Revision 1.24  2004/01/24 04:35:53  calderon
00064  * The last commits using BFChain broke StAssociator.  The histograms were
00065  * not found and nothing was plotted.  Revert back to the usual mode (this
00066  * code is not meant to run in bfc anyway, it's meant as a collection of
00067  * examples).
00068  *
00069  * Revision 1.23  2004/01/17 19:12:28  fisyak
00070  * Add protection when StMcAnalysisMaker is not running within StBFChain and/or TFile is not coming from StBFChain parameters
00071  *
00072  * Revision 1.22  2004/01/13 21:06:04  fisyak
00073  * Add TpcHitNtuple usind IdTruth info
00074  *
00075  * Revision 1.21  2003/09/02 17:58:41  perev
00076  * gcc 3.2 updates + WarnOff
00077  *
00078  * Revision 1.20  2000/05/11 21:57:34  calderon
00079  * Book the histograms before creating the file so that Kathy's macros
00080  * can pick them up.
00081  *
00082  * Revision 1.19  2000/05/11 16:21:26  calderon
00083  * Write only one V0 matched pair to supress the screen output.
00084  *
00085  * Revision 1.18  2000/04/20 21:31:52  calderon
00086  * More checks for the cases where a track has no partner, Thanks Janet.
00087  *
00088  * Revision 1.17  2000/04/20 16:59:47  calderon
00089  * Pick up the makers with the new names
00090  * Change the name from "McAnalysis" to "StMcAnalysisMaker"
00091  * No longer use the helix, use the primary track momentum when available
00092  * (also avoids dependency on StEventSummary)
00093  * Denominator for mom. resolution histograms is now simulated momentum.
00094  *
00095  * Revision 1.16  2000/04/12 21:33:57  calderon
00096  * check case where no multimaps are found
00097  *
00098  * Revision 1.15  2000/04/04 23:18:12  calderon
00099  * B -> B * kilogauss
00100  *
00101  * Revision 1.14  2000/03/28 02:28:04  calderon
00102  * Return to calculating momentum at the origin, instead of taking the
00103  * momentum from the first point.
00104  *
00105  * Revision 1.13  2000/03/06 21:47:56  calderon
00106  * Add Lee's V0 example
00107  *
00108  * Revision 1.12  2000/02/07 16:43:37  calderon
00109  * Find the first hit wherever it may be, instead of assuming a sector and padrow.
00110  *
00111  * Revision 1.11  2000/01/24 22:22:24  calderon
00112  * use delete [] for the array of floats
00113  *
00114  * Revision 1.10  1999/12/14 07:08:48  calderon
00115  * First version to work with new StEvent, StMcEvent & StAssociationMaker.
00116  * Need to add more examples of all the new maps.
00117  *
00118  * Revision 1.9  1999/10/01 14:11:15  calderon
00119  * Chech to see whether StEvent has primary vertex
00120  * before trying to use it.  If no primary vertex is
00121  * found, assume its position is (0,0,0).
00122  *
00123  * Revision 1.8  1999/09/28 15:03:29  didenko
00124  * Cleanup dependencies on non existing h-files
00125  *
00126  * Revision 1.7  1999/09/23 21:25:48  calderon
00127  * Added Log & Id
00128  * Modified includes according to Yuri
00129  *
00130  * Revision 1.6  1999/09/10 19:11:15  calderon
00131  * Write the Ntuple in StMcAnalysisMaker into a file.
00132  * This way it can be accessed after the macro finishes,
00133  * otherwise it gets deleted.
00134  *
00135  * Revision 1.5  1999/09/09 23:59:37  calderon
00136  * Made the following changes:
00137  *
00138  * -book histograms and ntuple in Init()
00139  * -do not delete the histograms, they are supposed to be
00140  *  deleted automatically by the chain
00141  * -don't create the canvas here, now done in the macro
00142  *
00143  * Revision 1.4  1999/07/30 16:19:19  calderon
00144  * Use value_type typedef for inserting pairs in multimaps, Victor corrected iterators on HP in SL99h, Improved use of const for HP compilation
00145  *
00146  * Revision 1.3  1999/07/29 15:08:33  calderon
00147  * Include Mom. Resolution example (Histograms & Ntuple)
00148  *
00149  * Revision 1.2  1999/07/28 20:27:30  calderon
00150  * Version with SL99f libraries
00151  *
00152  *
00153  * Examples that use the structures of
00154  * StMcEvent and StAssociationMaker
00155  *
00156  *************************************************/
00157 #include <assert.h>
00158 #include <Stiostream.h>
00159 #include <stdlib.h>
00160 #include <string>
00161 #include <vector>
00162 #include <set>
00163 #include <map>
00164 #include <cmath>
00165 
00166 #include "TStyle.h"
00167 #include "TCanvas.h"
00168 #include "TH1.h"
00169 #include "TH2.h"
00170 #include "TNtuple.h"
00171 #include "TFile.h"
00172 
00173 #include "StMcAnalysisMaker.h"
00174 #include "PhysicalConstants.h"
00175 #include "SystemOfUnits.h"
00176 
00177 #include "StMessMgr.h"
00178 
00179 #include "StAssociationMaker/StAssociationMaker.h"
00180 #include "StAssociationMaker/StTrackPairInfo.hh"
00181 
00182 #include "StThreeVectorF.hh"
00183 
00184 #include "StEventTypes.h"
00185 
00186 #include "StMcEventTypes.hh"
00187 
00188 #include "StMcEvent.hh"
00189 
00190 // Define data Members for the histograms
00191 const Int_t   StMcAnalysisMaker::mNumDeltaX = 50;
00192 const Int_t   StMcAnalysisMaker::mNumDeltaZ = 50;
00193 const Float_t StMcAnalysisMaker::mMinDeltaX = -0.52;
00194 const Float_t StMcAnalysisMaker::mMaxDeltaX =  0.52;
00195 const Float_t StMcAnalysisMaker::mMinDeltaZ = -0.52;
00196 const Float_t StMcAnalysisMaker::mMaxDeltaZ =  0.52;
00197 struct TpcHitMRPair_t {
00198   Float_t sector, row, isDet,
00199     xM, yM, zM, pxM, pyM, pzM, dEM, dSM, nM,
00200     xR, yR, zR, dER, IdM, IdR, qR, nR;
00201 };
00202 static const Char_t *vTpcHitMRPair = "sector:row:isDet:xM:yM:zM:pxM:pyM:pzM:dEM:dSM:nM:xR:yR:zR:dER:IdM:IdR:qR:nR";
00203 static TpcHitMRPair_t TpcHitMRPair;
00204 struct SvtHitMRPair_t {
00205   Float_t barrel, layer, ladder, wafer,  hybrid, index,
00206     xM, yM, zM, pxM, pyM, pzM, dEM, dSM, nM,
00207     xR, yR, zR, dER, IdM, IdR, qR, nR;
00208 };
00209 static const Char_t *vSvtHitMRPair = "barrel:layer:ladder:wafer:hybrid:index:xM:yM:zM:pxM:pyM:pzM:dEM:dSM:nM:xR:yR:zR:dER:IdM:IdR:qR:nR";
00210 static SvtHitMRPair_t SvtHitMRPair;
00211 struct SsdHitMRPair_t {
00212   Float_t ladder, wafer,
00213     xM, yM, zM, pxM, pyM, pzM, dEM, dSM, nM,
00214     xR, yR, zR, dER, IdM, IdR, qR, nR;
00215 };
00216 static const Char_t *vSsdHitMRPair = "ladder:wafer:xM:yM:zM:pxM:pyM:pzM:dEM:dSM:nM:xR:yR:zR:dER:IdM:IdR:qR:nR";
00217 static SsdHitMRPair_t SsdHitMRPair;
00218 ClassImp(StMcAnalysisMaker)
00219 
00220 //_________________________________________________
00221 StMcAnalysisMaker::StMcAnalysisMaker(const char *name, const char *title):
00222   StMaker(name,title), mNtupleFile(0)
00223 {
00224     //  StMcAnalysisMaker Constructor
00225     // - zero all pointers defined in the header file
00226     mAssociationCanvas = 0;
00227     mMomResolution  = 0;
00228     mHitResolution  = 0;   
00229     mSvtHitResolution  = 0;   
00230     mSsdHitResolution  = 0;   
00231     coordRec        = 0;  
00232     coordMcPartner  = 0;
00233     mTrackNtuple    = 0;
00234     mTpcHitNtuple   = 0;
00235     mSvtHitNtuple   = 0;
00236     mSsdHitNtuple   = 0;
00237 }
00238 
00239 //_________________________________________________
00240 StMcAnalysisMaker::~StMcAnalysisMaker()
00241 {
00242     //  StMcAnalysisMaker Destructor
00243     //  delete the histograms
00244     cout << "Inside StMcAnalysisMaker Destructor" << endl;
00245     SafeDelete(mAssociationCanvas);
00246 //     SafeDelete(mMomResolution);
00247 //     SafeDelete(mHitResolution);
00248 //     SafeDelete(coordRec);
00249 //     SafeDelete(coordMcPartner);
00250 //     SafeDelete(mTrackNtuple);
00251 }
00252 
00253 //_____________________________________________________________________________
00254 
00255 void StMcAnalysisMaker::Clear(const char*)
00256 {
00257     // StMcAnalysisMaker - Clear,
00258     // Don't delete the canvas, so that it stays if needed
00259     
00260     StMaker::Clear();
00261 }
00262 
00263 //_________________________________________________
00264 Int_t StMcAnalysisMaker::Finish()
00265 {
00266   if (mNtupleFile) {
00267     mNtupleFile->Write();
00268     mNtupleFile->Close();
00269   }
00270     return StMaker::Finish();
00271 }
00272 
00273 
00274 //_________________________________________________
00275 Int_t StMcAnalysisMaker::Init()
00276 {
00277     // StMcAnalysisMaker - Init
00278     SetZones();  // This is my method to set the zones for the canvas.
00279 
00280     mNtupleFile = GetTFile();
00281     if (mNtupleFile) {mNtupleFile->cd(); mNtupleFile = 0;}
00282     else {mNtupleFile = new TFile("TrackMapNtuple.root","RECREATE","Track Ntuple");}
00283     // Book Histograms Here so they can be found and deleted by Victor's chain (I hope).
00284     mHitResolution = new TH2F("hitRes","Delta Z Vs Delta X for Hits",
00285                              mNumDeltaX,mMinDeltaX,mMaxDeltaX,mNumDeltaZ,mMinDeltaZ,mMaxDeltaZ);
00286     mHitResolution->SetXTitle("Delta X (cm)");
00287     mHitResolution->SetYTitle("Delta Z (cm)");
00288     mSvtHitResolution = new TH2F("SvtHitRes","Delta Z Vs Delta X for SvtHits",
00289                              mNumDeltaX,mMinDeltaX,mMaxDeltaX,mNumDeltaZ,mMinDeltaZ,mMaxDeltaZ);
00290     mSvtHitResolution->SetXTitle("Delta X (cm)");
00291     mSvtHitResolution->SetYTitle("Delta Z (cm)");
00292 
00293     mSsdHitResolution = new TH2F("SsdHitRes","Delta Z Vs Delta X for SsdHits",
00294                              mNumDeltaX,mMinDeltaX,mMaxDeltaX,mNumDeltaZ,mMinDeltaZ,mMaxDeltaZ);
00295     mSsdHitResolution->SetXTitle("Delta X (cm)");
00296     mSsdHitResolution->SetYTitle("Delta Z (cm)");
00297 
00298     mMomResolution = new TH1F("momRes","(|p| - |pmc|)/|pmc|",100,-1.,1.);
00299     mMomResolution->SetXTitle("Resolution (%)");
00300 
00301     coordRec = new TH2F("coordRc","X vs Y pos. of Hits", 100, -200, 200, 100, -200, 200);
00302     coordRec->SetXTitle("X (cm)");
00303     coordRec->SetYTitle("Y (cm)");
00304     
00305     coordMcPartner = new TH2F("coordMc","X vs Y pos. of Hits", 100, -200, 200, 100, -200, 200);
00306     coordMcPartner->SetXTitle("X (cm)");
00307     coordMcPartner->SetYTitle("Y (cm)");
00308 
00309     // Define the file for the Ntuple, otherwise it won't be available later.
00310     // one must define the file _after_ the histograms are booked, otherwise they are
00311     // not owned by the maker, but are stored in the file, breaking the code in StAssociator.
00312 
00313     mNtupleFile = new TFile("TrackMapNtuple.root","RECREATE","Track Ntuple");
00314     
00315     const char* vars = "px:py:pz:p:pxrec:pyrec:pzrec:prec:commTpcHits:hitDiffX:hitDiffY:hitDiffZ:mcTrkId:mostCommIdTruth:nHitsIdTruth:nMcHits:nFitPts:nDetPts:quality";
00316     mTrackNtuple = new TNtuple("TrackNtuple","Track Pair Info",vars);
00317     mTrackNtuple->SetAutoSave(100000000);
00318     if (! m_Mode || m_Mode & 0x1) {
00319       mTpcHitNtuple = new TNtuple("TpcHitNtuple","the TPC hit pairs Info",vTpcHitMRPair);
00320       mTpcHitNtuple->SetAutoSave(100000000);
00321     }
00322     if (! m_Mode || m_Mode & 0x2) {
00323       mSvtHitNtuple = new TNtuple("SvtHitNtuple","the SVT hit pairs Info",vSvtHitMRPair);
00324       mSvtHitNtuple->SetAutoSave(100000000);
00325     }
00326     if (! m_Mode || m_Mode & 0x4) {
00327       mSsdHitNtuple = new TNtuple("SsdHitNtuple","the SSD hit pairs Info",vSsdHitMRPair);
00328       mSsdHitNtuple->SetAutoSave(100000000);
00329     }
00330     return StMaker::Init();
00331 }
00332 //_________________________________________________
00333 Int_t StMcAnalysisMaker::Make()
00334 {
00335   // Get the pointers we need, we have to use the titles we gave them in the
00336   // macro.  I just used the defaults.
00337   
00338   // StEvent
00339   StEvent* rEvent =  (StEvent*) GetInputDS("StEvent");
00340   const StMcTrack  *mTrack = 0;
00341   
00342   // StMcEvent
00343   StMcEvent* mEvent = (StMcEvent*) GetDataSet("StMcEvent");
00344   
00345   // StAssociationMaker
00346   StAssociationMaker* assoc = 0;
00347   assoc = (StAssociationMaker*) GetMaker("StAssociationMaker");
00348   
00349   // the Multimaps...
00350   rcTpcHitMapType* theHitMap   = assoc->rcTpcHitMap();
00351   mcTpcHitMapType* theMcHitMap = assoc->mcTpcHitMap();
00352   rcSvtHitMapType* svtHitMap   = assoc->rcSvtHitMap();
00353   mcSvtHitMapType* svtMcHitMap = assoc->mcSvtHitMap();
00354   rcSsdHitMapType* ssdHitMap   = assoc->rcSsdHitMap();
00355   mcSsdHitMapType* ssdMcHitMap = assoc->mcSsdHitMap();
00356   rcTrackMapType*  theTrackMap = assoc->rcTrackMap();
00357   mcV0MapType* theMcV0Map      = assoc->mcV0Map(); 
00358   
00359   if (!theHitMap) {
00360     gMessMgr->Warning() << "----------WARNING----------\n"
00361                         << "No Hit Map found for this event!" << endm;
00362     return kStWarn;
00363   }
00364   // Example: look at the position of the primary vertex
00365   //          Map is not needed for this, but it's a good check,
00366   //          tracking will not be good if primary vertex was not well placed.
00367   
00368   // First check whether the Primary Vertex is there at all.
00369   StThreeVectorD VertexPos(0,0,0);
00370   if (rEvent->primaryVertex()) {
00371     VertexPos = rEvent->primaryVertex()->position();
00372     cout << "Position of Primary Vertex from StEvent:" << endl;
00373     cout << VertexPos << endl;
00374   }
00375   else {
00376     cout << "----------- WARNING ------------" << endl;
00377     cout << "No primary vertex from StEvent!!" << endl;
00378     cout << "Assume it is at (0,0,0)         " << endl;
00379     
00380   }
00381   cout << "Position of Primary Vertex from StMcEvent:" << endl;
00382   cout << mEvent->primaryVertex()->position() << endl;
00383   
00384   // Example: look at hits associated with 1st REC hit in Tpc Hit collection.
00385   
00386   StTpcHit*     firstHit;
00387   Bool_t        gotOneHit;
00388   StTpcHitCollection* tpcColl = rEvent->tpcHitCollection();
00389   unsigned int j,k, nhits;
00390   nhits = tpcColl->numberOfHits();
00391   if (tpcColl && nhits) {
00392     gotOneHit = kFALSE;
00393     memset (&TpcHitMRPair, 0, sizeof(TpcHitMRPair_t));
00394     for (k=0; !gotOneHit && k<tpcColl->numberOfSectors(); k++)
00395       for (j=0; !gotOneHit && j<tpcColl->sector(k)->numberOfPadrows(); j++)
00396         if (tpcColl->sector(k)->padrow(j)->hits().size()) {
00397           
00398           firstHit = tpcColl->sector(k)->padrow(j)->hits()[0];
00399           cout << "First Hit Found in Sector "
00400                << firstHit->sector() << " Padrow "
00401                << firstHit->padrow() << endl;
00402           gotOneHit = kTRUE;
00403         }
00404     cout << "This hit has " <<  theHitMap->count(firstHit) << " MC Hits associated with it."<< endl;
00405     // To get the associated hits of the first hit we use equal_range(key), which returns
00406     // 2 iterators, the lower bound and upper bound, so that then we can loop over them.
00407     
00408     cout << "Position of First Rec. Hit and Associated (if any) MC Hit:" << endl;
00409     pair<rcTpcHitMapIter,rcTpcHitMapIter> hitBounds = theHitMap->equal_range(firstHit);    
00410     for (rcTpcHitMapIter it=hitBounds.first; it!=hitBounds.second; ++it) 
00411       cout << "[" << (*it).first->position() << ", " << (*it).second->position() << "]" << endl;
00412   }
00413   else {
00414     cout << "There are no reconstructed TPC Hits in this event" << endl;
00415   }
00416   // Example: Make a histogram using the Hit Map.
00417   // Fill histogram from map
00418   
00419   float DeltaX;
00420   float DeltaZ;
00421   if (theHitMap && mTpcHitNtuple) {// TPC
00422     StTpcHitCollection* recHits = rEvent->tpcHitCollection();
00423     StMcTpcHitCollection* mcHits = mEvent->tpcHitCollection();
00424     assert (recHits || mcHits);
00425     cout << "Making Hit Resolution Histogram..." << endl;
00426     // Loop over Rec Hits
00427     
00428     for (unsigned int iSector=0; iSector< recHits->numberOfSectors(); iSector++) {
00429       for (unsigned int iPadrow=0; iPadrow<recHits->sector(iSector)->numberOfPadrows();
00430            iPadrow++) {
00431         for (StTpcHitIterator iter = recHits->sector(iSector)->padrow(iPadrow)->hits().begin();
00432              iter != recHits->sector(iSector)->padrow(iPadrow)->hits().end();
00433              iter++) {
00434           const StTpcHit   *rhit = dynamic_cast<const StTpcHit   *> (*iter);
00435           assert(rhit);
00436           if (rhit->TestBit(StMcHit::kMatched)) 
00437           {
00438             pair<rcTpcHitMapIter,rcTpcHitMapIter>
00439               recBounds = theHitMap->equal_range(rhit);
00440             
00441             for (rcTpcHitMapIter it2=recBounds.first; it2!=recBounds.second; ++it2){
00442               
00443               const StMcTpcHit *mhit = dynamic_cast<const StMcTpcHit *> ((*it2).second);
00444               assert ( mhit);
00445               DeltaX = rhit->position().x() - mhit->position().x();
00446               DeltaZ = rhit->position().z() - mhit->position().z();
00447               
00448               mHitResolution->Fill(DeltaX,DeltaZ);
00449               
00450               memset (&TpcHitMRPair, 0, sizeof(TpcHitMRPair));
00451               TpcHitMRPair.sector   = rhit->sector();
00452               TpcHitMRPair.row      = rhit->padrow();
00453               TpcHitMRPair.xR       = rhit->position().x();
00454               TpcHitMRPair.yR       = rhit->position().y();
00455               TpcHitMRPair.zR       = rhit->position().z();
00456               TpcHitMRPair.dER      = rhit->charge();
00457               TpcHitMRPair.IdR      = rhit->idTruth();
00458               TpcHitMRPair.qR       = rhit->qaTruth();
00459               TpcHitMRPair.nR       = theHitMap->count(rhit);
00460               TpcHitMRPair.isDet    = mhit->isDet();
00461               TpcHitMRPair.xM       = mhit->position().x();
00462               TpcHitMRPair.yM       = mhit->position().y();
00463               TpcHitMRPair.zM       = mhit->position().z();
00464               TpcHitMRPair.pxM      = mhit->localMomentum().x();
00465               TpcHitMRPair.pyM      = mhit->localMomentum().y();
00466               TpcHitMRPair.pzM      = mhit->localMomentum().z();
00467               TpcHitMRPair.dEM      = mhit->dE();
00468               TpcHitMRPair.dSM      = mhit->dS();
00469               mTrack     = mhit->parentTrack();
00470               if (mTrack) TpcHitMRPair.IdM = mTrack->key();
00471               else        TpcHitMRPair.IdM = 0;
00472               TpcHitMRPair.nM     = theMcHitMap->count(mhit);
00473               mTpcHitNtuple->Fill(&TpcHitMRPair.sector);
00474             }
00475           }
00476           else {
00477             memset (&TpcHitMRPair, 0, sizeof(TpcHitMRPair));
00478             TpcHitMRPair.sector   = rhit->sector();
00479             TpcHitMRPair.row      = rhit->padrow();
00480             TpcHitMRPair.xR       = rhit->position().x();
00481             TpcHitMRPair.yR       = rhit->position().y();
00482             TpcHitMRPair.zR       = rhit->position().z();
00483             TpcHitMRPair.dER      = rhit->charge();
00484             TpcHitMRPair.IdR      = rhit->idTruth();
00485             TpcHitMRPair.qR       = rhit->qaTruth();
00486             TpcHitMRPair.nR       = theHitMap->count(rhit);
00487             mTpcHitNtuple->Fill(&TpcHitMRPair.sector);
00488           }
00489         }
00490       }
00491     }
00492     for (unsigned int iSector=0;
00493          iSector<mcHits->numberOfSectors(); iSector++) {
00494       
00495       if (Debug()) {cout << iSector + 1 << " "; flush(cout);}
00496       StMcTpcSectorHitCollection* tpcSectHitColl = mcHits->sector(iSector);
00497       for (unsigned int iPadrow=0;
00498            iPadrow<tpcSectHitColl->numberOfPadrows();
00499            iPadrow++) {
00500         StMcTpcPadrowHitCollection* tpcPadRowHitColl = tpcSectHitColl->padrow(iPadrow);
00501         for (unsigned int iHit=0;
00502              iHit<tpcPadRowHitColl->hits().size();
00503              iHit++){
00504           const StMcTpcHit *mhit = dynamic_cast<const StMcTpcHit *> (tpcPadRowHitColl->hits()[iHit]);
00505           assert (mhit);
00506           if (mhit->TestBit(StMcHit::kMatched)) continue;
00507           memset (&TpcHitMRPair, 0, sizeof(TpcHitMRPair));
00508           TpcHitMRPair.sector   = mhit->sector();
00509           TpcHitMRPair.row      = mhit->padrow();
00510           TpcHitMRPair.isDet    = mhit->isDet();
00511           TpcHitMRPair.xM       = mhit->position().x();
00512           TpcHitMRPair.yM       = mhit->position().y();
00513           TpcHitMRPair.zM       = mhit->position().z();
00514           TpcHitMRPair.pxM      = mhit->localMomentum().x();
00515           TpcHitMRPair.pyM      = mhit->localMomentum().y();
00516           TpcHitMRPair.pzM      = mhit->localMomentum().z();
00517           TpcHitMRPair.dEM      = mhit->dE();
00518           TpcHitMRPair.dSM      = mhit->dS();
00519           mTrack     = mhit->parentTrack();
00520           if (mTrack) TpcHitMRPair.IdM = mTrack->key();
00521           else        TpcHitMRPair.IdM = 0;
00522           mTpcHitNtuple->Fill(&TpcHitMRPair.sector);
00523         }
00524       }
00525     }
00526   }
00527   if (svtHitMap && svtMcHitMap && mSvtHitNtuple) {  // svt hits
00528     StSvtHitCollection* recHits = rEvent->svtHitCollection();
00529     StMcSvtHitCollection* mcHits = mEvent->svtHitCollection();
00530     if (recHits && mcHits) {
00531       for (unsigned int iBarrel=0; iBarrel< recHits->numberOfBarrels(); iBarrel++) {
00532         for (unsigned int iLadder=0; iLadder<recHits->barrel(iBarrel)->numberOfLadders(); iLadder++) {
00533           for (unsigned int iWafer = 0; iWafer < recHits->barrel(iBarrel)->ladder(iLadder)->numberOfWafers(); iWafer++) {
00534             for (StSvtHitIterator iter = recHits->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().begin();
00535                iter != recHits->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().end();
00536                  iter++) {
00537               const StSvtHit   *rhit = dynamic_cast<const StSvtHit   *> (*iter);
00538               assert(rhit);
00539               if (rhit->TestBit(StMcHit::kMatched)) 
00540               {
00541                 pair<rcSvtHitMapIter,rcSvtHitMapIter> recBounds = svtHitMap->equal_range(rhit);
00542                 for (rcSvtHitMapIter it2=recBounds.first; it2!=recBounds.second; ++it2){
00543                   const StMcSvtHit *mhit = dynamic_cast<const StMcSvtHit *> ((*it2).second);
00544                   assert ( mhit);
00545                   DeltaX = rhit->position().x() - mhit->position().x();
00546                   DeltaZ = rhit->position().z() - mhit->position().z();
00547                   mSvtHitResolution->Fill(DeltaX,DeltaZ);
00548                   memset (&SvtHitMRPair, 0, sizeof(SvtHitMRPair));
00549                   SvtHitMRPair.barrel   = rhit->barrel();
00550                   SvtHitMRPair.ladder   = rhit->ladder();
00551                   SvtHitMRPair.layer    = rhit->layer();
00552                   SvtHitMRPair.wafer    = rhit->wafer();
00553                   SvtHitMRPair.hybrid   = rhit->hybrid();
00554                   SvtHitMRPair.index    = rhit->index();
00555                   SvtHitMRPair.xR       = rhit->position().x();
00556                   SvtHitMRPair.yR       = rhit->position().y();
00557                   SvtHitMRPair.zR       = rhit->position().z();
00558                   SvtHitMRPair.dER      = rhit->charge();
00559                   SvtHitMRPair.IdR      = rhit->idTruth();
00560                   SvtHitMRPair.qR       = rhit->qaTruth();
00561                   SvtHitMRPair.nR       = svtHitMap->count(rhit);
00562                   SvtHitMRPair.xM       = mhit->position().x();
00563                   SvtHitMRPair.yM       = mhit->position().y();
00564                   SvtHitMRPair.zM       = mhit->position().z();
00565                   SvtHitMRPair.pxM      = mhit->localMomentum().x();
00566                   SvtHitMRPair.pyM      = mhit->localMomentum().y();
00567                   SvtHitMRPair.pzM      = mhit->localMomentum().z();
00568                   SvtHitMRPair.dEM      = mhit->dE();
00569                   SvtHitMRPair.dSM      = mhit->dS();
00570                   mTrack     = mhit->parentTrack();
00571                   if (mTrack) SvtHitMRPair.IdM = mTrack->key();
00572                   else        SvtHitMRPair.IdM = 0;
00573                   SvtHitMRPair.nM     = svtMcHitMap->count(mhit);
00574                   mSvtHitNtuple->Fill(&SvtHitMRPair.barrel);
00575                 }
00576               }
00577               else {
00578                 memset (&SvtHitMRPair, 0, sizeof(SvtHitMRPair));
00579                 SvtHitMRPair.barrel   = rhit->barrel();
00580                 SvtHitMRPair.ladder   = rhit->ladder();
00581                 SvtHitMRPair.layer    = rhit->layer();
00582                 SvtHitMRPair.wafer    = rhit->wafer();
00583                 SvtHitMRPair.hybrid   = rhit->hybrid();
00584                 SvtHitMRPair.index    = rhit->index();
00585                 SvtHitMRPair.xR       = rhit->position().x();
00586                 SvtHitMRPair.yR       = rhit->position().y();
00587                 SvtHitMRPair.zR       = rhit->position().z();
00588                 SvtHitMRPair.dER      = rhit->charge();
00589                 SvtHitMRPair.IdR      = rhit->idTruth();
00590                 SvtHitMRPair.qR       = rhit->qaTruth();
00591                 SvtHitMRPair.nR       = svtHitMap->count(rhit);
00592                 mSvtHitNtuple->Fill(&SvtHitMRPair.barrel);
00593               }
00594             }
00595           }
00596         }
00597       }
00598       for (unsigned int iBarrel=0; iBarrel< mcHits->numberOfBarrels(); iBarrel++) {
00599         for (unsigned int iLadder=0; iLadder<mcHits->barrel(iBarrel)->numberOfLadders(); iLadder++) {
00600           for (unsigned int iWafer = 0; iWafer < mcHits->barrel(iBarrel)->ladder(iLadder)->numberOfWafers(); iWafer++) {
00601             for (StMcSvtHitIterator iter = mcHits->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().begin();
00602                iter != mcHits->barrel(iBarrel)->ladder(iLadder)->wafer(iWafer)->hits().end();
00603                  iter++) {
00604               const StMcSvtHit   *mhit = dynamic_cast<const StMcSvtHit   *> (*iter);
00605               assert (mhit);
00606               if (mhit->TestBit(StMcHit::kMatched)) continue;
00607               memset (&SvtHitMRPair, 0, sizeof(SvtHitMRPair));
00608               SvtHitMRPair.barrel   = mhit->barrel();
00609               SvtHitMRPair.ladder   = mhit->ladder();
00610               SvtHitMRPair.layer    = mhit->layer();
00611               SvtHitMRPair.wafer    = mhit->wafer();
00612               SvtHitMRPair.hybrid   = mhit->hybrid();
00613               SvtHitMRPair.index    = -1;
00614               SvtHitMRPair.barrel   = mhit->barrel();
00615               SvtHitMRPair.ladder   = mhit->ladder();
00616               SvtHitMRPair.xM       = mhit->position().x();
00617               SvtHitMRPair.yM       = mhit->position().y();
00618               SvtHitMRPair.zM       = mhit->position().z();
00619               SvtHitMRPair.pxM      = mhit->localMomentum().x();
00620               SvtHitMRPair.pyM      = mhit->localMomentum().y();
00621               SvtHitMRPair.pzM      = mhit->localMomentum().z();
00622               SvtHitMRPair.dEM      = mhit->dE();
00623               SvtHitMRPair.dSM      = mhit->dS();
00624               mTrack     = mhit->parentTrack();
00625               if (mTrack) SvtHitMRPair.IdM = mTrack->key();
00626               else        SvtHitMRPair.IdM = 0;
00627               mSvtHitNtuple->Fill(&SvtHitMRPair.barrel);
00628             }
00629           }
00630         }
00631       }
00632     }
00633   }
00634   if (ssdHitMap && ssdMcHitMap && mSsdHitNtuple) {  // ssd hits
00635     StSsdHitCollection* recHits = rEvent->ssdHitCollection();
00636     StMcSsdHitCollection* mcHits = mEvent->ssdHitCollection();
00637     if (recHits && mcHits) {
00638       for (unsigned int iLadder=0; iLadder<recHits->numberOfLadders(); iLadder++) {
00639         for (unsigned int iWafer = 0; iWafer < recHits->ladder(iLadder)->numberOfWafers(); iWafer++) {
00640           for (StSsdHitIterator iter = recHits->ladder(iLadder)->wafer(iWafer)->hits().begin();
00641                iter != recHits->ladder(iLadder)->wafer(iWafer)->hits().end();
00642                iter++) {
00643             const StSsdHit   *rhit = dynamic_cast<const StSsdHit   *> (*iter);
00644             assert(rhit);
00645             if (rhit->TestBit(StMcHit::kMatched)) 
00646               {
00647                 pair<rcSsdHitMapIter,rcSsdHitMapIter> recBounds = ssdHitMap->equal_range(rhit);
00648                 for (rcSsdHitMapIter it2=recBounds.first; it2!=recBounds.second; ++it2){
00649                   const StMcSsdHit *mhit = dynamic_cast<const StMcSsdHit *> ((*it2).second);
00650                   assert ( mhit);
00651                   DeltaX = rhit->position().x() - mhit->position().x();
00652                   DeltaZ = rhit->position().z() - mhit->position().z();
00653                   mSsdHitResolution->Fill(DeltaX,DeltaZ);
00654                   memset (&SsdHitMRPair, 0, sizeof(SsdHitMRPair));
00655                   SsdHitMRPair.ladder   = rhit->ladder();
00656                   SsdHitMRPair.wafer    = rhit->wafer();
00657                   SsdHitMRPair.xR       = rhit->position().x();
00658                   SsdHitMRPair.yR       = rhit->position().y();
00659                   SsdHitMRPair.zR       = rhit->position().z();
00660                   SsdHitMRPair.dER      = rhit->charge();
00661                   SsdHitMRPair.IdR      = rhit->idTruth();
00662                   SsdHitMRPair.qR       = rhit->qaTruth();
00663                   SsdHitMRPair.nR       = ssdHitMap->count(rhit);
00664                   SsdHitMRPair.xM       = mhit->position().x();
00665                   SsdHitMRPair.yM       = mhit->position().y();
00666                   SsdHitMRPair.zM       = mhit->position().z();
00667                   SsdHitMRPair.pxM      = mhit->localMomentum().x();
00668                   SsdHitMRPair.pyM      = mhit->localMomentum().y();
00669                   SsdHitMRPair.pzM      = mhit->localMomentum().z();
00670                   SsdHitMRPair.dEM      = mhit->dE();
00671                   SsdHitMRPair.dSM      = mhit->dS();
00672                   mTrack     = mhit->parentTrack();
00673                   if (mTrack) SsdHitMRPair.IdM = mTrack->key();
00674                   else        SsdHitMRPair.IdM = 0;
00675                   SsdHitMRPair.nM     = ssdMcHitMap->count(mhit);
00676                   mSsdHitNtuple->Fill(&SsdHitMRPair.ladder);
00677                 }
00678               }
00679             else {
00680               memset (&SsdHitMRPair, 0, sizeof(SsdHitMRPair));
00681               SsdHitMRPair.ladder   = rhit->ladder();
00682               SsdHitMRPair.wafer    = rhit->wafer();
00683               SsdHitMRPair.xR       = rhit->position().x();
00684               SsdHitMRPair.yR       = rhit->position().y();
00685               SsdHitMRPair.zR       = rhit->position().z();
00686               SsdHitMRPair.dER      = rhit->charge();
00687               SsdHitMRPair.IdR      = rhit->idTruth();
00688               SsdHitMRPair.qR       = rhit->qaTruth();
00689               SsdHitMRPair.nR       = ssdHitMap->count(rhit);
00690               mSsdHitNtuple->Fill(&SsdHitMRPair.ladder);
00691             }
00692           }
00693         }
00694       }
00695       for (unsigned int iLadder=0; iLadder<mcHits->numberOfLadders(); iLadder++) {
00696         for (unsigned int iWafer = 0; iWafer < mcHits->ladder(iLadder)->numberOfWafers(); iWafer++) {
00697           for (StMcSsdHitIterator iter = mcHits->ladder(iLadder)->wafer(iWafer)->hits().begin();
00698                iter != mcHits->ladder(iLadder)->wafer(iWafer)->hits().end();
00699                iter++) {
00700             const StMcSsdHit   *mhit = dynamic_cast<const StMcSsdHit   *> (*iter);
00701             assert (mhit);
00702             if (mhit->TestBit(StMcHit::kMatched)) continue;
00703             memset (&SsdHitMRPair, 0, sizeof(SsdHitMRPair));
00704             SsdHitMRPair.ladder   = mhit->ladder();
00705             SsdHitMRPair.wafer    = mhit->wafer();
00706             SsdHitMRPair.ladder   = mhit->ladder();
00707             SsdHitMRPair.xM       = mhit->position().x();
00708             SsdHitMRPair.yM       = mhit->position().y();
00709             SsdHitMRPair.zM       = mhit->position().z();
00710             SsdHitMRPair.pxM      = mhit->localMomentum().x();
00711             SsdHitMRPair.pyM      = mhit->localMomentum().y();
00712             SsdHitMRPair.pzM      = mhit->localMomentum().z();
00713             SsdHitMRPair.dEM      = mhit->dE();
00714             SsdHitMRPair.dSM      = mhit->dS();
00715             mTrack     = mhit->parentTrack();
00716             if (mTrack) SsdHitMRPair.IdM = mTrack->key();
00717             else        SsdHitMRPair.IdM = 0;
00718             mSsdHitNtuple->Fill(&SsdHitMRPair.ladder);
00719           }
00720         }
00721       }
00722     }
00723   }
00724   cout << "Finished Making Histogram." << endl;
00725   
00726   if (!theTrackMap) {
00727     gMessMgr->Warning() << "----------WARNING----------\n"
00728                         << "No Track Map found for this event!" << endm;
00729     return kStWarn;
00730   }
00731   // Example: look at the magnitude of the momentum of
00732   //          the MC track associated with first track in track Collection
00733   
00734   StSPtrVecTrackNode& rcTrackNodes = rEvent->trackNodes();
00735   StTrackNode*        firstTrackNode = 0;
00736   const StGlobalTrack*      firstTrack = 0;
00737   if (! rcTrackNodes.size()) {
00738     cout << "Track Nodes List is empty!" << endl;
00739     return kStWarn;
00740   }
00741   firstTrackNode = *(rcTrackNodes.begin());
00742   if (! firstTrackNode) {
00743     cout << "No tracks in Track Nodes List!" << endl;
00744     return kStWarn;
00745   }
00746   firstTrack = dynamic_cast<const StGlobalTrack*>(firstTrackNode->track(global));
00747   if (! firstTrack) {
00748     cout << "GlobalTrack for first Track Nodehas not been found!" << endl;
00749     return kStWarn;
00750   }
00751   cout << "MC Tracks associated with first Track in collection: " << theTrackMap->count(firstTrack) << endl;
00752   if (!theTrackMap->count(firstTrack) && theTrackMap->size()) {
00753       cout << "First track in track node container was not associated.  Pick first track in map." << endl;
00754       firstTrack = theTrackMap->begin()->first; //first entry in the map is a pair, and first entry of pair is the global track
00755   }
00756   pair<rcTrackMapIter,rcTrackMapIter> trackBounds = theTrackMap->equal_range(firstTrack);
00757   cout << "Momentum of First Track and of Associated Tracks (if there are any):" << endl;
00758   // Get the momentum of the track and compare it to MC Track.
00759   // Use primary track if available
00760   const StPrimaryTrack* pTrk = dynamic_cast<const StPrimaryTrack*>(firstTrack->node()->track(primary));
00761   StThreeVectorD recMom(0,0,0);
00762   if (pTrk)
00763     recMom = pTrk->geometry()->momentum();
00764   else
00765     recMom = firstTrack->geometry()->momentum();
00766   for (rcTrackMapIter trkIt=trackBounds.first; trkIt!=trackBounds.second; ++trkIt) { 
00767     const StMcTrack*   origTrk = (*trkIt).second->partnerMcTrack();
00768     cout << "[" << abs(recMom) << ", ";
00769     cout << abs((*trkIt).second->partnerMcTrack()->momentum()) << "]" << endl;
00770     cout << "These tracks have : \n";
00771     cout << (*trkIt).second->commonTpcHits() << " TPC  hits in common, out of " << origTrk->tpcHits().size() << endl;
00772     cout << (*trkIt).second->commonSvtHits() << " SVT  hits in common, out of " << origTrk->svtHits().size() << endl;
00773     cout << (*trkIt).second->commonFtpcHits() <<" FTPC hits in common, out of " << origTrk->ftpcHits().size() << endl;
00774   }
00775   // Example: Make a histogram of the momentum resolution of the event
00776   //          Make an Ntuple with rec & monte carlo mom, mean hit difference, and # of common hits
00777   const StGlobalTrack* recTrack;
00778   const StPrimaryTrack* primTrk;
00779   const StMcTrack*     mcTrack;
00780   StThreeVectorD p(0,0,0);
00781   StThreeVectorD pmc(0,0,0);
00782   float diff =0;
00783   
00784   float* values = new float[19];
00785   
00786   for (rcTrackMapIter tIter=theTrackMap->begin();
00787        tIter!=theTrackMap->end(); ++tIter){
00788     
00789     recTrack = (*tIter).first;
00790     //yf    if ((*tIter).second->commonTpcHits()<10) continue;
00791     mcTrack = (*tIter).second->partnerMcTrack();
00792     pmc = mcTrack->momentum();
00793     for (int k=0; k<3; k++) values[k] = pmc[k];          
00794     values[3]=pmc.mag();         
00795     
00796     primTrk = dynamic_cast<const StPrimaryTrack*>(recTrack->node()->track(primary));     
00797     if (primTrk)         
00798         p = primTrk->geometry()->momentum();     
00799     else         
00800         p = recTrack->geometry()->momentum();    
00801     
00802     for (int j=0; j<3; j++) values[j+4] = p[j];          
00803     values[7]=p.mag();   
00804     values[8]=(*tIter).second->commonTpcHits();
00805     // Fill 1d Mom. resolution Histogram
00806     
00807     diff = (p.mag() - pmc.mag())/pmc.mag();
00808     mMomResolution->Fill(diff,1.);
00809     // Loop to get Mean hit position diff.
00810     StThreeVectorF rHitPos(0,0,0);
00811     StThreeVectorF mHitPos(0,0,0);
00812     
00813     StPtrVecHit recTpcHits = recTrack->detectorInfo()->hits(kTpcId);
00814     multimap<int,float> idTruths;
00815     set<int> uniqueIdTruths;
00816     for (StHitIterator hi=recTpcHits.begin();
00817          hi!=recTpcHits.end(); hi++) {
00818       StHit* hit = *hi;
00819       StTpcHit* rHit = dynamic_cast<StTpcHit*>(hit);
00820       if (!rHit) { cout << "This Hit is not a TPC Hit"<< endl; continue;}
00821       idTruths.insert( multimap<int,float>::value_type(rHit->idTruth(),rHit->qaTruth()));
00822       uniqueIdTruths.insert(static_cast<int>(rHit->idTruth()));
00823       pair<rcTpcHitMapIter,rcTpcHitMapIter> rBounds = theHitMap->equal_range(rHit);
00824       for (rcTpcHitMapIter hIter=rBounds.first; hIter!=rBounds.second; hIter++) {
00825         const StMcTpcHit* mHit = (*hIter).second;
00826         if (mHit->parentTrack() != mcTrack) continue;
00827         rHitPos += rHit->position();
00828         mHitPos += mHit->position();
00829       }// Associated Hits Loop.
00830       
00831     } // Hits of rec. Track Loop
00832     
00833     rHitPos /=(float) (*tIter).second->commonTpcHits();
00834     mHitPos /=(float) (*tIter).second->commonTpcHits();
00835     for (int jj=0; jj<3; jj++) values[9+jj] = rHitPos[jj] - mHitPos[jj];
00836     values[12] = mcTrack->key();
00837     // Figure out the most common IdTruth; the dominatrix track!
00838     int mostCommonIdTruth = -9; 
00839     int cachedNHitsIdTruth = 0;
00840     for (set<int>::iterator si=uniqueIdTruths.begin(); si!=uniqueIdTruths.end(); ++si) {
00841         int currentNHitsIdTruth = idTruths.count(static_cast<int>(*si));
00842         if (currentNHitsIdTruth>cachedNHitsIdTruth) {
00843             mostCommonIdTruth = *si; 
00844             cachedNHitsIdTruth = currentNHitsIdTruth;
00845         }
00846     }
00847     // at this point we know the most common IdTruth,
00848     // now calculate the track "quality" for this track, averaging
00849     // the hit qualities
00850     float idQuality = 0;
00851     
00852     pair<multimap<int,float>::iterator,multimap<int,float>::iterator> mostCommRange = idTruths.equal_range(mostCommonIdTruth);
00853     for (multimap<int,float>::iterator mi=mostCommRange.first; mi!=mostCommRange.second; ++mi) {
00854         idQuality+=mi->second;
00855     }
00856     idQuality/=cachedNHitsIdTruth;
00857     
00858     values[13] = mostCommonIdTruth;
00859     values[14] = cachedNHitsIdTruth;
00860     values[15] = mcTrack->tpcHits().size();
00861     values[16] = recTrack->fitTraits().numberOfFitPoints(kTpcId);
00862     values[17] = recTpcHits.size();
00863     values[18] = idQuality;
00864     mTrackNtuple->Fill(values);
00865   } // Tracks in Map Loop
00866   cout << "Finished Track Loop, Made Ntuple" << endl;
00867   //delete vars;
00868   delete [] values;
00869   //mNtupleFile->Write(); // Write the Ntuple to the File in Finish(), not here
00870   
00871   // Example: Make 2 Histograms
00872   // - x and y positions of the hits from the reconstructed track.
00873   // - x and y positions of the hits from the  Monte Carlo  track.
00874   
00875   unsigned int maxCommonTpcHits = 0;
00876   const StMcTrack* partner = 0;
00877   trackBounds = theTrackMap->equal_range(firstTrack);
00878   for (rcTrackMapIter rcIt = trackBounds.first;
00879        rcIt != trackBounds.second;
00880        ++rcIt) {
00881     if ((*rcIt).second->commonTpcHits() >  maxCommonTpcHits) {
00882       partner = (*rcIt).second->partnerMcTrack();
00883       maxCommonTpcHits = (*rcIt).second->commonTpcHits();
00884     }
00885   }
00886   StHitIterator rcHitIt;
00887   StMcTpcHitIterator mcHitIt;
00888   StMcSvtHitIterator mcSitIt;
00889   StPtrVecHit theHits = firstTrack->detectorInfo()->hits();
00890   for (rcHitIt  = theHits.begin();
00891        rcHitIt != theHits.end();
00892        rcHitIt++) coordRec->Fill((*rcHitIt)->position().x(),(*rcHitIt)->position().y());
00893   if (partner) {
00894     for (mcHitIt  = ((std::vector<StMcTpcHit*> *)&partner->tpcHits() )->begin();
00895          mcHitIt != partner->tpcHits().end();
00896          mcHitIt++) coordMcPartner->Fill((*mcHitIt)->position().x(),(*mcHitIt)->position().y());
00897     for (mcSitIt  = ((std::vector<StMcSvtHit*> *)&partner->svtHits())->begin();
00898          mcSitIt != partner->svtHits().end();
00899          mcSitIt++) coordMcPartner->Fill((*mcSitIt)->position().x(),(*mcSitIt)->position().y());
00900     
00901   }
00902   if (!theMcV0Map) {
00903     gMessMgr->Warning() << "----------WARNING----------\n"
00904                         << "No V0 Map found for this event!" << endm;
00905     return kStWarn;
00906   }
00907   //Example: Print out position of V0 vertices that have been associated.
00908   //         (LSB)
00909   StSPtrVecMcVertex& mcVertices = mEvent->vertices();
00910   const StV0Vertex* rcV0Partner;
00911   StMcVertexIterator mcVertexIt;
00912   
00913   //Loop over all MC vertices
00914   bool foundV0Pair = false;
00915   for (mcVertexIt = mcVertices.begin(); mcVertexIt != mcVertices.end();
00916        mcVertexIt++){
00917     // Get the upper and lower bounds.
00918     pair<mcV0MapIter,mcV0MapIter> mcV0Bounds = theMcV0Map->equal_range(*mcVertexIt);
00919     
00920     // Print out MC vertex position if there is an associated V0.
00921     
00922     if (mcV0Bounds.first != mcV0Bounds.second) {
00923       cout << "Printing Position of a V0 pair:\n";
00924       cout << "Position of MC V0 vertex: " << (*mcVertexIt)->position() << endl;
00925       foundV0Pair = true;
00926     }
00927     //Now loop over the bounds      
00928     for(mcV0MapIter mcV0MapIt = mcV0Bounds.first;
00929         mcV0MapIt != mcV0Bounds.second; ++mcV0MapIt){
00930       rcV0Partner = (*mcV0MapIt).second;
00931       cout << "Position of rc V0 vertex: " << rcV0Partner->position() << endl;
00932       
00933     }
00934     if (foundV0Pair) break; // Only print the information of 1 reconstructed vertex, to avoid a lot of output.
00935   }
00936   
00937   
00938   
00939   //mAssociationCanvas = new TCanvas("mAssociationCanvas", "Histograms",200,10,900,500);
00940   
00941   return kStOK;
00942 }