//----------------------------------------------------------- // File and Version Information: // $Id$ // // Description: // TPC Helitron/Plawa Matching Routine // // // Environment: // Software developed for the Prototype Detector at FOPI // // Author List: // Paul Buehler SMI // Philipp Muellner SMI // //----------------------------------------------------------- // This Class' Header ------------------ #include "TpcHelMatchingTaskPhilipp.h" // C++ headers // Collaborating Class Headers -------- #include "FairRun.h" #include "FairRunAna.h" #include "TpcDigiPar.h" #include "TpcCluster.h" #include "FairRuntimeDb.h" #include "TVectorD.h" #include "TH2.h" #include "TCanvas.h" #include "TGraph.h" #include "TMultiGraph.h" #include "TLatex.h" #include "TLegend.h" #include "RKTrackRep.h" #include "GFTrack.h" #include "GFTrackCand.h" #include "RecoHits/GFAbsRecoHit.h" #include "PseudoSpacePoint.h" #include "GFKalman.h" #include "FOPIField.h" #include "FopiEvent.h" #include "HeliHit.h" #include "HeliTrack.h" #include "PlawaTrack.h" #include "FopiForwardHit.h" #include "GFException.h" #define DEBUG 0 //----------------------------------------------------------- TpcHelMatchingTaskPhilipp::TpcHelMatchingTaskPhilipp() : fFopiEvBranchName("FopiEvent"), fTpcTrackFitBranchName("TpcTrackPostFit"), fHelHitAllBranchName("HeliHitAll"), fHelHitBranchName("HeliHit"), fHelTrackBranchName("HeliTrack"), fPlaTrackBranchName("PlawaTrack"), fTPCHelMatchParsBranchName("TpcHelMatchPars"), fTPCHelTrackBranchName("TpcHelPreFit"), fPersistence(kTRUE), fHHitScale(kTRUE) { fRad = 3.1415926 / 180.; fPosErr = TVector3(1., 1., 0.1); fMomErr = TVector3(0.1, 0.1, 0.1); fResScaleFac = 1; fHHScale = 1; fHResX = -1.; fHResY = -1.; fHResZ = 1; fminPmatch = 0; fPResX = 1.0; fPResY = 1.0; fPResZ = 1.0; } //----------------------------------------------------------- TpcHelMatchingTaskPhilipp::~TpcHelMatchingTaskPhilipp() { delete fPar; } //----------------------------------------------------------- InitStatus TpcHelMatchingTaskPhilipp::Init() { //Get ROOT Manager ------------------ FairRootManager *ioman = FairRootManager::Instance(); if (ioman == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "RootManager not instantiated!"); return kERROR; } // Get input ------------------------ fFopiEventArray = (TClonesArray *) ioman->GetObject(fFopiEvBranchName); if (fFopiEventArray == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "FOPI event not found!"); return kERROR; } fTpcTrackFitArray = (TClonesArray *) ioman->GetObject(fTpcTrackFitBranchName); if (fTpcTrackFitArray == 0) { fprintf(stderr, "GFTrack array not found!\n"); Error("TpcHelMatchingTaskPhilipp::Init", "GFTrack array not found!"); return kERROR; } fHelHitAllArray = (TClonesArray *) ioman->GetObject(fHelHitAllBranchName); if (fHelHitAllArray == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "Helitron Hits All array not found!"); return kERROR; } fHelHitArray = (TClonesArray *) ioman->GetObject(fHelHitBranchName); if (fHelHitArray == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "Helitron Hits array not found!"); return kERROR; } fHelTrackArray = (TClonesArray *) ioman->GetObject(fHelTrackBranchName); if (fHelTrackArray == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "Helitron track array not found!"); return kERROR; } fPlaTrackArray = (TClonesArray *) ioman->GetObject(fPlaTrackBranchName); if (fPlaTrackArray == 0) { Error("TpcHelMatchingTaskPhilipp::Init", "Plawa track array not found!"); return kERROR; } //register out branch // prepare TPC/Helitron/PLAWA Matching tree branch for output MatchingPars = new Double_t[23]; MatchingParsVec = new TVectorD(23); fMatchingPars = new TClonesArray("TVectorD"); ioman->Register(fTPCHelMatchParsBranchName, "Tpc", fMatchingPars, fPersistence); fTPCHelTracks = new TClonesArray("GFTrack"); ioman->Register(fTPCHelTrackBranchName, "Tpc", fTPCHelTracks, fPersistence); /* // prepare Helitron-hit-angles tree branch for output HelHitAngles = new Double_t[1]; HelHitAnglesVec = new TVectorD(1); fHelHitAngles = new TClonesArray("TVectorD"); ioman->Register("HelHitAngles","HelHitAngles",fHelHitAngles,fPersistence); */ return kSUCCESS; } //----------------------------------------------------------- void TpcHelMatchingTaskPhilipp::SetParContainers() { // std::cout<<"TpcResidualTask::SetParContainers"<GetRuntimeDb(); if (!db) Fatal("SetParContainers", "No runtime database"); // Get Tpc digitisation parameter container // fPar = (TpcDigiPar *) db->getContainer("TpcDigiPar"); // if (!fPar) // Fatal("SetParContainers", "TpcDigiPar not found"); } //----------------------------------------------------------- void TpcHelMatchingTaskPhilipp::Exec(Option_t * opt) { fprintf(stderr, "TpcHelMatchingTaskPhilipp:Exec\n"); // Definition of Helitron and Plawa hit covariance matrices TMatrixD HHcov = TMatrixD(3, 3); Double_t HHcovvals[9] = { 0., 0., 0., 0., 0., 0., 0., 0., pow(fHResZ, 2) }; TMatrixD PHcov = TMatrixD(3, 3); Double_t PHcovvals[9] = { pow(fPResX, 2), 0., 0., 0., pow(fPResY, 2), 0., 0., 0., pow(fPResZ, 2) }; PHcov.SetMatrixArray(PHcovvals, "F"); // clean up result buffer if (fTPCHelTracks == 0) Fatal("TpcHelMatchingTaskPhilipp::Exec()", "No fTPCHelTracks"); fTPCHelTracks->Delete(); if (fMatchingPars == 0) Fatal("TpcHelMatchingTaskPhilipp::Exec()", "No fMatchingPars"); fMatchingPars->Delete(); /* if(fHelHitAngles==0) Fatal("TpcHelMatchingTaskPhilipp::Exec()","No fHelHitAngles"); fHelHitAngles->Delete(); */ // FairRunAna *run = FairRunAna::Instance(); // FOPIField *fMagField = (FOPIField *) run->GetField(); // Double_t Bz = fMagField->GetBz(0., 0., 0.); // fprintf(stderr,"Magfield: %f\n",Bz); // get FOPI event FopiEvent *fopiEv = (FopiEvent *) fFopiEventArray->At(0); if (fopiEv == 0) return; Int_t nrun = fopiEv->GetRunNum(); Int_t nspill = fopiEv->GetSpillNum(); Int_t nevent = fopiEv->GetEvtNum(); /* printf("run: %i\n", nrun); printf("spill: %i\n", nspill); printf("event: %i\n", nevent); */ // transformation Helitron-to-beam Float_t taroff = fopiEv->GetTaroff(); Float_t phi0 = fopiEv->GetPhi0(); TVector3 hoff = fopiEv->GetHoff(); TVector3 hrot = fopiEv->GetHrot(); // define globhit and lochit FopiForwardHit globhit = FopiForwardHit(1, 0, taroff, hoff, hrot); FopiForwardHit lochit = FopiForwardHit(0, 0, taroff, hoff, hrot); // how many tracks are there Int_t NumTpcTracks = fTpcTrackFitArray->GetEntriesFast(); Int_t NumHelTracks = fopiEv->GetHmul(); ////////////////////////////////////////// // TPC // ////////////////////////////////////////// // prepare Int_t *NumTpcTrackHits = new Int_t[NumTpcTracks]; int tcounter = 0; TVector3 TPCmom, pos, deltapos; Double_t *theTPC = new Double_t[3 * NumTpcTracks]; //TPC track Theta Double_t *phiTPC = new Double_t[3 * NumTpcTracks]; //TPC track Phi Double_t *momTPC = new Double_t[NumTpcTracks]; //TPC track Momentum Double_t *chaTPC = new Double_t[NumTpcTracks]; //TPC track Charge //analyse TPC tracks for (Int_t ttrack = 0; ttrack < NumTpcTracks; ttrack++) { GFTrack *tpcsingletrack = (GFTrack *) fTpcTrackFitArray->At(ttrack); NumTpcTrackHits[ttrack] = tpcsingletrack->getNumHits(); momTPC[ttrack] = tpcsingletrack->getMom().Mag(); chaTPC[ttrack] = tpcsingletrack->getCharge(); TPCmom = tpcsingletrack->getMom(); //angles from the track fit theTPC[3 * ttrack] = TPCmom.Theta() / fRad; phiTPC[3 * ttrack] = TPCmom.Phi() / fRad; if (phiTPC[3 * ttrack] < 0) phiTPC[3 * ttrack] += 360; // estimate theta and phi using hit points Float_t xmean = 0; Float_t ymean = 0; theTPC[3 * ttrack + 1] = 0; phiTPC[3 * ttrack + 1] = 0; if (NumTpcTrackHits[ttrack] > 1) { for (Int_t ii = 0; ii < NumTpcTrackHits[ttrack]; ii++) { GFAbsRecoHit *tsinglehit = tpcsingletrack->getHit(ii); TVectorD raw(3); raw = tsinglehit->getRawHitCoord(); pos.SetXYZ(raw[0],raw[1],raw[2]); Double_t posPerp = pos.Perp(); theTPC[3 * ttrack + 1] += pos.Theta(); xmean += pos.X() / posPerp; ymean += pos.Y() / posPerp; } theTPC[3 * ttrack + 1] = theTPC[3 * ttrack + 1] / NumTpcTrackHits[ttrack] / fRad; phiTPC[3 * ttrack + 1] = atan2(ymean, xmean) / fRad; if (phiTPC[3 * ttrack + 1] < 0) phiTPC[3 * ttrack + 1] += 360; } else { theTPC[3 * ttrack + 1] = -360; phiTPC[3 * ttrack + 1] = -360; } //estimate theta and phi using hit pairs theTPC[3 * ttrack + 2] = 0; phiTPC[3 * ttrack + 2] = 0; xmean = 0; ymean = 0; tcounter = 0; if (NumTpcTrackHits[ttrack] > 1) { for (Int_t ii = 0; ii < NumTpcTrackHits[ttrack]; ii++) { GFAbsRecoHit *tsinglehit = tpcsingletrack->getHit(ii); TVectorD raw(3); raw = tsinglehit->getRawHitCoord(); pos.SetXYZ(raw[0],raw[1],raw[2]); for (int jj = ii + 1; jj < NumTpcTrackHits[ttrack]; jj++) { GFAbsRecoHit *tsinglehit2 = tpcsingletrack->getHit(jj); TVectorD raw2(3); raw2 = tsinglehit2->getRawHitCoord(); deltapos.SetXYZ(raw2[0],raw2[1],raw2[2]); deltapos = deltapos - pos; Double_t deltaposPerp = deltapos.Perp(); // theta/phi calculation theTPC[3 * ttrack + 2] += deltapos.Theta(); xmean += deltapos.X() / deltaposPerp; ymean += deltapos.Y() / deltaposPerp; tcounter++; } } theTPC[3 * ttrack + 2] = theTPC[3 * ttrack + 2] / tcounter / fRad; phiTPC[3 * ttrack + 2] = atan2(ymean, xmean) / fRad; if (phiTPC[3 * ttrack + 2] < 0) phiTPC[3 * ttrack + 2] += 360; } else { theTPC[3 * ttrack + 2] = -360; phiTPC[3 * ttrack + 2] = -360; } } ////////////////////////////////////////// // Helitron // ////////////////////////////////////////// //prepare Double_t hrsec; Int_t hhrc = -1; Int_t *NumHelTrackHits = new Int_t[NumHelTracks]; TVector3 Hpos, Hposall; TVector3 *Ppos = new TVector3[NumHelTracks]; Double_t *theHel = new Double_t[2 * NumHelTracks]; //Hel track Theta Double_t *phiHel = new Double_t[2 * NumHelTracks]; //Hel track Phi Double_t *momHel = new Double_t[NumHelTracks]; //Hel track momentum Double_t *charHel = new Double_t[NumHelTracks]; //Hel track charge Double_t *massHel = new Double_t[NumHelTracks]; //Hel track mass Double_t *dedxHel = new Double_t[NumHelTracks]; //Hel track dEdx Double_t *phmatch = new Double_t[NumHelTracks]; //Quality of matched PLAWA hit Double_t *PHvel = new Double_t[NumHelTracks]; //PLAWA velocity (cm/ns) // analyse Helitron tracks Int_t counter = -1; for (Int_t htrack = 0; htrack < NumHelTracks; htrack++) { HeliTrack *hsingletrack = (HeliTrack *) fHelTrackArray->At(htrack); PlawaTrack *psingletrack = (PlawaTrack *) fPlaTrackArray->At(htrack); NumHelTrackHits[htrack] = hsingletrack->GetHhmul(); // momentum and Helitron/PLAWA matching theHel[2 * htrack] = hsingletrack->GetHthe(); phiHel[2 * htrack] = hsingletrack->GetHphi(); momHel[htrack] = hsingletrack->GetHmom(); massHel[htrack] = hsingletrack->GetHmass(); charHel[htrack] = hsingletrack->GetHchar(); dedxHel[htrack] = hsingletrack->GetHdedx(); phmatch[htrack] = hsingletrack->GetPHmatch(); PHvel[htrack] = psingletrack->GetPhvel(); // reference sector hrsec = hsingletrack->GetHrsec(); lochit.SetHrsec(hrsec); // compute theta and phi from track hits Float_t xmean = 0; Float_t ymean = 0; theHel[2 * htrack + 1] = 0; phiHel[2 * htrack + 1] = 0; for (Int_t ii = 0; ii < NumHelTrackHits[htrack]; ii++) { hhrc++; HeliHit *hsinglehit = (HeliHit *) fHelHitArray->At(hhrc); // transform into global coordinate system lochit.SetPos(0, hsinglehit->GetHhx(), hsinglehit->GetHhy(), hsinglehit->GetHhz()); //get position of Helitron track hitpoints Hpos = TVector3(lochit.GetGlobalPos().X(), lochit.GetGlobalPos().Y(), lochit.GetGlobalPos().Z()); // Hpos.Print(); Double_t HposPerp = Hpos.Perp(); //angles theHel[2 * htrack + 1] += Hpos.Theta(); xmean += Hpos.X() / HposPerp; ymean += Hpos.Y() / HposPerp; } theHel[2 * htrack + 1] = theHel[2 * htrack + 1] / NumHelTrackHits[htrack] / fRad; phiHel[2 * htrack + 1] = atan2(ymean, xmean) / fRad; if (phiHel[2 * htrack + 1] < 0) phiHel[2 * htrack + 1] += 360; // get matching PLAWA hit if (phmatch[htrack] > fminPmatch) { Double_t phthe = fRad * psingletrack->GetPhthe(); Double_t phphi = fRad * psingletrack->GetPhphi(); Double_t pz = psingletrack->GetPhz(); Double_t pr = pz * tan(phthe); Ppos[htrack] = TVector3(pr * cos(phphi), pr * sin(phphi), pz); } } //get phi-angle of all Heli-Hits /* Int_t NumHelHitAll = fHelHitAllArray->GetEntriesFast(); for(Int_t jj=0; jjAt(jj); hrsec = allhelihit->GetHrsec(); lochit.SetHrsec(hrsec); // transform into global coordinate system lochit.SetPos(0,allhelihit->GetHhx(), allhelihit->GetHhy(), allhelihit->GetHhz()); Hposall=TVector3(lochit.GetGlobalPos().X(), lochit.GetGlobalPos().Y(), lochit.GetGlobalPos().Z()); //angles AHHA = Hposall.Phi()/fRad; if(AHHA<0)AHHA+=360; HelHitAngles[0] = AHHA; counter++; HelHitAnglesVec->SetElements(HelHitAngles); new((*fHelHitAngles)[counter]) TVectorD(*HelHitAnglesVec); } */ ////////////////////////////////////////////////////////// // Combining of TPC and Helitron Hits // ////////////////////////////////////////////////////////// Int_t hitcounter = -1; counter = -1; Int_t numcomb = NumTpcTracks*NumHelTracks; if (numcomb>0) { for (Int_t ttrack = 0; ttrack < NumTpcTracks; ttrack++) { GFTrack *tpcsingletrack = (GFTrack *) fTpcTrackFitArray->At(ttrack); hhrc = -1; for (Int_t htrack = 0; htrack < NumHelTracks; htrack++) { counter++; //definition of representation TVector3 pos0, mom0, poserr, momerr; Int_t pdg1 = 211; //pion Int_t pdg2 = 2212; //proton GFAbsRecoHit *firsthit = tpcsingletrack->getHit(0); TVectorD raw(3); raw = firsthit->getRawHitCoord(); pos0.SetXYZ(raw[0],raw[1],raw[2]); mom0 = TVector3(tpcsingletrack->getMom()); mom0.SetMagThetaPhi(1, fRad * theTPC[3 * ttrack + 2], fRad * phiTPC[3 * ttrack + 2]); poserr = fPosErr; momerr = TVector3(mom0.X() * fMomErr.X(), mom0.Y() * fMomErr.Y(), mom0.Z() * fMomErr.Z()); if (phiTPC[3 * ttrack + 2] - phiHel[2 * htrack] < 0) { pdg1 = -1 * pdg1; pdg2 = -1 * pdg1; } RKTrackRep* rep1 = new RKTrackRep(pos0, mom0, poserr, momerr, pdg1); RKTrackRep* rep2 = new RKTrackRep(pos0, mom0, poserr, momerr, pdg2); // fprintf(stderr,"rep[%i] created\n",counter); //define new track GFTrack* trk = new GFTrack(rep1); trk->addTrackRep(rep2); // fprintf(stderr,"trk[%i] cretated\n",counter); hitcounter = -1; HeliTrack *hsingletrack = (HeliTrack *) fHelTrackArray->At(htrack); PlawaTrack *psingletrack = (PlawaTrack *) fPlaTrackArray->At(htrack); hrsec = hsingletrack->GetHrsec(); lochit.SetHrsec(hrsec); //TPC clusters for (Int_t ii = 0; ii < NumTpcTrackHits[ttrack]; ii++) { hitcounter++; GFAbsRecoHit *TPCHit = tpcsingletrack->getHit(ii); trk->addHit(TPCHit); } // fprintf(stderr,"%i TPC hit points added\n",NumTpcTrackHits[ttrack]); //Set Cov-Matrix for Hel-Tracks Double_t HresX, HresY, HresZ; if (fHResX < 0.) { HresX = hsingletrack->GetHdx(); } else { HresX = fHResX; } if (fHResY < 0.) { HresY = hsingletrack->GetHdy(); } else { HresY = fHResY; } HresZ = fHResZ; // if (fHResZ<0.) {HresZ = hsingletrack->GetHdz();} // else {HresZ = fHResZ;} HHcovvals[0] = pow(fResScaleFac * fHHScale * HresX, 2); HHcovvals[4] = pow(fResScaleFac * fHHScale * HresY, 2); HHcovvals[8] = pow(fResScaleFac * HresZ, 2); HHcov.SetMatrixArray(HHcovvals, "F"); lochit.SetCov(0, HHcov); //Helitron hitpoints PseudoSpacePoint *fwdhits = NULL; if (NumHelTrackHits[htrack] > 0) { Int_t numfwdhits = NumHelTrackHits[htrack]; if (phmatch[htrack] > fminPmatch) numfwdhits++; fwdhits = new PseudoSpacePoint [numfwdhits]; for (Int_t jj = 0; jj < NumHelTrackHits[htrack]; jj++) { hhrc++; HeliHit *hsinglehit = (HeliHit *) fHelHitArray->At(hhrc); // transform into global coordinate system lochit.SetPos(0, hsinglehit->GetHhx(), hsinglehit->GetHhy(), hsinglehit->GetHhz()); fwdhits[jj] = PseudoSpacePoint(&lochit); trk->addHit(&fwdhits[jj]); } // fprintf(stderr,"%i Helitron hit points added\n",NumHelTrackHits[htrack]); // Matched PLAWA hitpoint // fprintf(stderr,"phmatch[%i]: %f\n", htrack, phmatch[htrack]); if (phmatch[htrack] > fminPmatch) { globhit.SetPos(1,Ppos[htrack].X(),Ppos[htrack].Y(),Ppos[htrack].Z()); globhit.SetCov(0,PHcov); fwdhits[NumHelTrackHits[htrack]] = PseudoSpacePoint(&globhit); trk->addHit(&fwdhits[NumHelTrackHits[htrack]]); // fprintf(stderr,"Plawa hit point added\n"); } } new((*fTPCHelTracks)[counter]) GFTrack(*trk); // fprintf(stderr,"TPCHelitrk[%i] saved\n",counter); // delete trk; trk = NULL; delete rep1; rep1 = NULL; delete rep2; rep2 = NULL; delete [] fwdhits; fwdhits = NULL; } } } //MatchingPars Branch filling counter = -1; if (NumTpcTracks > 0 && NumHelTracks > 0) { for (Int_t ttrack = 0; ttrack < NumTpcTracks; ttrack++) { MatchingPars[0] = nrun; //run number MatchingPars[1] = nspill; //spill number MatchingPars[2] = nevent; //event number MatchingPars[3] = ttrack; //TPC track number MatchingPars[4] = theTPC[3 * ttrack]; //Theta of TPC track fit MatchingPars[5] = theTPC[3 * ttrack + 1]; //Theta of TPC track clusters (with vertex) MatchingPars[6] = theTPC[3 * ttrack + 2]; //Theta of TPC track cluster pairs (w/o vertex) MatchingPars[7] = phiTPC[3 * ttrack]; //Phi of TPC track fit MatchingPars[8] = phiTPC[3 * ttrack + 1]; //Phi of TPC track clusters (with vertex) MatchingPars[9] = phiTPC[3 * ttrack + 2]; //Phi of TPC track cluster pairs (w/o vertex) MatchingPars[10] = momTPC[ttrack]; //Momentum of TPC tracks MatchingPars[11] = chaTPC[ttrack]; //Charge of TPC tracks for (Int_t htrack = 0; htrack < NumHelTracks; htrack++) { MatchingPars[12] = htrack; //Hel track number MatchingPars[13] = theHel[2 * htrack]; //Hel track Theta from fit MatchingPars[14] = theHel[2 * htrack + 1]; //Hel track Theta from track points MatchingPars[15] = phiHel[2 * htrack]; //Hel track Phi from fit MatchingPars[16] = phiHel[2 * htrack + 1]; //Hel track Phi from track points MatchingPars[17] = momHel[htrack]; //Momentum of Helitron tracks MatchingPars[18] = dedxHel[htrack]; //dEdx of Helitron tracks MatchingPars[19] = phmatch[htrack]; //Matched Plawa hit MatchingPars[20] = charHel[htrack]; //Char of Helitron tracks MatchingPars[21] = massHel[htrack]; //Mass of Helitron tracks MatchingPars[22] = PHvel[htrack]; //Velocity of PLAWA hit counter++; MatchingParsVec->SetElements(MatchingPars); new((*fMatchingPars)[counter]) TVectorD(*MatchingParsVec); } } } else if (NumTpcTracks > 0 && NumHelTracks == 0) { for (Int_t ttrack = 0; ttrack < NumTpcTracks; ttrack++) { MatchingPars[0] = nrun; //run number MatchingPars[1] = nspill; //spill number MatchingPars[2] = nevent; //event number MatchingPars[3] = ttrack; //TPC track number MatchingPars[4] = theTPC[3 * ttrack]; //Theta of TPC track fit MatchingPars[5] = theTPC[3 * ttrack + 1]; //Theta of TPC track clusters (with vertex) MatchingPars[6] = theTPC[3 * ttrack + 2]; //Theta of TPC track cluster pairs (w/o vertex) MatchingPars[7] = phiTPC[3 * ttrack]; //Phi of TPC track fit MatchingPars[8] = phiTPC[3 * ttrack + 1]; //Phi of TPC track clusters (with vertex) MatchingPars[9] = phiTPC[3 * ttrack + 2]; //Phi of TPC track cluster pairs (w/o vertex) MatchingPars[10] = momTPC[ttrack]; //Momentum of TPC tracks MatchingPars[11] = chaTPC[ttrack]; //Charge of TPC tracks MatchingPars[12] = -100; //Hel track number MatchingPars[13] = -100; //Hel track Theta from fit MatchingPars[14] = -100; //Hel track Theta from track points MatchingPars[15] = -100; //Hel track Phi from fit MatchingPars[16] = -100; //Hel track Phi from track points MatchingPars[17] = -100; //Momentum of Helitron tracks MatchingPars[18] = -100; //dEdx of Helitron tracks MatchingPars[19] = -100; //Matched Plawa hit MatchingPars[20] = -100; //Char of Helitron tracks MatchingPars[21] = -100; //Mass of Helitron tracks MatchingPars[22] = -100; //Velocity of PLAWA hit counter++; MatchingParsVec->SetElements(MatchingPars); new((*fMatchingPars)[counter]) TVectorD(*MatchingParsVec); } } else if (NumTpcTracks == 0 && NumHelTracks > 0) { MatchingPars[0] = nrun; //run number MatchingPars[1] = nspill; //spill number MatchingPars[2] = nevent; //event number MatchingPars[3] = -100; //TPC track number MatchingPars[4] = -100; //Theta of TPC track fit MatchingPars[5] = -100; //Theta of TPC track clusters (with vertex) MatchingPars[6] = -100; //Theta of TPC track cluster pairs (w/o vertex) MatchingPars[7] = -100; //Phi of TPC track fit MatchingPars[8] = -100; //Phi of TPC track clusters (with vertex) MatchingPars[9] = -100; //Phi of TPC track cluster pairs (w/o vertex) MatchingPars[10] = -100; //Momentum of TPC tracks MatchingPars[11] = -100; //Charge of TPC tracks for (Int_t htrack = 0; htrack < NumHelTracks; htrack++) { MatchingPars[12] = htrack; //Hel track number MatchingPars[13] = theHel[2 * htrack]; //Hel track Theta from fit MatchingPars[14] = theHel[2 * htrack + 1]; //Hel track Theta from track points MatchingPars[15] = phiHel[2 * htrack]; //Hel track Phi from fit MatchingPars[16] = phiHel[2 * htrack + 1]; //Hel track Phi from track points MatchingPars[17] = momHel[htrack]; //Momentum of Helitron tracks MatchingPars[18] = dedxHel[htrack]; //dEdx of Helitron tracks MatchingPars[19] = phmatch[htrack]; //Matched Plawa hit MatchingPars[20] = charHel[htrack]; //Char of Helitron tracks MatchingPars[21] = massHel[htrack]; //Mass of Helitron tracks MatchingPars[22] = PHvel[htrack]; //Velocity of PLAWA hit counter++; MatchingParsVec->SetElements(MatchingPars); new((*fMatchingPars)[counter]) TVectorD(*MatchingParsVec); } } delete [] NumTpcTrackHits; delete [] NumHelTrackHits; delete [] theTPC; delete [] phiTPC; delete [] momTPC; delete [] chaTPC; delete [] theHel; delete [] phiHel; delete [] momHel; delete [] charHel; delete [] massHel; delete [] dedxHel; delete [] phmatch; delete [] PHvel; delete [] Ppos; } ClassImp(TpcHelMatchingTaskPhilipp) //-----------------------------------------------------------