123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237 |
- //--------------------------------------------------------------------------
- //----- MpdTgemHitProducer -----
- //----- Created 8/12/09 by A. Zinchenko -----
- //--------------------------------------------------------------------------
- /** MpdTgemHitProducer.cxx
- *@author A.Zinchenko <Alexander.Zinchenko@jinr.ru>
- **
- ** MPD End-Cap Tracker (ECT) hit producer
- **/
- #include "MpdTgemHitProducer.h"
- #include "MpdTgemHit.h"
- #include "MpdTgemPoint.h"
- #include "FairRun.h"
- #include "FairRuntimeDb.h"
- #include "FairRootManager.h"
- #include "FairDetector.h"
- #include <TClonesArray.h>
- #include "TGeoManager.h"
- #include "TGeoVolume.h"
- #include "TGeoNode.h"
- #include "TGeoMatrix.h"
- #include <TMath.h>
- #include <TRandom.h>
- #include "TVector3.h"
- #include <iostream>
- // ----- Default constructor -------------------------------------------
- MpdTgemHitProducer::MpdTgemHitProducer(const char* fileGeo) :
- FairTask("Ideal Tgem hit Producer") {
- fFileGeo = fileGeo;
- eneThr = 0.001; // Energy threshold for Tgem
- //= hZ = new TH1F("hZ","This is the Z distribution",100,0,300)
- }
- // -------------------------------------------------------------------------
- // ----- Destructor ----------------------------------------------------
- MpdTgemHitProducer::~MpdTgemHitProducer() { }
- // -------------------------------------------------------------------------
- // ----- Public method Init --------------------------------------------
- InitStatus MpdTgemHitProducer::Init() {
-
- cout << " INITIALIZATION *********************" << endl;
-
- //FairDetector::Initialize();
- //FairRun* sim = FairRun::Instance();
- //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
-
- // Get RootManager
- FairRootManager* ioman = FairRootManager::Instance();
- if ( ! ioman ) {
- cout << "-E- MpdTgemHitProducer::Init: "
- << "RootManager not instantiated!" << endl;
- return kFATAL;
- }
-
- // Get input array
- fPointArray = (TClonesArray*) ioman->GetObject("TgemPoint");
- if ( ! fPointArray ) {
- cout << "-W- MpdTgemHitProducer::Init: "
- << "No TgemPoint array!" << endl;
- return kERROR;
- }
-
- // Create and register output array
- fDigiArray = new TClonesArray("MpdTgemHit");
-
- ioman->Register("TgemHit","TgemHits",fDigiArray,kTRUE);
- CreateStructure();
-
- cout << "-I- MpdTgemHitProducer: Intialization successfull" << endl;
-
- return kSUCCESS;
- }
- // -------------------------------------------------------------------------
- // ----- Public method Exec --------------------------------------------
- void MpdTgemHitProducer::Exec(Option_t* opt) {
-
- //cout << " DIGI EXECUTION *********************" << endl;
- static int eventCounter = 0;
- cout << " Event " << ++eventCounter << endl;
- // Reset output array
- if ( ! fDigiArray ) Fatal("Exec", "No DigiArray");
-
- fDigiArray->Clear();
-
- /*
- // Declare some variables
- MpdTgemPoint* point = NULL;
- map<Int_t, Float_t> fTrackEnergy;
- fTrackEnergy.clear();
- map<Int_t, Float_t>::const_iterator p;
-
- // Loop over ZdcPoints
- Int_t nPoints = fPointArray->GetEntriesFast();
- for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
- point = (MpdTgemPoint*) fPointArray->At(iPoint);
- fTrackEnergy[point->GetDetectorID()] += point->GetEnergyLoss();
- }
- // Loop over ZdcPoint
-
- // Loop to register ZdcHit
- for(p=fTrackEnergy.begin(); p!=fTrackEnergy.end(); ++p) {
- if ((*p).second>eneThr)
- AddHit(1, (*p).first, (*p).second);
- }
- */
- Bool_t mirror = kFALSE; //kTRUE; // add mirror hits if TRUE
- Int_t nPoints = fPointArray->GetEntriesFast(), lay = 0, nKH = 0, itube = 0;
- //++ cout<< "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! nPoints = "<< nPoints <<endl;
- //========================= test =============================
- //=+ MpdTgemPoint *h = (MpdTgemPoint*) fPointArray->UncheckedAt(2);
- //=+ cout<< "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! GetZ = "<< h->GetZ() <<" "<<"h->GetDetectorID()="<<h->GetDetectorID()<<endl;
- //============================================================
- /* ==
- Double_t phi, r, coord, errR = 0.2, errRphi = 0.02; // 2mm in R, 200um in R-Phi
- //Double_t phi, r, coord, errR = 0.01, errRphi = 0.01; // 100um in R, 100um in R-Phi
- Double_t firstHit[100][1000] = {{0},{0}};
- Int_t indxTube[100][1000];
- for (Int_t i = 0; i < 100; ++i) {
- for (Int_t j = 0; j < 1000; ++j) indxTube[i][j] = -1;
- }
- TVector3 pos, dpos(0.,0.,0.);
- for (Int_t ih = 0; ih < nPoints; ++ih) {
- MpdTgemPoint *h = (MpdTgemPoint*) fPointArray->UncheckedAt(ih);
- if (h->GetZ() < 0) continue; // keep only Z>0
- //if (TMath::Abs(TMath::ATan2(h->GetTrackY(),h->GetTrackX())*TMath::RadToDeg()-97.5) > 7.5) continue; // !!!
- //if (h->GetTrackID() != 1427) continue; // !!!
- phi = h->GetPhi(); // tube Phi
- lay = h->GetDetectorID() / 1000; // + 1;
- itube = h->GetDetectorID() % 1000; // tube number
- // Extrapolate track to Z = Ztube
- Double_t dZ = h->GetZ() - h->GetTrackZ();
- Double_t dt = dZ / h->GetPz();
- if (TMath::Abs(h->GetPz()) > 1.e-6 && h->GetPz() * h->GetZ() > 0) dt = dZ / h->GetPz();
- Double_t xNew = h->GetTrackX() + dt * h->GetPx();
- Double_t yNew = h->GetTrackY() + dt * h->GetPy();
- Double_t zNew = h->GetTrackZ() + dt * h->GetPz();
- pos.SetXYZ(xNew,yNew,zNew);
- //cout << dZ << " " << h->GetTrackX() << " " << xNew << " " << h->GetTrackY() << " " << yNew << " " << lay << endl;
- // Transform to the rotated coordinate system
- Double_t cosPhi = TMath::Cos(phi);
- Double_t sinPhi = TMath::Sin(phi);
- //Double_t xLoc = h->GetX() * cosPhi + h->GetY() * sinPhi; // cross-check
- //Double_t yLoc = -h->GetX() * sinPhi + h->GetY() * cosPhi;
- Double_t xRot = xNew * cosPhi + yNew * sinPhi;
- Double_t yRot = -xNew * sinPhi + yNew * cosPhi;
- //Double_t xLoc = (xNew - h->GetX()) * cosPhi + (yNew - h->GetY()) * sinPhi;
- //Double_t yLoc = -(xNew - h->GetX()) * sinPhi + (yNew - h->GetY()) * cosPhi;
- //r = xNew * xNew + yNew * yNew;
- //r = TMath::Sqrt (r);
- //r = TMath::Abs(xLoc);
- r = xRot;
- //cout << xRot << " " << yRot << " " << r << " " << h->GetPz() << endl;
- coord = yRot;
- Double_t yWire = -h->GetX() * sinPhi + h->GetY() * cosPhi;
- Double_t dY = TMath::Abs (yWire - coord);
- // Add error
- Double_t dRphi = 0, dR = 0;
- gRandom->Rannor(dRphi,dR);
- MpdTgemHit *hit = new ((*fDigiArray)[nKH]) MpdTgemHit(h->GetDetectorID(), pos, dpos, ih, 0);
- hit->SetMeas(coord+dRphi*errRphi);
- hit->SetError(errRphi);
- hit->SetPhi(phi);
- hit->SetIndex(ih);
- //MpdKalmanHitZ *hit = new ((*fKHits)[nKH]) MpdKalmanHitZ(r+dR*errR, phi, coord+dRphi*errRphi,
- // h->GetTrackZ(), errR, errRphi, lay, ih);
- //hit->SetXY(h->GetX(), h->GetY());
- // Add second measurement - just for test at the moment
- //!!!
- //hit->SetNofDim(2);
- //!!!
- ++nKH;
- }
- */
- //= fDigiArray->Compress();
- //= fDigiArray->Sort(); // in ascending order in abs(Z)
- //= cout << " Max layer = " << ((MpdTgemHit*)fDigiArray->Last())->GetLayer() << " " << fDigiArray->GetEntriesFast() << endl;
- }
- // -------------------------------------------------------------------------
- // ----- Public method Create Structure --------------------------------
- void MpdTgemHitProducer::CreateStructure() {
-
- /*
- TString work = getenv("VMCWORKDIR");
- work = work + "/geometry/" + fFileGeo;
- cout << "-I- <MpdTgemHitProducer::CreateStructure> Zdc geometry loaded from: "
- << work << endl;
- Int_t detId = -1;
- MpdTgemReader read(work);
-
- for(Int_t module=1; module<=read.GetMaxModules(); module++)
- for(Int_t row=1; row<=read.GetMaxRows(module); row++)
- for(Int_t crystal=1; crystal<=read.GetMaxCrystals(module,row); crystal++) {
- DataG4 data = read.GetData(module,row,crystal);
- for(Int_t copy=1; copy<=20; copy++) {
- detId = module*100000000 + row*1000000 + copy*10000 + crystal;
- emcX[detId] = data.posX; emcY[detId] = data.posY; emcZ[detId] = data.posZ;
- emcTheta[detId] = data.theta; emcTau[detId] = data.tau;
- if (module==3)
- emcPhi[detId] = fmod(data.phi+90.*(copy-1),360);
- else
- emcPhi[detId] = fmod(data.phi+22.5*(copy-1),360);
- }
- }
- */
- }
- // ----- Private method AddDigi --------------------------------------------
- MpdTgemHit* MpdTgemHitProducer::AddHit(Int_t trackID,Int_t detID, Float_t energy){
- // It fills the MpdTgemHit category
- // cout << "MpdTgemHitProducer: track " << trackID << " evt " << eventID << " sec " << sec << " plane " << pla << " strip " << strip << "box " << box << " tube " << tub << endl;
- TClonesArray& clref = *fDigiArray;
- Int_t size = clref.GetEntriesFast();
- return new(clref[size]) MpdTgemHit(); // FIXME: real hit info needed here
- }
- // ----
- ClassImp(MpdTgemHitProducer)
|