|
//____________________________________________________________________ // // // //____________________________________________________________________ // $Id: BrC1DigModule.cxx,v 1.1.1.1 2001/06/21 14:55:06 hagel Exp $ // $Author: hagel $ // $Date: 2001/06/21 14:55:06 $ // $Copyright: 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov> // #ifndef BRAT_BrC1DigModule #include "BrC1DigModule.h" #endif #ifndef ROOT_TMath #include "TMath.h" #endif #ifndef ROOT_TDirectory #include "TDirectory.h" #endif #ifndef ROOT_TRandom #include "TRandom.h" #endif #ifndef BRAT_BrCoordinateSystem #include "BrCoordinateSystem.h" #endif #ifndef BRAT_BrRotMatrix #include "BrRotMatrix.h" #endif #ifndef BRAT_BrGeometryDbManager #include "BrGeometryDbManager.h" #endif #ifndef BRAT_BrParameterDbManager #include "BrParameterDbManager.h" #endif #ifndef BRAT_BrGeantHit #include "BrGeantHit.h" #endif #ifndef BRAT_BrGeantTrack #include "BrGeantTrack.h" #endif #ifndef ROOT_TH1 #include "TH1.h" #endif #ifndef ROOT_TH2 #include "TH2.h" #endif //____________________________________________________________________ ClassImp(BrC1DigModule); //____________________________________________________________________ BrC1DigModule::BrC1DigModule() { // Default constructor. DO NOT USE SetState(kSetup); } //____________________________________________________________________ BrC1DigModule::BrC1DigModule(const Char_t* name, const Char_t* title) : BrModule(name, title) { // Named constructor SetState(kSetup); fMeanNoPhotons = 0; fNoPhotons = 0; fCosCkvAngle = 0; fTopMirrorPlane = 0; fBotMirrorPlane = 0; SetGasRefractionIndex(); SetLightCollectEfficiency(); SetQuantumEfficiency(); SetMirrorReflectivity(); SetLambdaMin(); SetLambdaMax(); fNoPmts = 32; // should get it from Db or DetectorParameters.txt } //____________________________________________________________________ void BrC1DigModule::DefineHistograms() { // Initialiser at job level if (GetState() != kInit) { Stop("DefineHistograms", "Must be called after Init"); return; } TDirectory* saveDir = gDirectory; TDirectory* histDir = gDirectory->mkdir(GetName()); histDir->cd(); // Make histograms here fTracks = new TList(); fHits = new TList(); fCkv = new TList(); fTracks->Add(new TH1F("pi", "Track Momentum (GeV/c)", 150, 0, 15)); fTracks->Add(new TH1F("po", "Track Momentum above threshold(GeV/c)", 150, 0, 15)); fTracks->Add(new TH1F("mi", "Track Mass (GeV/c^2)", 200, 0, 2.5)); fTracks->Add(new TH1F("mo", "Track Mass (remaining) (GeV/c^2)", 200, 0, 2.5)); fTracks->Add(new TH1F("L", "Track Length (cm)", 400, 0, 100)); fHits->Add(new TH1F("hs", "hit statistics", 51, -0.5, 50.5)); fHits->Add(new TH1F("xi", "Local Xin Hit position", 200, -30, 30)); fHits->Add(new TH1F("yi", "Local Yin Hit position", 200, -30, 30)); fHits->Add(new TH1F("zi", "Local Zin Hit position", 400, -100, 100)); fHits->Add(new TH1F("xo", "Local Xout Hit position", 200, -30, 30)); fHits->Add(new TH1F("yo", "Local Yout Hit position", 200, -30, 30)); fHits->Add(new TH1F("zo", "Local Zout Hit position", 400, -100, 100)); fCkv->Add(new TH1F("mgp", "Mean generated photons", 100, 0, 100)); fCkv->Add(new TH1F("pgp", "Generated photons (Poisson)", 101, -0.5, 100.5)); fCkv->Add(new TH1F("ca", "Cherenkov angle", 200, 0, 10)); fCkv->Add(new TH2F("cab", "Ckv angle vs Beta", 200, 0.5, 1.1, 200, 0, 10)); fCkv->Add(new TH2F("cap", "Ckv angle vs P", 200, 0, 20, 200, 0, 10)); fCkv->Add(new TH2F("gpp", "Number of photons vs P", 150, 0, 15, 101, -0.5, 100.5)); fCkv->Add(new TH2F("cam", "Ckv angle vs Mass", 200, 0, 2.5, 200, 0, 10)); fCkv->Add(new TH1F("mir", "Mirror statistics", 31, -0.5, 2.5)); fCkv->Add(new TH1F("pmt", "PMT statistics", fNoPmts, 0.5, fNoPmts+0.5)); Double_t c1Width = fC1Volume->GetSize()[0]; fCkv->Add(new TH2F("tpmt2d", "Top PMT 2D statistics", 200, -c1Width/2, c1Width/2, 200, 20, 40)); fCkv->Add(new TH2F("bpmt2d", "Bot PMT 2D statistics", 200, -c1Width/2, c1Width/2, 200, 20, 40)); fCkv->Add(new TH2F("pxy", "Photon X vs Y ", 200, -1.5, 1.5, 200, -1.5, 1.5)); for (Int_t i = 1; i <= fNoPmts; i++) { fCkv->Add(new TH1F(Form("cp%d", i), Form("Photons collected by PMT %d", i), 51, -0.5, 50.5)); fCkv->Add(new TH2F(Form("tcp%d", i), Form("PMT photons %d vs momentum", i), 200, 0, 20, 51, -0.5, 50.5)); } gDirectory = saveDir; } //____________________________________________________________________ void BrC1DigModule::Init() { // Initialiser at job level SetState(kInit); SetGeometry(); } //____________________________________________________________________ void BrC1DigModule::SetGeometry() { // all things are in the C1 coordinate system BrGeometryDbManager* geoDb = BrGeometryDbManager::Instance(); fC1Volume = (BrDetectorVolume*)geoDb->GetDetectorVolume("BrDetectorVolume","C1"); // ---------------------------------------- // PMT geometry // set X, Z for PMTs (for the moment, don't need Y positions) // this piece of code is stolen from gbrahms (eg_c1.F) Float_t fPmtSpacing = 6.51; Float_t fC1Height = fC1Volume->GetSize()[1]; Float_t fC1Length = fC1Volume->GetSize()[2]; Float_t fC1Width = fC1Volume->GetSize()[0]; for (Int_t j = 1; j <=2; j++) { Float_t z = fC1Length/2. - (2*j - 3)*fPmtSpacing/2. - 9.; Float_t y = 0; for (Int_t i = 1; i <=8; i++) { Float_t x = (i - 4.5)*fPmtSpacing; Int_t n = (j - 1)*8 + i; BrVector3D topPos(x, y, z); BrVector3D botPos(x, -y, z); fPmtPosition[n-1] = topPos; // 0 to 15 fPmtPosition[n-1 + 16] = botPos; // 16 to 31 } } for (Int_t i = 0; i < fNoPmts; i++) cout << " PMT " << i+1 << ": " << fPmtPosition[i] << endl; fPmtRadius = 2.54; // stolen from bd_c1.F //---------------------------- // mirror stuff // set mirror center positions Float_t fMirrorOffset = 1.; Float_t fMirrorThickness = 0.6; Float_t x = 0; Float_t y = fC1Height/4.; Float_t z = fC1Length/2. - fMirrorOffset - (fC1Height/4. + fMirrorThickness/TMath::Sqrt(2.)); // 1st mirror (top: +y) BrVector3D pt1(x, y, z - fMirrorThickness/TMath::Sqrt(2.)); BrVector3D pt2(x + 10, y, z - fMirrorThickness/TMath::Sqrt(2.)); BrVector3D pt3(x, y + 10, z + 10 - fMirrorThickness/TMath::Sqrt(2.)); // 2nd mirror (top: -y) BrVector3D pt4(x, -y, z - fMirrorThickness/TMath::Sqrt(2.)); BrVector3D pt5(x + 10, -y, z - fMirrorThickness/TMath::Sqrt(2.)); BrVector3D pt6(x, -y - 10, z + 10 - fMirrorThickness/TMath::Sqrt(2.)); fTopMirrorPlane = new BrPlane3D(pt1, pt2, pt3); fBotMirrorPlane = new BrPlane3D(pt4, pt5, pt6); // C1 planes // top BrVector3D pt7(0, fC1Height/2, fC1Length/2); BrVector3D pt8(0, fC1Height/2, 0); BrVector3D pt9(fC1Width/2, fC1Height/2, fC1Length/2); // bot BrVector3D pt10(0, -fC1Height/2, fC1Length/2); BrVector3D pt11(0, -fC1Height/2, 0); BrVector3D pt12(fC1Width/2, -fC1Height/2, fC1Length/2); fTopC1Plane = new BrPlane3D(pt7, pt8, pt9); fBotC1Plane = new BrPlane3D(pt10, pt11, pt12); } //____________________________________________________________________ void BrC1DigModule::Begin() { // Initialiser at run level SetState(kBegin); } //____________________________________________________________________ void BrC1DigModule::Event(BrEventNode* inNode, BrEventNode* outNode) { // Main method. Process one event. SetState(kEvent); BrDataTable* c1Hits = inNode->GetDataTable("GeantHits C1VA"); if (!c1Hits) { if (DebugLevel()) Warning("BrC1DigModule::Event()", "No C1 hits in this event"); return; } // will check if hit has a momentum below Cerenkov threshold // I'm using Pthreshold = mass / sqrt(n^2 - 1) // derived from 1 = 1/(beta*n) // I prefer using the momentum for relativistic particles Int_t nhits = c1Hits->GetEntries(); if (HistOn()) ((TH1F*)fHits->FindObject("hs"))->Fill(nhits); for (Int_t i = 0; i < nhits; i++) { BrGeantHit* hit = (BrGeantHit*)c1Hits->At(i); Double_t P = hit->GetTrack()->P(); Double_t m = hit->GetTrack()->Mass(); if (HistOn()) { ((TH1F*)fTracks->FindObject("pi"))->Fill(P); ((TH1F*)fTracks->FindObject("mi"))->Fill(m); } Double_t PThreshold = m/ TMath::Sqrt(fGasRefractionIndex*fGasRefractionIndex - 1); if (P < PThreshold) continue; if (HistOn()) { ((TH1F*)fTracks->FindObject("po"))->Fill(P); ((TH1F*)fTracks->FindObject("mo"))->Fill(m); } GeneratePhotons(hit); PropagatePhotons(); CollectPmtLight(); // diagnostic histos if (HistOn()) { ((TH2F*)fCkv->FindObject("gpp"))->Fill(P, fNoPhotons); for (Int_t i = 0; i < fNoPmts; i++) { Int_t pmt = i + 1; ((TH1F*)fCkv->FindObject(Form("cp%d", pmt)))-> Fill(fNoCollectedPhotons[i]); ((TH2F*)fCkv->FindObject(Form("tcp%d", pmt)))-> Fill(P, fNoCollectedPhotons[i]); ((TH1F*)fCkv->FindObject("pmt"))-> Fill(pmt, (Stat_t)fNoCollectedPhotons[i]); } } delete [] fPhotonLine; delete [] fStruckMirror; } } //____________________________________________________________________ void BrC1DigModule::GeneratePhotons(BrGeantHit* hit) { // private method // generation of cerenkov photons // here, I assume that the charge is Z = +/- 1 // in this case, the number of photons is: // _ _ // dN | 1 | 1 // ------- = C * L * | 1 - ---------- | * -------- // dlambda |_ beta^2*n^2 _| lambda^2 // C = 2*Pi*alpha*Z^2 = 2*Pi / 137 = 4.586267e-2 // this is for the moment an extremely simple model, highly // unrealistic since I don't know the effective path length of the particle // through the gas. I only know that between the entrance point and // exit point of hit, there was the same track. I'll use these // points as a start // L is the path length of the particle through the gaz // beta is the velocity of the particle // n is the refraction index of the gaz. // // Note: in principle, n depends on lambda but I will assume that in // the lambda range I'm considering, n is constant. // ----------------------------------------------------------- // evaluate path length of the particle (vectors are in C1 c.s.) BrVector3D hitPosIn(hit->LocalPosIn()); BrVector3D hitPosOut(hit->LocalPosOut()); fPartPath = hitPosOut - hitPosIn; Double_t L = fPartPath.Norm(); // in cm if (HistOn()) { ((TH1F*)fTracks->FindObject("L"))->Fill(L); ((TH1F*)fHits->FindObject("xi"))->Fill(hitPosIn.GetX()); ((TH1F*)fHits->FindObject("yi"))->Fill(hitPosIn.GetY()); ((TH1F*)fHits->FindObject("zi"))->Fill(hitPosIn.GetZ()); ((TH1F*)fHits->FindObject("xo"))->Fill(hitPosOut.GetX()); ((TH1F*)fHits->FindObject("yo"))->Fill(hitPosOut.GetY()); ((TH1F*)fHits->FindObject("zo"))->Fill(hitPosOut.GetZ()); } // some constants Double_t C = 2*TMath::Pi()/137.; Double_t beta = hit->GetTrack()->Beta(); Double_t P = hit->GetTrack()->P(); // Cerenkov angle fCosCkvAngle = 1/(beta*fGasRefractionIndex); Double_t sinCkv = TMath::Sqrt(1 - fCosCkvAngle*fCosCkvAngle); if (HistOn()) { ((TH1F*)fCkv->FindObject("ca"))-> Fill(TMath::ACos(fCosCkvAngle)*180/TMath::Pi()); ((TH2F*)fCkv->FindObject("cab"))-> Fill(beta, TMath::ACos(fCosCkvAngle)*180/TMath::Pi()); ((TH2F*)fCkv->FindObject("cap"))-> Fill(P, TMath::ACos(fCosCkvAngle)*180/TMath::Pi()); Double_t mass = hit->GetTrack()->Mass(); ((TH2F*)fCkv->FindObject("cam"))-> Fill(mass, TMath::ACos(fCosCkvAngle)*180/TMath::Pi()); } // mean number of photons in lambda range fMeanNoPhotons = C*L*(1 - fCosCkvAngle*fCosCkvAngle)* (1/fLambdaMin - 1/fLambdaMax); // number of photons Poisson distributed fNoPhotons = gRandom->Poisson(fMeanNoPhotons); if (HistOn()) { ((TH1F*)fCkv->FindObject("mgp"))->Fill(fMeanNoPhotons); ((TH1F*)fCkv->FindObject("pgp"))->Fill(fNoPhotons); } // create photon lines BrVector3D null; fPhotonLine = new BrLine3D [fNoPhotons] (null, null); fStruckMirror = new Int_t [fNoPhotons]; // distribute photons uniformly randomly along particle path // Note: realistic ?? // line between entrance and exit points: BrLine3D line(hitPosIn, fPartPath); for (Int_t i = 0; i < fNoPhotons; i++) { fStruckMirror[i] = 0; // set origin point of the photon Double_t t = gRandom->Uniform(1); BrVector3D origin(fPartPath.GetX()*t + hitPosIn.GetX(), fPartPath.GetY()*t + hitPosIn.GetY(), fPartPath.GetZ()*t + hitPosIn.GetZ()); // set direction of the line (with a random angle since all cone // lines are equivalent) // direction of cone: BrVector3D dcone(fPartPath.GetX()/fPartPath.Norm(), fPartPath.GetY()/fPartPath.Norm(), fPartPath.GetZ()/fPartPath.Norm()); // define rotation matrix so that the cone direction is the new X axis // needs to be checked!! I might screw up the Euler angles Double_t thetaE = TMath::ATan(dcone.GetY()/dcone.GetX())*180./TMath::Pi(); Double_t phiE = 90.; Double_t psiE = TMath::ASin(-sinCkv)*180./TMath::Pi(); BrRotMatrix rot(thetaE, phiE, psiE); // now, we can work in this new frame and go back to the old BrCoordinateSystem coneCS(origin, rot); Double_t phi = gRandom->Uniform(TMath::Pi()); // random dir on cone BrVector3D photonDir(1, // X in cone CS sinCkv/fCosCkvAngle*TMath::Cos(phi), sinCkv/fCosCkvAngle*TMath::Sin(phi)); BrLine3D phl(origin, coneCS.TransformToMaster(photonDir)); fPhotonLine[i] = phl; // hopefully right! if (HistOn()) ((TH2F*)fCkv->FindObject("pxy"))->Fill(phl.GetDirection().GetX(), phl.GetDirection().GetY()); } } //____________________________________________________________________ void BrC1DigModule::PropagatePhotons() { // private method // photon are propagated to the mirror planes // ------------------------------------------- // intersect photon lines with mirror planes // Note: ignoring case where the geant hit exits behind the mirrors // there might be a small contamination of photons created behind // (cf histograms) but the reflection should propagate them in a // zone without PMTs...maybe for (Int_t i = 0; i < fNoPhotons; i++) { fStruckMirror[i] = 0; BrVector3D intersTop = fTopMirrorPlane->GetIntersectionWithLine(fPhotonLine[i]); BrVector3D intersBot = fBotMirrorPlane->GetIntersectionWithLine(fPhotonLine[i]); //--------------------------------------------- // check if intersections are on mirror surface //--------------------------------------------- // mirror height (Y coordinate in C1 c.s): C1 height/4 // mirror width (X coordinate in C1 c.s): C1 width Float_t mHeight = fC1Volume->GetSize()[1]/2; Float_t mWidth = fC1Volume->GetSize()[0]/2; Float_t mYThick = 0.6/TMath::Sqrt(2.); // stolen from geant if (intersTop.GetX() > -mWidth && intersTop.GetX() < mWidth && intersTop.GetY() > mYThick && intersTop.GetY() < mHeight) { fStruckMirror[i] = 1; fPhotonLine[i].SetOrigin(intersTop); } if (intersBot.GetX() > -mWidth && intersBot.GetX() < mWidth && intersBot.GetY() < -mYThick && intersBot.GetY() > -mHeight) { fStruckMirror[i] = 2; fPhotonLine[i].SetOrigin(intersBot); } if (HistOn()) ((TH1F*)fCkv->FindObject("mir"))->Fill(fStruckMirror[i]); if (fStruckMirror[i] != 0) ReflectPhoton(i); // making a real reflection } } //____________________________________________________________________ void BrC1DigModule::ReflectPhoton(Int_t p) { // private method // reflect a photon on a plane given a incoming direction and // intersection point on this plane // maybe it's easier to ignore this and consider the reflection of // the winston cones instead... // but since the mirrors are at 45 deg., the reflection of the // photon line is very easy (if I don't make mistake) BrVector3D dir = fPhotonLine[p].GetDirection(); // the direction in X remains the same Double_t y = dir.GetY(); Double_t z = dir.GetZ(); if (fStruckMirror[p] == 1) { dir.SetZ(y); dir.SetY(z); } else { dir.SetZ(-y); dir.SetY(-z); } fPhotonLine[p].SetDirection(dir); } //____________________________________________________________________ void BrC1DigModule::CollectPmtLight() { // private method // Collect photons on PMT. // For the moment, I simply project the photon lines on top and // bottom planes of C1 and check which PMT might be hit // will refine this with real cell that should be considered // if a photon doesn't reach a Pmt, I ignore it //-------- for (Int_t i = 0; i < fNoPmts; i++) fNoCollectedPhotons[i] = 0; for (Int_t p = 0; p < fNoPhotons; p++) { BrVector3D inters; if (fStruckMirror[p] == 1) // top mirror & PMTs inters = fTopC1Plane->GetIntersectionWithLine(fPhotonLine[p]); else if (fStruckMirror[p] == 2) // bot mirror & PMTs inters = fBotC1Plane->GetIntersectionWithLine(fPhotonLine[p]); else continue; Int_t pi = 0; Int_t pf = 16; TString s = "tpmt2d"; if (fStruckMirror[p] == 2) { pi = 16; pf = 32; s = "bpmt2d"; } // check if this point is in pmt disk (should be changed since // I think we have to deal with squares...not sure) for (Int_t i = pi; i < pf; i++) { Double_t x = inters.GetX() - fPmtPosition[i].GetX(); Double_t z = inters.GetZ() - fPmtPosition[i].GetZ(); if (x*x + z*z < fPmtRadius*fPmtRadius) { fNoCollectedPhotons[i]++; if (HistOn()) ((TH2F*)fCkv->FindObject(s.Data()))-> Fill(inters.GetX(), inters.GetZ()); break; } } } } //____________________________________________________________________ void BrC1DigModule::End() { // Run-level finaliser. SetState(kEnd); } //____________________________________________________________________ void BrC1DigModule::Finish() { // Run-level finaliser. SetState(kFinish); } //____________________________________________________________________ void BrC1DigModule::Print(Option_t* option) const { // Module Information method // // Options: (see also BrModule::Print) // TString opt(option); opt.ToLower(); BrModule::Print(option); if (opt.Contains("d")) cout << endl << " Original author: Djamel Ouerdane" << endl << " Revisted by: $Author: hagel $" << endl << " Revision date: $Date: 2001/06/21 14:55:06 $" << endl << " Revision number: $Revision: 1.1.1.1 $ " << endl << endl << "-------------------------------------------------" << endl; } //____________________________________________________________________ // // $Log: BrC1DigModule.cxx,v $ // Revision 1.1.1.1 2001/06/21 14:55:06 hagel // Initial revision of brat2 // // Revision 1.1 2001/06/19 12:55:46 ouerdane // Added class BrC1DigModule before the transition to brat2. // It still needs some debugging with the geometry manipulation but // compiles and works fine otherwise with geant simulation (cdat files) // // |
||||||
|
This page automatically generated by script docBrat by Christian Holm |
Copyright ; 2002 BRAHMS Collaboration
<brahmlib@rcf.rhic.bnl.gov>
|