|
//____________________________________________________________________ // // RICH pid module - makes particle identification from RICH data. // // Find the Bfs track table and the RichRdo object in the inNode. // Calculate the track lines reflection on the RICH mirror and find // the intersection with the image plane - this should be the center // of a ring generated by the particle. // // Find a ring by looking at the number of hits in delta radius areas // around the center. // //____________________________________________________________________ // // $Id: BrRichPidModule.cxx,v 1.15 2002/06/10 13:30:57 ekman Exp $ // $Author: ekman $ // $Date: 2002/06/10 13:30:57 $ // $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrRichPidModule #include "BrRichPidModule.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrChkvRdo #include "BrChkvRdo.h" #endif #ifndef BRAT_BrRichPid #include "BrRichPid.h" #endif #ifndef BRAT_BrSpectrometerTracks #include "BrSpectrometerTracks.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef ROOT_TF1 #include "TF1.h" #endif #ifndef ROOT_TProfile #include "TProfile.h" #endif #ifndef ROOT_TVector3 #include "TVector3.h" #endif #ifndef WIN32 #include <iostream> #else #include <iostream.h> #endif #ifndef __IOMANIP__ #include <iomanip> #endif //____________________________________________________________________ ClassImp(BrRichPidModule); //____________________________________________________________________ BrRichPidModule::BrRichPidModule() { // Default constructor. DO NOT USE SetState(kSetup); } //____________________________________________________________________ BrRichPidModule::BrRichPidModule(const Char_t* name, const Char_t* title) : BrModule(name, title) { // Named Constructor SetState(kSetup); SetDRWindow(); SetMinHitsInRing(); SetCheckGeometry(); SetMaxRadius(); SetRefractiveIndex(); SetUseD3Momentum(); SetUseD4Momentum(); SetImagePlaneOffset(); fNSnapshots = 0; fRichVolume = 0; fT5Volume = 0; fRichParams = 0; fNDEDR = 0; } //____________________________________________________________________ void BrRichPidModule::DefineHistograms() { // Define histograms. if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } TDirectory* saveDir = gDirectory; fHistDir = gDirectory->mkdir(GetName()); fHistDir->cd(); // definition of functions Char_t* formula_e = Form("%f*tan(acos((1/(%f))*sqrt(1+(%f/(x*x)))))", fMirrorFocalLength,fRefractiveIndex,0.0005*0.0005); Char_t* formula_mu = Form("%f*tan(acos((1/(%f))*sqrt(1+(%f/(x*x)))))", fMirrorFocalLength,fRefractiveIndex,0.105*0.105); Char_t* formula_pi = Form("%f*tan(acos((1/(%f))*sqrt(1+(%f/(x*x)))))", fMirrorFocalLength,fRefractiveIndex,0.139*0.139); Char_t* formula_k = Form("%f*tan(acos((1/(%f))*sqrt(1+(%f/(x*x)))))", fMirrorFocalLength,fRefractiveIndex,0.494*0.494); Char_t* formula_p = Form("%f*tan(acos((1/(%f))*sqrt(1+(%f/(x*x)))))", fMirrorFocalLength,fRefractiveIndex,0.938*0.938); fElectronRvsP = new TF1("theoretical_e",formula_e,-30,30); fElectronRvsP->SetNpx(100000); fMuonRvsP = new TF1("theoretical_mu",formula_mu,-30,30); fPionRvsP = new TF1("theoretical_pi",formula_pi,-30,30); fKaonRvsP = new TF1("theoratical_k",formula_k,-30,30); fProtonRvsP = new TF1("theoretical_p",formula_p,-30,30); fHistDir->Add(fElectronRvsP); fHistDir->Add(fMuonRvsP); fHistDir->Add(fPionRvsP); fHistDir->Add(fKaonRvsP); fHistDir->Add(fProtonRvsP); // definitions of histograms fRingCenters = new TH2F("ringCenters", "Calculated center of rings on the image plane", 100,-13.666,13.666,100,-11.266,11.266); fDRVsPhi = new TH2F("dRVsPhi", "(rTube - rAverage) versus phi", 50,-3.2,3.2,100,-3,3); fDRVsPhiAcc = new TH2F("dRVsPhiAcc", "(rTube - rAverage) versus phi (accepted)", 50,-3.2,3.2,100,-3,3); fRing = new TH2F("ring","average ring",200,-2,2,200,-2,2); fRingAcc = new TH2F("ringAcc","average ring accepted hits",200,-2,2,200,-2,2); fMass = new TH1F("mass","mass spectrum", 750, -1.5, 1.5); fRadVsP = new TH2F("radiusVsMomentum","Radius versus momentum",250,-25,25,100,-1,15); fRadVsP1 = new TH2F("radiusVsMomentum_nTracks=1","Radius versus momentum, nTracks=1",650,-25,25,320,-1,15); fRadVsP2 = new TH2F("radiusVsMomentum_nTracks>1","Radius versus momentum, nTracks>1",650,-25,25,320,-1,15); fHitsInRing = new TH1F("hitsInRing","Number of hits per ring", 50,0.5,50.5); fHitsVsRadius= new TH2F("hitsVsRadius","Number of ring hits vs r", 100,-1,15,50,0.5,50.5 ); fHitsVsPhotons = new TH2F("hitsVsPhotons", "Hits versus number of photons", 51,-0.5,50.5,101,-0.5,100.5); fHitsVsEstPhotons = new TH2F("hitsVsEstPhotons", "Hits versus estimated number of photons", 51,-0.5,50.5,51,-0.5,50.5); fMassVsP = new TH2F("massVsMomentum","massVsMomentum",500,-25,25,750,-1.5,1.5); fMassVsRad = new TH2F("massVsRadius","massVsRadius",320,-1,15,750,-1.5,1.5); fSnapshots = new TList(); fHistDir->mkdir("Snapshots"); fHistDir->cd("Snapshots"); for (Int_t i = 0; i <= fMaxNSnapshots; i++) fSnapshots->Add(new TH2F(Form("Snapshot%d", i), Form("Snapshot of image plane, no.%d", i), 22,-13.666,13.666,18,-11.266,11.266)); fHistDir->cd(); fDEDR = new TList(); fHistDir->mkdir("DEDR"); // setting titles etc. fMass->SetLineWidth(3); fMass->SetXTitle("mass [GeV/c^{2}]"); fRadVsP->SetMarkerSize(0.3); fRadVsP->SetMarkerStyle(20); fRadVsP->SetXTitle("reconstructed track momentum [GeV/c]"); fRadVsP->SetYTitle("RICH ring radius [cm]"); fRadVsP1->SetMarkerSize(0.3); fRadVsP1->SetMarkerStyle(20); fRadVsP1->SetXTitle("reconstructed track momentum [GeV/c]"); fRadVsP1->SetYTitle("RICH ring radius [cm]"); fRadVsP2->SetMarkerSize(0.3); fRadVsP2->SetMarkerStyle(20); fRadVsP2->SetXTitle("reconstructed track momentum [GeV/c]"); fRadVsP2->SetYTitle("RICH ring radius [cm]"); fHitsVsRadius->SetXTitle("radius [cm]"); fHitsVsRadius->SetYTitle("hits in ring"); fHitsVsEstPhotons->SetXTitle("estimated number of photons"); fHitsVsEstPhotons->SetYTitle("hits in ring"); fHitsVsPhotons->SetXTitle("summed energy signal (photons)"); fHitsVsPhotons->SetYTitle("hits in ring"); gDirectory = saveDir; } //____________________________________________________________________ void BrRichPidModule::Init() { // In this method the following is initialized: // // 1) Geometry and Parameter db managers // 2) Rich parameters (tube numbering and positions) // 3) The Geometry of the RICH and T5 // 4) Rich mirror information (focal length, mirror sphere radius...) // 5) Image plane // Job-level initialisation SetState(kInit); fGausFit = new TF1("Fit", "gaus", 0, 15); BrGeometryDbManager *geoDb = BrGeometryDbManager::Instance(); BrParameterDbManager *parDb = BrParameterDbManager::Instance(); if (!geoDb) { Failure("Init ", "Couldn't instantiate geometry manager"); return; } if (!parDb) { Failure("Init ", "Couldn't instantiate detector parameter manager"); return; } fRichParams = (BrChkvParameters*)parDb->GetDetectorParameters("BrChkvParameters", "RICH"); fRichVolume = (BrDetectorVolume*)geoDb->GetDetectorVolume("BrDetectorVolume","RICH"); fT5Volume = (BrDetectorVolume*)geoDb->GetDetectorVolume("BrDetectorVolume","T5"); //------------------------------------------------------------------ // Initialization of mirror geometry fMirrorFocalAngle = 9. / (180/TMath::Pi()); // in radians fMirrorSphereRadius = 300; // in cm fMirrorFocalLength = 0.5 * fMirrorSphereRadius; // in cm Float_t zOff = 5; // fixed offset Double_t x = 0; Double_t y = fMirrorSphereRadius * TMath::Sin(fMirrorFocalAngle); Double_t z = ( fRichVolume->GetSize()[2]/2 - zOff - fMirrorSphereRadius * TMath::Cos(fMirrorFocalAngle) ); fMirrorSphereCenter.SetX(x); fMirrorSphereCenter.SetY(y); fMirrorSphereCenter.SetZ(z); // Definition of image plane // // the image plane is situated ~150 cm from the center of the // mirror, in an angle of 2 * the focal angle (seen from the // z-axis). TVector3 mirrorCenter(0,0,(fMirrorFocalLength * 2 * TMath::Cos(fMirrorFocalAngle) + fMirrorSphereCenter.Z())); TVector3 focalLineR(0, 0, -1); focalLineR.Rotate(2*fMirrorFocalAngle, TVector3(1,0,0)); BrVector3D focalLine(focalLineR.X(),focalLineR.Y(),focalLineR.Z()); Float_t halfRadius = fMirrorSphereRadius/2; // offsets are in local coordinates Float_t offsety = fImagePlaneOffset.Y(); fImagePlaneOffset.SetY(offsety * TMath::Sin((TMath::Pi()/2) - 2* fMirrorFocalAngle)); fImagePlaneOffset.SetZ(offsety * TMath::Cos((TMath::Pi()/2) - 2* fMirrorFocalAngle)); fImagePlaneCenter = BrVector3D(mirrorCenter.X() + halfRadius*focalLine.X() + fImagePlaneOffset.X(), mirrorCenter.Y() + halfRadius*focalLine.Y() + fImagePlaneOffset.Y(), mirrorCenter.Z() + halfRadius*focalLine.Z() + fImagePlaneOffset.Z()); // I need one point and a normal vector to define the plane fImagePlane = BrPlane3D(fImagePlaneCenter,focalLine); if (DebugLevel()>15) { cout << "BrRichPidModule::Init: Image plane middle: " << fImagePlaneCenter.X() << " " << fImagePlaneCenter.Y() << " " << fImagePlaneCenter.Z() << endl; } } //____________________________________________________________________ void BrRichPidModule::Begin() { // Run-level initialisation SetState(kBegin); } //____________________________________________________________________ void BrRichPidModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Following is done in the event method: // // Per event method SetState(kEvent); // get bfs track table BrDataTable* bfsTrackTable = (BrDataTable*) inNode->GetDataTable("BfsTracks"); if (!bfsTrackTable) { bfsTrackTable = (BrDataTable*) outNode->GetDataTable("BfsTracks"); if (!bfsTrackTable) { if (Verbose()>15) Warning("Event","Could not find Bfs track table."); return; } } // get rhic rdo BrChkvRdo* richRdo = (BrChkvRdo*) inNode->GetObject("ChkvRdo RICH"); if (!richRdo) { richRdo = (BrChkvRdo*) outNode->GetObject("ChkvRdo RICH"); if (!richRdo) { if (Verbose()>15) Warning("Event","Could not find ChkvRdo Rich."); return; } } // create data table with BrRichPid objects BrDataTable* richPidTable = new BrDataTable("RichPid","Rich pid object table"); outNode->AddDataTable(richPidTable); //------------------------------------------------------------------------- //Now get the track info and the rich hit (tube) info: // get the number of tracks fNumberOfTracks = bfsTrackTable->GetEntries(); if ( (fNumberOfTracks == 0) || (fNumberOfTracks>3) ) return; //highest momentum first... if (fNumberOfTracks>1) { for (Int_t u = 0; u<fNumberOfTracks; u++) { Float_t p = ((BrBfsTrack*)bfsTrackTable->At(u))->GetMomentum(); Int_t m = u; for (Int_t t = u+1; t<fNumberOfTracks; t++) { Float_t p1 = ((BrBfsTrack*)bfsTrackTable->At(t))->GetMomentum(); if (p1 >= p) { p = p1; m = t; } } if (m != u) { BrBfsTrack* t = (BrBfsTrack*)bfsTrackTable->At(u); bfsTrackTable->AddAt(bfsTrackTable->At(m), u); bfsTrackTable->AddAt(t, m); } } } //copy tracks into track array TRACKHIT trackHits[fNumberOfTracks]; for (Int_t t = 0; t<fNumberOfTracks; t++) { BrBfsTrack* track = (BrBfsTrack*)bfsTrackTable->At(t); //finding the track projections to the image plane BrDetectorTrack* t5Track = track->GetBackDetectorTrack(); BrLine3D trackLineT5 = t5Track->GetTrackLine(); BrLine3D trackLineG = fT5Volume->LocalToGlobal(trackLineT5); BrLine3D trackLineRich = fRichVolume->GlobalToLocal(trackLineG); BrLine3D reflectedLine = ReflectTrackLineOnMirror(trackLineRich); BrVector3D ringCenter = fImagePlane.GetIntersectionWithLine(reflectedLine); // get it in ring center coordinates ringCenter = InImagePlaneCoordinates(ringCenter); trackHits[t].x = ringCenter.X(); trackHits[t].y = ringCenter.Y(); trackHits[t].nr = t; trackHits[t].id = track->GetId(); trackHits[t].nhits = -1; trackHits[t].e = -1; trackHits[t].p = track->GetMomentum(); if ((fUseD3Momentum)&&(!fUseD4Momentum)) trackHits[t].p = track->GetBfsFrontTrack()->GetMomentum(); if ((fUseD4Momentum)&&(!fUseD3Momentum)) trackHits[t].p = track->GetBfsBackTrack()->GetMomentum(); } //copy hits in rdo to hits in array fNumberOfHits = richRdo->GetEntries(); RICHHIT hits[fNumberOfHits]; for(Int_t i=0; i<fNumberOfHits; i++) { BrChkvRdo::BrChkvHit* hit = (BrChkvRdo::BrChkvHit*)richRdo->GetHit(i); hits[i].x = fRichParams->GetTubePosition(hit->GetTubeNum()).X(); hits[i].y = fRichParams->GetTubePosition(hit->GetTubeNum()).Y(); hits[i].e = hit->GetEnergy(); hits[i].used = kFALSE; hits[i].tracknr = -1; hits[i].r = new Float_t[fNumberOfTracks]; for (Int_t t = 0; t<fNumberOfTracks; t++){ Float_t dX = hits[i].x - trackHits[t].x; Float_t dY = hits[i].y - trackHits[t].y; hits[i].r[t] = TMath::Sqrt(dX*dX + dY*dY); } } //---------------------------------------------------------------- // loop over tracks for (Int_t t=0; t<fNumberOfTracks; t++) { Float_t radius = FindRing(hits, trackHits[t]); Float_t momentum = trackHits[t].p; Float_t mass=-1; if (radius>1) mass = CalculateMass(momentum, radius); // create Rich Pid BrRichPid* richPid = new BrRichPid(); richPid->SetMass(mass); richPid->SetRingRadius(radius); richPid->SetRingHits(trackHits[t].nhits); // Here the summed energy signal should be added ... richPid->SetTrackId(trackHits[t].id); richPidTable->Add(richPid); if (HistOn()) { Snapshot(hits, trackHits[t]); fMass->Fill(mass); fRingCenters->Fill(trackHits[t].x,trackHits[t].y); fRadVsP->Fill(momentum, radius); fMassVsP->Fill(momentum, mass); fMassVsRad->Fill(radius, mass); if (fNumberOfTracks == 1) fRadVsP1->Fill(momentum, radius); else fRadVsP2->Fill(momentum, radius); } } //delete r array for(Int_t i=0; i<fNumberOfHits; i++) if (hits[i].r) delete [] hits[i].r; } //____________________________________________________________________ BrLine3D BrRichPidModule::ReflectTrackLineOnMirror(const BrLine3D& trackLine) { // Reflects the track line on the focal mirror. // // 1) find the reflection axis (the axis defined by the center of // the mirror sphere and the track-mirror intersection point). // // 2) find angle between track line and reflection axis // // 3) find the rotation axis (perpendicular to trackLine and // reflection axis) // // 4) rotate trackline with 2*angle around reflection axis BrVector3D intWithMirror = IntersectionWithMirror(trackLine); BrVector3D refAxis = intWithMirror - fMirrorSphereCenter; refAxis = refAxis.Unit(); BrVector3D trackDir = trackLine.GetDirection(); trackDir = trackDir.Unit(); Double_t angle = TMath::ACos(refAxis.Dot(trackDir)); BrVector3D rotAxis = trackDir.Cross(refAxis); // ok, there is no rotate method for BrVector3D, so I translate to // TVector3, make the rotation and translate back. TVector3 rotAxisR(rotAxis.X(),rotAxis.Y(),rotAxis.Z()); TVector3 newDirR(trackDir.X(), trackDir.Y(), trackDir.Z()); newDirR.Rotate(2 * angle, rotAxisR); BrVector3D newDir(newDirR.X(),newDirR.Y(),newDirR.Z()); return BrLine3D(intWithMirror,newDir); } //____________________________________________________________________ BrVector3D BrRichPidModule::IntersectionWithMirror(const BrLine3D& trackLine) { // Calculates the intersection between the trackLine and the focal // mirror of the RICH (Djamels calculations). // // find intersection point with sphere: // // line equations: sphere: // x = ax*t + xo (x - xs)^2 + (y - ys)^2 + (z - zs)^2 = r^2 // y = ay*t + yo // z = az*t + zo // // Equation to solve: // // (axt + xo - xs)^2 + (ayt + yo - ys)^2 + (azt + zo - zs)^2 - r^2 = 0 // // t^2(ax^2 + ay^2 + az^2) // + 2t(ax*(xo - xs) + ay*(yo - ys) + az*(zo - zs)) // + ((xo - xs)^2 + (yo - ys)^2 + (zo - zs)^2 - r^2) = 0 // // c1 t^2 + c2 t + c3 = 0; // new definitions: Double_t Xo = trackLine.GetOrigin().X() - fMirrorSphereCenter.X(); // = (xo - xs) Double_t Yo = trackLine.GetOrigin().Y() - fMirrorSphereCenter.Y(); // = (yo - ys) Double_t Zo = trackLine.GetOrigin().Z() - fMirrorSphereCenter.Z(); // = (zo - zs) Double_t ax = trackLine.GetDirection().X(); Double_t ay = trackLine.GetDirection().Y(); Double_t az = trackLine.GetDirection().Z(); Double_t radius = fMirrorSphereRadius; // ______________ // solve 2nd deg. equation : t = (-c2 + V c2^2 - 4c1c3 )/2c1 // (the negative solution is not good) Double_t c[3] = { ax*ax + ay*ay + az*az, 2*(ax*Xo + ay*Yo + az*Zo), Xo*Xo + Yo*Yo + Zo*Zo - radius*radius }; Double_t delta = c[1]*c[1] - 4*c[0]*c[2]; Double_t t = (-c[1] + TMath::Sqrt(delta))/2/c[0]; // intersection point: BrVector3D inters(ax*t + trackLine.GetOrigin().X(), ay*t + trackLine.GetOrigin().Y(), az*t + trackLine.GetOrigin().Z()); if (DebugLevel()>15) { cout << "BrRichPidModule::IntersectionWithmirror: " << endl << "Track intersection with mirror: " << endl << " x = " << inters.X() << endl << " y = " << inters.Y() << endl << " z = " << inters.Z() << endl; } return inters; } //____________________________________________________________________ BrVector3D BrRichPidModule::InImagePlaneCoordinates(const BrVector3D& hit) { // Hit must be on image plane. Return hit coordinates (x,y,0) on // Rich image plane. if (fImagePlane.GetMinimumDistanceFromPlane(hit)>0.01) { Warning("InImagePlaneCoordinates", "The hit is not on the imageplane - something is screwy!"); return BrVector3D(0,0,-1); } Float_t x = fImagePlaneCenter.X() - hit.X(); Float_t dy = hit.Y() - fImagePlaneCenter.Y(); Float_t dz = hit.Z() - fImagePlaneCenter.Z(); Float_t y; if (dy>0) y = TMath::Sqrt(dy*dy + dz*dz); else y = - TMath::Sqrt(dy*dy + dz*dz); return BrVector3D(x,y,0); } //____________________________________________________________________ Float_t BrRichPidModule::GuessRadius(RICHHIT* richHits, Int_t trackNumber) { // Tries to find the radius of the ring. // // Seach for hits with distance to the ring center in a dR window and // fills the number of hits into a histogram. // // Fits this histogram with a gaus to find the approximate maximum // (radius) , if the maximum of the histogram is too far from the // mean of the fit, simply use the max as radius. if (fNumberOfHits<fMinHitsInRing) return 0; TH1F* histE = new TH1F(Form("DEDR%d",fNDEDR), Form("DEDR, no.%d", fNDEDR), 40,-1,12); // loop over Rich hits for(Int_t i=0; i< fNumberOfHits; i++) { Float_t r = richHits[i].r[trackNumber]; if (r<fMaxRadius) { for (Int_t b=0; b<histE->GetNbinsX(); b++) { if (TMath::Abs(r - histE->GetBinCenter(b)) < fDRWindow/2) histE->Fill(histE->GetBinCenter(b),richHits[i].e); } } } // fit around maximum Float_t max = histE->GetBinCenter(histE->GetMaximumBin()); fGausFit->SetParameter(0, histE->GetMaximum()); fGausFit->SetParameter(1, max); fGausFit->SetParLimits(1, max-2,max+2); histE->Fit(fGausFit,"Q0","",max-2,max+2); Float_t radius = fGausFit->GetParameter(1); if ((TMath::Abs(radius-max)>2.)||(radius<1)) radius = max; if ((radius<0)||(radius>fMaxRadius)) radius = 0; if (HistOn()) { fNDEDR++; if (fNDEDR<200) { fHistDir->cd("DEDR"); fDEDR->Add(histE->Clone()); } } delete histE; return radius; } //____________________________________________________________________ Float_t BrRichPidModule::FindRing(RICHHIT* richHits, TRACKHIT track) { // Rind a ring with center at (cX, cY) from hits. // // 1) Requires nHits >= fMinHitsInRing // 2) Guesses a radius // 3) Checks that there is a reasonable number of hits in the ring // 4) If not, remove hits in ring from the hits array and try again // 5) Does this 5 times. If no rings found returns 0. // 6) If ring is found, refit the radius. Bool_t used[fNumberOfHits]; for (Int_t i=0; i<fNumberOfHits; i++) used[i] = richHits[i].used; if (fNumberOfHits<fMinHitsInRing) return 0; Float_t radiusSum = 0; Float_t radiusGuess = 0; Int_t nRingHits = 0; Float_t nPhotons = 0; Float_t estNPhotons = 0; Int_t counter = 0; Bool_t ringFound = kFALSE; // loop until ring found or counter > 5 while((!ringFound)&&(counter++<5)) { // guess radius, find corresponding hits and the estimated number // of photons for a ring of that radius. radiusGuess = GuessRadius(richHits, track.nr); for (Int_t i=0; i<fNumberOfHits; i++) { if ((richHits[i].used == kFALSE)&&(TMath::Abs(richHits[i].r[track.nr] - radiusGuess) < fDRWindow/2)) { richHits[i].used = kTRUE; nRingHits++; nPhotons = nPhotons + richHits[i].e; radiusSum = radiusSum + richHits[i].r[track.nr]; } } estNPhotons = EstimateNumberOfPhotons(radiusGuess); // if there is enough hits in ring accept it // due to saturation ,the rings cannot be expected to have more // than 8 hits. if (nRingHits>8) ringFound = kTRUE; else if ((nRingHits>(0.4*estNPhotons)&&(nRingHits>=fMinHitsInRing))) ringFound = kTRUE; else { nRingHits=0; nPhotons =0; radiusSum=0; } } // if no ring found clear the ringHits array and return 0 if (!ringFound) return 0; // taking the average radius as radius Float_t radius=0; if (nRingHits>0) radius = radiusSum/(Float_t)nRingHits; // Float_t phiRad = 0; // loop over hits for (Int_t i=0; i<fNumberOfHits; i++) { // restore old used hits richHits[i].used = used[i]; // geo hist if (fCheckGeometry && HistOn()) { phiRad = TMath::ATan2((richHits[i].y-track.y),(richHits[i].x-track.x)); fDRVsPhi->Fill(phiRad, (richHits[i].r[track.nr] - radiusGuess), richHits[i].e); if (radiusGuess!=0) fRing->Fill(TMath::Cos(phiRad)*richHits[i].r[track.nr]/radiusGuess, TMath::Sin(phiRad)*richHits[i].r[track.nr]/radiusGuess, richHits[i].e); } // find out what hits were used now and add them to used if (TMath::Abs(richHits[i].r[track.nr] - radiusGuess) < fDRWindow/2) { richHits[i].used = kTRUE; richHits[i].tracknr = track.nr; if (fCheckGeometry && HistOn()) { fDRVsPhiAcc->Fill(phiRad, (richHits[i].r[track.nr] - radiusGuess), richHits[i].e); if (radiusGuess!=0) fRingAcc->Fill(TMath::Cos(phiRad)*richHits[i].r[track.nr]/radiusGuess, TMath::Sin(phiRad)*richHits[i].r[track.nr]/radiusGuess , richHits[i].e); } } } track.nhits = nRingHits; track.e = nPhotons; if (HistOn()) { fHitsInRing->Fill(nRingHits); fHitsVsEstPhotons->Fill(estNPhotons,Float_t(nRingHits)); fHitsVsPhotons->Fill(nPhotons,Float_t(nRingHits)); fHitsVsRadius->Fill(radius,nRingHits); } return radius; } //____________________________________________________________________ Float_t BrRichPidModule::EstimateNumberOfPhotons(Float_t radius) { // The number of detected cherenkov photons is given by // // N = No * L sin^2(theta) // // where No is a detector quality parameter, L is the path length // and theta is the cherenkov angle Float_t sinTheta = TMath::Sin(TMath::ATan(radius/fMirrorFocalLength)); Float_t N = 55 * 150 * sinTheta * sinTheta; return N; } //____________________________________________________________________ Float_t BrRichPidModule::CalculateMass(Float_t p, Float_t radius) { // Calculates the theoretical mass from the momentum and the ring // radius. // ____________ // we have : 1/beta = V 1 + (m/p)^2 and cos(theta) = 1/(n*beta) // _____________________ // => m = p V (n*cos(theta))^2 - 1 // Double_t theta = TMath::ATan(radius/fMirrorFocalLength); Double_t ncostheta = fRefractiveIndex * TMath::Cos(theta); Double_t mass = TMath::Abs(p) * TMath::Sqrt(TMath::Abs(ncostheta*ncostheta -1)); return mass; } //____________________________________________________________________ void BrRichPidModule::Snapshot(RICHHIT* richHits, TRACKHIT track) { if (fNSnapshots<fMaxNSnapshots) { TDirectory* saveDir = gDirectory; fHistDir->cd(); //make the snapshot TH2F* hist = (TH2F*)fSnapshots->FindObject(Form("Snapshot%d",fNSnapshots)); fNSnapshots++; for(Int_t i=0; i<fNumberOfHits; i++) { hist->Fill(richHits[i].x,richHits[i].y,richHits[i].e); const Char_t* title = Form("Snapshot%03d (x=%03f,y=%06f,p=%06f)",fNSnapshots,track.x,track.y,track.p); hist->SetTitle(title); } gDirectory=saveDir; } else { if (Verbose()>15) Warning("Snapshot","NSnapshots > MaxNSnapshot. Image not created!"); } } //____________________________________________________________________ void BrRichPidModule::End() { // Run-level finalisation SetState(kEnd); } //____________________________________________________________________ void BrRichPidModule::Finish() { // Job-level finalisation SetState(kFinish); if (!HistOn()) return; TDirectory* saveDir = gDirectory; fHistDir->cd(); if (fCheckGeometry) { TProfile* profile = fDRVsPhiAcc->ProfileX("dRVsPhiAccProfile"); TF1* fit = new TF1("fit","[0]*cos(x +[1])"); fit->SetParameter(0,profile->GetMaximum()/2); profile->Fit(fit,"Q0","Q0",-3.14,3.14); Float_t xoff = fit->GetParameter(0) * TMath::Cos(fit->GetParameter(1)); Float_t yoff = fit->GetParameter(0) * TMath::Sin(fit->GetParameter(1)); if (Verbose()>4) { cout << "BrRichPidModule summary:" << endl; cout << "Rings off by approx " << setprecision(2) << fit->GetParameter(0) << " cm" << endl << "in direction " << setprecision(2) << fit->GetParameter(1) << " (rad)." << endl; cout << "Try new SetImagePlaneOffset(" << setprecision(3) << fImagePlaneOffset.GetX() + xoff << "," << setprecision(3) << fImagePlaneOffset.GetY() + yoff << ");" << endl; } } gDirectory=saveDir; } //____________________________________________________________________ void BrRichPidModule::Print(Option_t* option) const { // Print module information // See BrModule::Print for options. // In addition this module defines the Option: // <fill in here> TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Claus O. E. Jorgensen" << endl << " Last Modifications: " << endl << " $Author: ekman $" << endl << " $Date: 2002/06/10 13:30:57 $" << endl << " $Revision: 1.15 $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrRichPidModule.cxx,v $ // Revision 1.15 2002/06/10 13:30:57 ekman // Small fix to keep calculated masses < 0. It should be changed to m2 soon! // // Revision 1.14 2002/04/03 11:32:20 ekman // Added dR fit in finish method. When fCheckGeometry=kTRUE new imageplane // offsets are suggested. // // Revision 1.13 2002/03/06 21:11:26 ekman // Fixed bug in momentum sorting algorithm. // // Revision 1.12 2002/02/13 11:52:35 ekman // Added new geo histograms. Introduced sorting of tracks by momentum & other minor changes. // // Revision 1.11 2002/02/06 13:41:08 ekman // small bug fix // // Revision 1.10 2002/02/06 12:10:26 ekman // Added new geo (ring hist) // // Revision 1.9 2002/01/31 15:26:09 ouerdane // increased binning in histograms ring radius vs momentum and mass // // Revision 1.8 2001/11/19 15:03:22 ouerdane // check in Finish if histograms are booked before any actions on histograms, added TDirectory member in the C1 module // // Revision 1.7 2001/11/19 10:22:41 ekman // Removed Stop in event method (if no bfs track). // // Revision 1.6 2001/11/05 07:04:04 ouerdane // changed BFS for Bfs // // Revision 1.5 2001/10/27 04:14:30 ekman // Major changes. Arrays of tubehits and trackhits introduced. New geocheck histograms added. // // Revision 1.2 2001/10/02 01:54:49 ouerdane // removed named contructor calls to Pid data classes // // Revision 1.1 2001/09/24 13:15:03 ekman // Added BrRichPidModule. // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|