StRoot: StLaserAnalysisMaker/LaserEvent.cxx Source File

00001 //$Id: LaserEvent.cxx,v 1.7 2011/01/10 20:36:12 fisyak Exp $
00002 // $Log: LaserEvent.cxx,v $
00003 // Revision 1.7  2011/01/10 20:36:12  fisyak
00004 // Use sector/padrow in global => local transformation
00005 //
00006 // Revision 1.6  2008/06/02 13:48:02  fisyak
00007 // Add  t0 handlers for Tpx/Tpc time offsets
00008 //
00009 // Revision 1.5  2007/12/10 19:54:02  fisyak
00010 // Add Id and Log, correct spelling error in README
00011 //
00012 #include <assert.h>
00013 #include "TRandom.h"
00014 #include "TDirectory.h"
00015 #include "LaserEvent.h"
00016 #if 1
00017 #include "StProbPidTraits.h"
00018 #include "StDedxPidTraits.h"
00019 #endif
00020 #include "StVertex.h"
00021 #include "StPrimaryVertex.h"
00022 #include "StTrack.h"
00023 #include "StPrimaryTrack.h"
00024 #include "StGlobalTrack.h"
00025 #include "StTrackNode.h"
00026 #include "StTpcDb/StTpcDb.h"
00027 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00028 #include "StEventTypes.h"
00029 #include "TGeoMatrix.h"
00030 ClassImp(LaserRaft);
00031 ClassImp(EventHeader);
00032 ClassImp(LaserEvent);
00033 ClassImp(Vertex);
00034 ClassImp(Track);
00035 ClassImp(Hit);
00036 TClonesArray *LaserEvent::fgVertices = 0;
00037 TClonesArray *LaserEvent::fgTracks = 0;
00038 TClonesArray *LaserEvent::fgHits = 0;
00039 TClonesArray *LaserEvent::fgFit = 0;
00040 //______________________________________________________________________________
00041 LaserB::LaserB(const LaserRaft &laser) : Sector(laser.Sector), Raft(laser.Raft), Bundle(laser.Bundle), Mirror(laser.Mirror),
00042                                          XyzL(laser.XyzL), XyzU(laser.XyzU), XyzB(laser.XyzB), 
00043                                          dirL(laser.dirL), dirU(laser.dirU),dirB(laser.dirB),
00044                                          Theta(laser.Theta), Phi(laser.Phi){
00045 }
00046 //______________________________________________________________________________
00047 LaserEvent::LaserEvent()
00048 {
00049    // Create an LaserEvent object.
00050    // When the constructor is invoked for the first time, the class static
00051    // variable fgTracks is 0 and the TClonesArray fgTracks is created.
00052 
00053    if (!fgVertices) fgVertices = new TClonesArray("Vertex", 1000);
00054    fVertices = fgVertices;
00055    fNvertex = 0;
00056    if (!fgTracks) fgTracks = new TClonesArray("Track", 1000);
00057    fTracks = fgTracks;
00058    fNtrack = 0;
00059 
00060    if (!fgHits) fgHits = new TClonesArray("Hit", 1000);
00061    fHits = fgHits;
00062    fNhit = 0;
00063 
00064    if (!fgFit) fgFit = new TClonesArray("FitDV", 1000);
00065    fgFit->ExpandCreate(12);
00066    fFit = fgFit;
00067 }
00068 
00069 //______________________________________________________________________________
00070 LaserEvent::~LaserEvent()
00071 {
00072    Clear("C");
00073 }
00074 
00075 //______________________________________________________________________________
00076 Vertex *LaserEvent::AddVertex(StPrimaryVertex *vertex) {
00077   if (! vertex) return 0;
00078   TClonesArray &vertices = *fVertices;
00079   Vertex *vx = new(vertices[fNvertex++]) Vertex(vertex);
00080   const TGeoHMatrix &Tpc2Global = gStTpcDb->Tpc2GlobalMatrix();
00081   Tpc2Global.MasterToLocal(vx->Xyz.xyz(),vx->XyzL.xyz());
00082   return vx;
00083 }
00084 //______________________________________________________________________________
00085 Track *LaserEvent::AddTrack(Int_t sector, StTrack *track, LaserB *laser, Double_t z) {
00086   if (! track) return 0;
00087   TClonesArray &tracks = *fTracks;
00088   Track *t = new(tracks[fNtrack++]) Track(sector, track, laser, z);
00089   return t;
00090 }
00091 //______________________________________________________________________________
00092 Hit *LaserEvent::AddHit(StTpcHit *hit) {
00093   if (! hit) return 0;
00094   TClonesArray &hits = *fHits;
00095   Hit *t = new(hits[fNhit++]) Hit(hit);
00096   return t;
00097 }
00098 //________________________________________________________________________________
00099 void LaserEvent::AddTrackFit(Track *t) {
00100   if (! t) return;
00101   Int_t sector = t->Laser.Sector;
00102   Int_t bundle =  t->Laser.Bundle;
00103   Int_t mirror =  t->Laser.Mirror;
00104   Int_t s2 = (sector-1)/2;
00105   if (s2 >= 0 && s2 < 12) {
00106     TClonesArray &fits = *fFit;
00107     FitDV *fit = (FitDV *) fits[s2]; 
00108     fit->Sector = sector;
00109     Int_t N = fit->N;
00110     Double32_t x = t->Laser.XyzL.z();
00111     Double32_t y = t->XyzPL.z() - x;
00112     if (N < 42) {
00113       fit->X[N] = x;
00114       fit->Y[N] = y;
00115       fit->Bundle[N] =  bundle;
00116       fit->Mirror[N] =  mirror;
00117       fit->N   = N+1;
00118     }
00119   }
00120 }
00121 //______________________________________________________________________________
00122 void LaserEvent::Clear(Option_t *option) {
00123    fTracks->Clear(option);
00124    fHits->Clear(option);
00125    fVertices->Clear(option);
00126    fgFit->Clear(option);
00127 }
00128 //______________________________________________________________________________
00129 void LaserEvent::Reset() {
00130   // Static function to reset all static objects for this event
00131   SafeDelete(fgTracks);
00132   SafeDelete(fgHits);
00133   SafeDelete(fgVertices);
00134   SafeDelete(fgFit);
00135 }
00136 
00137 //______________________________________________________________________________
00138 void LaserEvent::SetHeader(Int_t i, Int_t run, Int_t date, Int_t time)
00139 {
00140    fNvertex = 0;
00141    fNtrack = 0;
00142    fNhit = 0;
00143    fEvtHdr.Set(i, run, date, time);
00144 }
00145 //______________________________________________________________________________
00146 void LaserEvent::SetHeader(Int_t i, Int_t run, Int_t date, Int_t time,
00147      Float_t tzero, Float_t drivel, Float_t clock)
00148 {
00149    SetHeader(i, run, date, time);
00150    fEvtHdr.SetE(tzero, drivel, clock);
00151 }
00152 //______________________________________________________________________________
00153 void LaserEvent::SetHeader(Int_t i, Int_t run, Int_t date, Int_t time,
00154      Float_t tzero, Float_t drivel, Float_t clock, Float_t trigger)
00155 {
00156    SetHeader(i, run, date, time);
00157    fEvtHdr.SetE(tzero, drivel, clock, trigger);
00158 }
00159 //______________________________________________________________________________
00160 Track::Track(Int_t sector, StTrack *track, LaserB *theLaser, Double_t z) : 
00161   Flag(0),mType(kUndefinedVtxId),  mSector(sector), 
00162   fpTInv(-999), thePath(0), dPhi(-999), dTheta(-999), zLastHit(z) {
00163   StTrackNode*   node = track->node();
00164   StGlobalTrack* gTrack = static_cast<StGlobalTrack*>(node->track(global));
00165   fgeoIn = *((StHelixModel *) gTrack->geometry());
00166   fgeoOut = *((StHelixModel *) gTrack->outerGeometry());
00167   StThreeVectorD g3 = fgeoOut.momentum();
00168   fpTInv = fgeoOut.charge()/g3.perp();
00169   fTheta = fgeoOut.momentum().theta();
00170   fPhi   = fgeoOut.momentum().phi();
00171   StPhysicalHelixD helixO = fgeoOut.helix();
00172   thePath = helixO.pathLength(Laser.XyzG.x(),Laser.XyzG.y());
00173   if (theLaser) Laser = *theLaser;
00174   if (gTrack) {
00175     StPrimaryTrack *pTrack =    static_cast<StPrimaryTrack*>(node->track(primary));
00176     if (pTrack) {
00177       StPrimaryVertex* vertex = (StPrimaryVertex*) pTrack->vertex();
00178       if (vertex) {
00179         mType = vertex->type();
00180         Vertex = vertex->position();
00181         Flag = 1;
00182       }
00183     } 
00184     mKey = gTrack->key();
00185     mFlag = gTrack->flag();
00186     mNumberOfPossiblePointsTpc = gTrack->numberOfPossiblePoints(kTpcId);
00187     mImpactParameter = gTrack->impactParameter();
00188     mLength = gTrack->length();
00189 #if 1
00190     StTrackFitTraits&  fitTraits =  gTrack->fitTraits();
00191     mNumberOfFitPointsTpc = fitTraits.numberOfFitPoints(kTpcId);
00192     mPrimaryVertexUsedInFit = fitTraits.primaryVertexUsedInFit();
00193 
00194     StSPtrVecTrackPidTraits &traits = gTrack->pidTraits();
00195     unsigned int size = traits.size();
00196     StDedxPidTraits *pid = 0;
00197     for (unsigned int i = 0; i < size; i++) {
00198       if (! traits[i]) continue;
00199       if ( traits[i]->IsZombie()) continue;
00200       pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
00201       if (pid && pid->detector() == kTpcId && pid->method() == kTruncatedMeanId) {
00202         fNdEdx = pid->numberOfPoints();      // Number of points used in dE/dx calc 
00203         fdEdx =  pid->mean();       
00204         break;
00205       }
00206     }
00207 #endif
00208   }
00209 }
00210 //______________________________________________________________________________
00211 void Track::SetPredictions(TGeoHMatrix *Raft2Tpc, TGeoHMatrix *Bundle2Tpc, TGeoHMatrix *Mirror2Tpc) {
00212   if ( Flag ) return;
00213   if ( Laser.Sector < 1 || Laser.Sector > 24) return;
00214   Flag   = 2;
00215   const TGeoHMatrix &Tpc2Global = gStTpcDb->Tpc2GlobalMatrix();
00216   StPhysicalHelixD helixO = fgeoOut.helix();
00217   thePath = helixO.pathLength(Laser.XyzG.x(),Laser.XyzG.y());
00218   XyzP   = helixO.at(thePath);
00219   Tpc2Global.MasterToLocal(XyzP.xyz(),XyzPL.xyz());
00220   StThreeVectorD g3 = fgeoOut.momentum();
00221   dirP = g3.unit();
00222   Tpc2Global.MasterToLocalVect(dirP.xyz(),dirPL.xyz());
00223   dU = StThreeVectorD(999.,999.,999.);
00224   if (Raft2Tpc) {
00225     Raft2Tpc->MasterToLocal(XyzPL.xyz(),XyzPU.xyz());
00226     Raft2Tpc->MasterToLocalVect(dirPL.xyz(),dirPU.xyz());
00227     dU = XyzPU - Laser.XyzU;
00228   }
00229   if (Bundle2Tpc) {
00230     Bundle2Tpc->MasterToLocal(XyzPL.xyz(),XyzPB.xyz());
00231     Bundle2Tpc->MasterToLocalVect(dirPL.xyz(),dirPB.xyz());
00232   }
00233   if (Mirror2Tpc) {
00234     Mirror2Tpc->MasterToLocal(XyzPL.xyz(),XyzPM.xyz());
00235     Mirror2Tpc->MasterToLocalVect(dirPL.xyz(),dirPM.xyz());
00236   }
00237   Flag += Matched();
00238 }
00239 //______________________________________________________________________________
00240 Vertex::Vertex(StPrimaryVertex *vertex) : mType(kUndefinedVtxId), WestEast(0), 
00241                                           Xyz(), numberOfDaughter(0){
00242   if (vertex) {
00243     mType = vertex->type();
00244     Xyz = vertex->position();
00245     numberOfDaughter = vertex->numberOfDaughters();
00246     if (numberOfDaughter > 24) {
00247       for (UInt_t i = 0; i < numberOfDaughter; i++) {
00248         StTrack *track = vertex->daughter(i);
00249         if (! track ) continue;
00250         StPtrVecHit hvec = track->detectorInfo()->hits();
00251         for (unsigned int j=0; j<hvec.size(); j++) {// hit loop
00252           if (hvec[j]->detector() != kTpcId) continue;
00253           StTpcHit *tpcHit = static_cast<StTpcHit *> (hvec[j]);
00254           Int_t sector = tpcHit->sector();
00255           Int_t WE = sector <= 12 ? 1 : 2;
00256           if (! WestEast) WestEast = WE;
00257           else {
00258             if (WE != WestEast) {
00259               WestEast = 0;
00260               break;
00261             }
00262           }
00263         }
00264       }
00265     }
00266   }
00267 }
00268 //________________________________________________________________________________
00269 Int_t Track::Matched() {
00270   Int_t iok = 0;
00271   dPhi = -999;
00272   dTheta = -999;
00273   iok = 1 << 10;
00274   if (Flag != 2) return iok; 
00275   iok = 10;
00276   if (Laser.Sector < 1 || Laser.Sector > 24 ||
00277       Laser.Mirror < 1 || Laser.Mirror > 7) return iok; 
00278   iok = 0;
00279   Int_t status = 0;
00280   if (TMath::Abs(dU.x()) > 0.05) status |= 1 << 3;
00281   if (TMath::Abs(dU.y()) > 0.05) status |= 1 << 4;
00282   if (thePath < 5 || thePath > 25)    status |= 1 << 5;
00283   if (mNumberOfFitPointsTpc  < 25)    status |= 1 << 6;
00284   dTheta = Laser.ThetaG-fgeoOut.dipAngle()-TMath::Pi()/2;
00285 #if 1
00286   if (TMath::Abs(dTheta) > 0.030)                  status |= 1 << 7;;
00287 #endif
00288   dPhi = Laser.PhiG - fgeoOut.psi();
00289   if (dPhi >=  TMath::Pi()) dPhi -= 2*TMath::Pi();
00290   if (dPhi <= -TMath::Pi()) dPhi += 2*TMath::Pi();
00291   //  if (SectorMirror[m][i].sigma > 0 && TMath::Abs(dPhi) > 5*SectorMirror[m][i].sigma) status |= 1 << 8;
00292 #if 1
00293   if ( TMath::Abs(dPhi) > 0.020)  status |= 1 << 9;
00294 #endif
00295   static const Double_t pTInv0 = 4.78815e-03;
00296   static const Double_t DpTInv0 = 9.75313e-03;
00297   if (TMath::Abs(fpTInv - pTInv0) > 3.0*DpTInv0) status |= 1 << 10;
00298   return iok + 10*status;
00299 }
00300 //________________________________________________________________________________
00301 Hit::Hit(StTpcHit *tpcHit)  : sector(0),row(0),charge(0),flag(0),usedInFit(0) {
00302   if (tpcHit) {
00303     sector = tpcHit->sector();
00304     row = tpcHit->padrow();;
00305     charge = tpcHit->charge();
00306     flag = tpcHit->flag();
00307     usedInFit = tpcHit->usedInFit();
00308     xyz = tpcHit->position(); // from StTpcHitMover
00309     static StTpcCoordinateTransform transform(gStTpcDb);
00310     StGlobalCoordinate glob(xyz.x(),xyz.y(),xyz.z());
00311 #if 0
00312     static StTpcLocalCoordinate  coorLT;
00313     transform(glob,coorLT);
00314     static StTpcLocalSectorAlignedCoordinate  coorLSA;
00315     transform(coorLT,coorLSA); //
00316     static StTpcLocalSectorCoordinate  local;
00317     transform(coorLSA,local); 
00318 #else
00319     static StTpcLocalSectorCoordinate  local;
00320     transform(glob,local,tpcHit->sector(),tpcHit->padrow()); 
00321 #endif
00322     xyzL = StThreeVectorF(local.position().x(),local.position().y(),local.position().z());
00323     static const Int_t NumberOfPadsAtRow[45] = {
00324       88, 96,104,112,118,126,134,142,150,158,
00325       166,174,182,  
00326       98,100,102,104,106,106,108,
00327       110,112,112,114,116,118,120,122,122,124,
00328       126,128,128,130,132,134,136,138,138,140,
00329       142,144,144,144,144
00330     };
00331     pad  = tpcHit->pad() - NumberOfPadsAtRow[row-1]/2 - 1;
00332     tbk  = tpcHit->timeBucket();
00333   }
00334 }