StRoot: StLaserAnalysisMaker/StLaserAnalysisMaker.cxx Source File

00001 // 
00002 // $Id: StLaserAnalysisMaker.cxx,v 1.16 2009/03/11 15:49:40 fisyak Exp $
00003 // $Log: StLaserAnalysisMaker.cxx,v $
00004 // Revision 1.16  2009/03/11 15:49:40  fisyak
00005 // Remove hits not beloging to primary tracks
00006 //
00007 // Revision 1.15  2008/06/02 13:48:03  fisyak
00008 // Add  t0 handlers for Tpx/Tpc time offsets
00009 //
00010 // Revision 1.14  2007/10/16 15:26:03  fisyak
00011 // Freeze version for Run VII
00012 //
00013 // Revision 1.13  2007/07/05 14:37:04  fisyak
00014 // Freeze version for run VII
00015 //
00016 // Revision 1.12  2007/05/09 13:36:44  fisyak
00017 // Freeze before intruducing Bundle coordinate system
00018 //
00019 // Revision 1.11  2007/04/17 05:12:06  perev
00020 // GetTFile()==>StMaker. Jerome request
00021 //
00022 // Revision 1.10  2007/03/06 16:31:55  fisyak
00023 // Before Selection of good runs
00024 //
00025 #define CORRECT_LASER_POSITIONS
00026 //#define __TRACKHITS__
00027 #ifndef __TRACKHITS__
00028 #define ADDPRIMTRACKHITS
00029 #endif
00030 #include <assert.h>
00031 #include "TFile.h"
00032 #include "StEventTypes.h"
00033 #include "StLaserAnalysisMaker.h"
00034 #include "laserino.h"
00035 #include "LaserEvent.h"
00036 #include "StTpcDb/StTpcDb.h"
00037 #include "TGeoMatrix.h"
00038 #include "TRVector.h"
00039 #include "TRMatrix.h"
00040 #include "TRSymMatrix.h"
00041 #include "StDetectorDbMaker/StDetectorDbClock.h"
00042 
00043 ClassImp(StLaserAnalysisMaker);
00044 static  LaserEvent             *event = 0;  
00045 static  const Int_t NS = 12;
00046 static  const Int_t NB =  6;
00047 static  const Int_t NM =  7;
00048 static LaserB *Lasers[NS][NB][NM];
00049 #if 1
00050 static Int_t NoBeams = 0;
00051 static LaserRaft *LaserBeams[NS*NB*NM];
00052 #endif
00053 static TGeoHMatrix *Traft[14];
00054 static TGeoHMatrix *Raft2Tpc[14];
00055 static TGeoHMatrix *Bundles2Tpc[14][6];
00056 static TGeoHMatrix *Mirrors2Tpc[14][6][7];
00057 //________________________________________________________________________________
00058 Int_t StLaserAnalysisMaker::Init(){
00059   event = new LaserEvent();
00060   TFile *f = GetTFile();
00061 #if 0
00062   if (! f) {
00063     f = new TFile("StLaserAnalysisMaker.root","recreate");
00064   }
00065 #endif
00066   assert(f);
00067   if (f) {
00068     f->cd();
00069     m_laser = new TTree("laser","Tpc laser track tree");
00070     m_laser->SetAutoSave(100000000); //Save every 100 MB
00071     Int_t bufsize= 64000;
00072     Int_t split = 99;
00073     if (split)  bufsize /= 4;
00074     m_laser->Branch("event", "LaserEvent",&event, bufsize, split);
00075   }
00076   const Int_t numSectors = 24;
00077   Double_t beta = 0;
00078   Double_t dBeta = 720./numSectors;
00079   Int_t sector = 0;
00080   Int_t half   = 0;
00081   TGeoHMatrix TpcHalf[2];
00082   Double_t rotHalfs[18] = {
00083     0, -1, 0, -1, 0, 0, 0, 0, -1,// sector  1-12
00084     0,  1, 0, -1, 0, 0, 0, 0,  1 // sector 13-24
00085   };
00086   for (half = 0; half <2; half++) TpcHalf[half].SetRotation(&rotHalfs[9*half]);
00087   TGeoHMatrix RotSec[24];
00088   for (sector = 1; sector <= numSectors; sector++) {
00089     if (sector > 12) beta = (numSectors-sector)*dBeta;
00090     else             beta = sector*dBeta;
00091     RotSec[sector-1].RotateZ(-beta);
00092   }
00093   memset (LaserBeams, 0, NS*NB*NM*sizeof(LaserRaft*));
00094   NoBeams = 0;
00095   memset(Traft, 0, 14*sizeof(TGeoHMatrix *));
00096   memset(Raft2Tpc, 0, 14*sizeof(TGeoHMatrix *));
00097   memset(Bundles2Tpc, 0,  14*6*sizeof(TGeoHMatrix *));
00098   memset(Mirrors2Tpc, 0,  14*6*7*sizeof(TGeoHMatrix *));
00099   for (Int_t i = 0; i < NoRaftPositions; i++) {
00100     if (! RaftPositions[i].Sector) continue;
00101     Int_t raft = RaftPositions[i].Raft;
00102     Traft[raft-1] = new TGeoHMatrix(Form("Raft%0i",raft));
00103     Traft[raft-1]->RotateX(RaftPositions[i].alpha*180/TMath::Pi());
00104     Traft[raft-1]->RotateY(RaftPositions[i].beta*180/TMath::Pi());
00105     Traft[raft-1]->RotateZ(RaftPositions[i].gamma*180/TMath::Pi());
00106     Traft[raft-1]->SetTranslation(&RaftPositions[i].X);
00107     if (Debug()) {
00108       RaftPositions[i].Print();
00109       Traft[raft-1]->Print();
00110     }
00111   }
00112   for (Int_t r = 0; r < 14; r++) {
00113     Int_t sector = Bundles[r][0].Sector;
00114     if (! sector) continue;
00115     Int_t raft = Bundles[r][0].Raft;
00116     Int_t half = 0;
00117     if (sector > 12) half = 1;
00118     Raft2Tpc[raft-1] = new TGeoHMatrix(Form("Raft%iToTpc",raft));
00119     *Raft2Tpc[raft-1] = RotSec[sector-1] * TpcHalf[half] * (*Traft[raft-1]);
00120     for (Int_t b = 0; b < 6; b++) {
00121       Int_t bundle = Bundles[r][b].Bundle;
00122       TGeoHMatrix rotB;
00123       rotB.SetTranslation(&Bundles[r][b].X);
00124       Bundles2Tpc[raft-1][bundle-1] = new TGeoHMatrix(Form("BundleR%iB%i",raft,bundle));
00125       *Bundles2Tpc[raft-1][bundle-1] = *Raft2Tpc[raft-1] * rotB;
00126       for (Int_t m = 0; m < 7; m++) {
00127         if (Mirrors[r][b][m].Sector < 2) continue;
00128         LASERINO_t *local = &Mirrors[r][b][m];
00129         if (! local) continue;
00130         Int_t mirror = Mirrors[r][b][m].Mirror;
00131         LaserBeams[NoBeams] = new LaserRaft();
00132         LaserBeams[NoBeams]->Sector = local->Sector;
00133         LaserBeams[NoBeams]->Raft   = local->Raft;
00134         LaserBeams[NoBeams]->Bundle = local->Bundle;
00135         LaserBeams[NoBeams]->Mirror = local->Mirror;
00136 #ifdef CORRECT_LASER_POSITIONS
00137        LaserBeams[NoBeams]->XyzB = StThreeVectorD(Mirrors[r][b][m].X+Mirrors[r][b][m].dX,
00138                                                   Mirrors[r][b][m].Y+Mirrors[r][b][m].dY,
00139                                                   Mirrors[r][b][m].Z+Mirrors[r][b][m].dZ);
00140 #else
00141        LaserBeams[NoBeams]->XyzB = StThreeVectorD(Mirrors[r][b][m].X,
00142                                                   Mirrors[r][b][m].Y,
00143                                                   Mirrors[r][b][m].Z);
00144 #endif
00145         Double_t theta = Bundles[r][b].ThetaZ + Mirrors[r][b][m].ThetaZ;
00146         Double_t phi   = Bundles[r][b].Phi    + Mirrors[r][b][m].Phi;
00147         TGeoHMatrix rotM;
00148         rotM.SetTranslation(LaserBeams[NoBeams]->XyzB.xyz());
00149         Mirrors2Tpc[raft-1][bundle-1][mirror-1] = new TGeoHMatrix(Form("MirrorR%iB%iM%i",raft,bundle,mirror));
00150         *Mirrors2Tpc[raft-1][bundle-1][mirror-1] = *Bundles2Tpc[raft-1][bundle-1] * rotM;
00151         
00152         LaserBeams[NoBeams]->dirB = StThreeVectorD(-TMath::Cos(phi)*TMath::Cos(theta), 
00153                                                    -TMath::Sin(phi)*TMath::Cos(theta), 
00154                                                    -TMath::Sin(theta));
00155         rotB.LocalToMaster(LaserBeams[NoBeams]->XyzB.xyz(), LaserBeams[NoBeams]->XyzU.xyz());
00156         rotB.LocalToMasterVect(LaserBeams[NoBeams]->dirB.xyz(), LaserBeams[NoBeams]->dirU.xyz());
00157         Raft2Tpc[raft-1]->LocalToMaster(LaserBeams[NoBeams]->XyzU.xyz(),LaserBeams[NoBeams]->XyzL.xyz());
00158         Raft2Tpc[raft-1]->LocalToMasterVect(LaserBeams[NoBeams]->dirU.xyz(),LaserBeams[NoBeams]->dirL.xyz());
00159         LaserBeams[NoBeams]->Theta = LaserBeams[NoBeams]->dirL.theta();
00160         LaserBeams[NoBeams]->Phi   = LaserBeams[NoBeams]->dirL.phi();
00161         if (Debug()) {
00162           local->Print();
00163           Raft2Tpc[raft-1]->Print();
00164           LaserBeams[NoBeams]->Print();
00165         }
00166         NoBeams++;
00167       }
00168     }
00169   }
00170   return StMaker::Init();
00171 }
00172 //_____________________________________________________________________________
00173 Int_t StLaserAnalysisMaker::InitRun(Int_t run){
00174   // average Z for membrane = -3.6 cm
00175   memset(Lasers, 0, NS*NB*NM*sizeof(LaserB *));
00176   const TGeoHMatrix &Tpc2Global = gStTpcDb->Tpc2GlobalMatrix();
00177   LaserB *theLaser = 0;
00178   for (Int_t i = 0; i < NoBeams; i++) {
00179     if (! LaserBeams[i]) continue;
00180     Int_t s = LaserBeams[i]->Sector/2 - 1;
00181     if (s < 0 || s >= NS) continue;
00182     Int_t b = LaserBeams[i]->Bundle - 1;
00183     if (b < 0 || b >= NB) continue;
00184     Int_t m = LaserBeams[i]->Mirror - 1;
00185     if (m < 0 || m >= NM) continue;
00186     theLaser = new LaserB(*LaserBeams[i]);
00187     Lasers[s][b][m] = theLaser;
00188     Tpc2Global.LocalToMaster(theLaser->XyzL.xyz(),theLaser->XyzG.xyz());
00189     Tpc2Global.LocalToMasterVect(theLaser->dirL.xyz(),theLaser->dirG.xyz());
00190     theLaser->PhiG = theLaser->dirG.phi();
00191     theLaser->ThetaG = theLaser->dirG.theta();
00192     theLaser->IsValid = 0;
00193 #if 0
00194     // From Laser Z
00195     if (theLaser->Sector ==  2 && theLaser->Bundle == 3 &&  theLaser->Mirror == 5) continue; // misplaced
00196     if (theLaser->Sector ==  4 && theLaser->Bundle == 4) continue; //dead
00197     if (theLaser->Sector ==  8) continue; // dead
00198     if (theLaser->Sector == 10) continue; // dead
00199     if (theLaser->Sector == 12 && theLaser->Bundle == 4 &&  theLaser->Mirror == 4) continue; // dead
00200     if (theLaser->Sector == 12 && theLaser->Bundle == 4 &&  theLaser->Mirror == 5) continue; // dead
00201     if (theLaser->Sector == 14 && theLaser->Bundle == 4 &&  theLaser->Mirror == 4) continue; // misplaced swap ?
00202     if (theLaser->Sector == 16 && theLaser->Bundle == 4 &&  theLaser->Mirror == 5) continue; // misplaced
00203     if (theLaser->Sector == 16 && theLaser->Bundle == 4 &&  theLaser->Mirror == 5) continue; // misplaced
00204     if (theLaser->Sector == 18 && theLaser->Bundle == 2) continue; // very strange pattern
00205     if (theLaser->Sector == 22 && theLaser->Bundle == 3) continue; // strange pattern
00206     if (theLaser->Sector == 22 && theLaser->Bundle == 4) continue; // missing
00207     if (theLaser->Sector == 24 && theLaser->Bundle == 2 && theLaser->Mirror == 6) continue; // missing
00208 #endif
00209     theLaser->IsValid = 1;
00210     if (Debug()) {
00211       LaserBeams[i]->Print();
00212       theLaser->Print();
00213     }
00214   }
00215   return kStOk;
00216 }
00217 //_____________________________________________________________________________
00218 void StLaserAnalysisMaker::Clear(const Option_t *option){
00219   if (event) event->Clear("C");
00220   StMaker::Clear(option);
00221 }
00222 //_____________________________________________________________________________
00223 Int_t StLaserAnalysisMaker::Make(){
00224   StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
00225     // Fill the event header.
00226   
00227   StEvtHddr *EvtHddr = (StEvtHddr*)GetDataSet("EvtHddr");
00228   if (! EvtHddr) return kStWarn;
00229   event->SetHeader(EvtHddr->GetEventNumber(), EvtHddr->GetRunNumber(), EvtHddr->GetDate(), EvtHddr->GetTime(),
00230                    gStTpcDb->Electronics()->tZero(), gStTpcDb->DriftVelocity(), gStTpcDb->Electronics()->samplingFrequency(), 
00231                    EvtHddr->GetInputTriggerMask());
00232   event->SetDVWest(gStTpcDb->DriftVelocity(1));
00233   event->SetDVEast(gStTpcDb->DriftVelocity(13));
00234   event->GetHeader()->SetDriftDistance(gStTpcDb->Dimensions()->gatingGridZ());
00235   event->GetHeader()->SetInnerSectorzOffset(gStTpcDb->Dimensions()->zInnerOffset());
00236   event->GetHeader()->SetOuterSectorzOffset(gStTpcDb->Dimensions()->zOuterOffset());
00237   event->GetHeader()->SettriggerTimeOffset(gStTpcDb->triggerTimeOffset());
00238   StDetectorDbClock* dbclock = StDetectorDbClock::instance();
00239   double freq = dbclock->getCurrentFrequency()/1000000.0;
00240   event->GetHeader()->SetOnlClock(freq);
00241   if (! pEvent) return kStWarn;
00242   UInt_t nvtx = pEvent->numberOfPrimaryVertices();
00243   for (UInt_t i = 0; i < nvtx; i++) {
00244     event->AddVertex(pEvent->primaryVertex(i));
00245   }
00246   StSPtrVecTrackNode& trackNode = pEvent->trackNodes();
00247   UInt_t nTracks = trackNode.size();
00248   StTrackNode *node=0;
00249   for (unsigned int i=0; i < nTracks; i++) {
00250     node = trackNode[i]; 
00251     if (!node) continue;
00252     StGlobalTrack  *gTrack = static_cast<StGlobalTrack *>(node->track(global));
00253     if (! gTrack) continue;
00254     //    if (gTrack->numberOfPossiblePoints(kTpcId) < 25) continue;
00255     StPrimaryTrack *pTrack =    static_cast<StPrimaryTrack*>(node->track(primary));
00256     StTrackFitTraits&  fitTraits =  gTrack->fitTraits();
00257     //    if (fitTraits.numberOfFitPoints(kTpcId) < 25) continue;
00258     StThreeVectorD g3 = gTrack->outerGeometry()->momentum();
00259     if (g3.mag() < 10) continue;
00260     StThreeVectorD unit = g3.unit();
00261     StPhysicalHelixD helixO = gTrack->outerGeometry()->helix();
00262     // find a match using last hit
00263     StPtrVecHit hvec = gTrack->detectorInfo()->hits();
00264     Double_t rMax = 0;
00265     Int_t jmax = -1;
00266     //    StThreeVectorD *pred = 0;
00267     StTpcHit *tpcHit;
00268     Int_t sector, s = -1;
00269     Int_t bundle = -1;
00270     Double_t dZ, dZmin  = 9999;
00271     Int_t b, m;
00272     Double_t dPhi, dPhimin  = 999;
00273     Int_t mmax = -1;
00274     StThreeVectorD Pred, Diff;
00275     Track *t = 0;
00276     for (unsigned int j=0; j<hvec.size(); j++) {// hit loop
00277       if (hvec[j]->detector() != kTpcId) continue;
00278       tpcHit = static_cast<StTpcHit *> (hvec[j]);
00279 #ifndef __TRACKHITS__
00280 #ifdef ADDPRIMTRACKHITS
00281       if (pTrack) 
00282         event->AddHit(tpcHit);
00283 #endif
00284 #else
00285       event->AddHit(tpcHit);
00286 #endif
00287       if (tpcHit->position().perp() > rMax) {
00288         rMax = tpcHit->position().perp();
00289         jmax = j;
00290       }
00291     }
00292     LaserB *theLaser = 0;
00293     if (jmax < 0 || rMax < 120) goto ADDTRACK;
00294     tpcHit = static_cast<StTpcHit *> (hvec[jmax]);
00295     // sector
00296     sector = tpcHit->sector();
00297     if (pTrack) goto ADDTRACK;
00298     if (2*(sector/2) != sector) goto ADDTRACK;
00299     s = sector/2 - 1;
00300     if (s < 0 || s >= NS) goto ADDTRACK;
00301     // bundle
00302     dZmin  = 9999;
00303     for (b = 0; b < NB; b++) {
00304       if (! Lasers[s][b][0]) continue;
00305       dZ = TMath::Abs(tpcHit->position().z() - Lasers[s][b][0]->XyzG.z());
00306       if (dZ < dZmin) {bundle = b; dZmin = dZ;}
00307     }
00308     if (bundle < 0 || dZmin > 15) goto ADDTRACK;
00309     // mirror
00310     // minimum distance in XY plane from laser spot
00311     dPhimin  = 999;
00312     for (m = 0; m < NM; m++) {
00313       if (! Lasers[s][bundle][m]) continue;
00314       dPhi = TMath::Abs(Lasers[s][bundle][m]->PhiG - g3.phi());
00315       if (dPhi < dPhimin) {
00316         dPhimin = dPhi;
00317         mmax = m;
00318       }
00319     }
00320     if (mmax < 0 || dPhimin > 0.1) goto ADDTRACK;
00321     theLaser = Lasers[s][bundle][mmax];
00322   ADDTRACK:
00323     t = event->AddTrack(sector,gTrack,theLaser,tpcHit->position().z());
00324     if (theLaser) {
00325       Int_t raft = theLaser->Raft;
00326       Int_t bundle = theLaser->Bundle;
00327       Int_t mirror = theLaser->Mirror;
00328       if (raft > 0 && raft <= 14 &&
00329           bundle > 0 && bundle <= 6 &&
00330           mirror > 0 && mirror <= 7) {
00331         t->SetPredictions(Raft2Tpc[raft-1], Bundles2Tpc[raft-1][bundle-1], Mirrors2Tpc[raft-1][bundle-1][mirror-1]);
00332         if (t->Flag == 2 || theLaser->IsValid)  event->AddTrackFit(t);
00333       }
00334     }
00335     if (Debug()) {
00336       t->Print();
00337     }
00338   }  
00339   static const Double_t sigma = 0.0385; // run 8090018 Full Field.
00340   for (Int_t k = 0; k < 11; k++) {
00341     FitDV *fit = (FitDV *) (*event->Fit())[k];
00342     Int_t N = fit->N;
00343     if (N > 2) {
00344       for (;;) {
00345         TRVector Y;
00346         TRMatrix A(0,2);
00347         for (Int_t i = 0; i < N; i++) {
00348           if (! fit->Flag[i]) {
00349             Double_t dev = fit->Y[i]/sigma;
00350             Y.AddRow(&dev);
00351             Double_t a[2] = {1./sigma, fit->X[i]/sigma};
00352             A.AddRow(a);
00353           }
00354         }
00355         Int_t ndf = A.GetNrows() - 2;
00356         if (ndf < 1) break;
00357         TRSymMatrix S(A,TRArray::kATxA);        if (Debug()) cout << "S: " << S << endl;
00358         TRVector    B(A,TRArray::kATxB,Y);      if (Debug()) cout << "B: " << B << endl;
00359         TRSymMatrix SInv(S,TRArray::kInverted); if (Debug()) cout << "SInv: " << SInv << endl;
00360         TRVector    X(SInv,TRArray::kSxA,B);    if (Debug()) cout << "X: " << X << endl;
00361         TRVector    R(Y);               
00362         R -= TRVector(A,TRArray::kAxB,X);       if (Debug()) cout << "R: " << R << endl;
00363         Double_t chisq = R*R;
00364         Double_t prob = TMath::Prob(chisq,ndf);
00365         if (prob > 0.01) {
00366           fit->offset  = X[0];
00367           fit->slope   = X[1];
00368           fit->chisq   = chisq;
00369           fit->dslope  = SInv[2];
00370           fit->doffset = SInv[0];
00371           fit->Prob = prob;
00372           fit->ndf = ndf;
00373           if (Debug()) fit->Print();
00374 #if 0
00375           TClonesArray &tracks = *event->Tracks();
00376           Int_t Ntrack = event->GetNtrack();
00377           for (Int_t i = 0; i < Ntrack; i++) {
00378             Track *t = (Track *) tracks[i];
00379             if (t->Flag == 2 && t->Laser.Sector == 2*(k+1)) {
00380               t->fit = *fit;
00381             }
00382           }
00383 #endif
00384           break;
00385         } else {
00386           Int_t j = -1;
00387           Int_t imax = -1;
00388           Double_t Rmax = 0;
00389           for (Int_t i = 0; i < N; i++) {
00390             if (! fit->Flag[i]) {
00391               j++;
00392               Double_t RR = TMath::Abs(R[j]);
00393               if (RR > Rmax) {
00394                 imax = i;
00395                 Rmax = RR;
00396               }
00397             }
00398           }
00399           if (imax < 0) break;
00400           fit->Flag[imax] = 1;
00401         }
00402       }
00403     }
00404   }
00405   m_laser->Fill(); //Fill the Tree
00406   return kStOK;
00407 }
00408 //________________________________________________________________________________
00409 Int_t StLaserAnalysisMaker::Finish() {
00410   for (Int_t s = 0; s < NS; s++) 
00411     for (Int_t b = 0; b < NB; b++) 
00412       for (Int_t m = 0; m < NM; m++) SafeDelete(Lasers[s][b][m]);
00413 #if 0
00414   if (! GetTFile()) {
00415     TSeqCollection *files = gROOT->GetListOfFiles();
00416     if (files) {
00417       TFile *f = 0;  
00418       TIter  next(files);
00419       while ((f = (TFile *) next())) {
00420         TString name(gSystem->BaseName(f->GetName()));
00421         if (name == "StLaserAnalysis.root") {
00422           delete f;
00423           break;
00424         }
00425       }
00426     }
00427   }
00428 #endif
00429   return StMaker::Finish();
00430 }