///////////////////////////////////////////////////////////// // CbmSttDigiProducer // // Class for digitalization for STT1 // // authors: Pablo Genova - Pavia University // Lia Lavezzi - Pavia University // ///////////////////////////////////////////////////////////// #include "TClonesArray.h" #include "CbmRootManager.h" #include "CbmSttDigiProducer.h" #include "CbmSttDigi.h" #include "CbmSttPoint.h" #include "TGeoManager.h" #include "TGeoVolume.h" #include "TGeoNode.h" #include "TGeoMatrix.h" #include "TVector3.h" #include "TStraw.h" // ----- Default constructor ------------------------------------------- CbmSttDigiProducer::CbmSttDigiProducer() : CbmTask("Ideal STT Digi Producer") { } // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- CbmSttDigiProducer::~CbmSttDigiProducer() { } // ------------------------------------------------------------------------- // ----- Public method Init -------------------------------------------- InitStatus CbmSttDigiProducer::Init() { // Get RootManager CbmRootManager* ioman = CbmRootManager::Instance(); if ( ! ioman ) { cout << "-E- CbmSttDigiProducer::Init: " << "RootManager not instantiated!" << endl; return kFATAL; } // Get input array fPointArray = (TClonesArray*) ioman->GetObject("SttPoint"); if ( ! fPointArray ) { cout << "-W- CbmSttDigiProducer::Init: " << "No SttPoint array!" << endl; return kERROR; } // Create and register output array fDigiArray = new TClonesArray("CbmSttDigi"); ioman->Register("SttDigi","Stt",fDigiArray,kTRUE); // Geometry loading TFile *tstfile=ioman->GetInFile(); TGeoManager *geoMan = (TGeoManager*) tstfile->Get("CBMGeom"); fVolumeArray = geoMan->GetListOfVolumes(); fVolumeArray = gGeoManager->GetListOfVolumes(); cout << "-I- CbmSttDigiProducer: Intialization successfull" << endl; return kSUCCESS; } // ------------------------------------------------------------------------- // ----- Public method Exec -------------------------------------------- void CbmSttDigiProducer::Exec(Option_t* opt) { Int_t evtn=0; CbmSttPoint *ptemp =(CbmSttPoint*) fPointArray->At(0); if(ptemp !=NULL) { evtn=ptemp->GetEventID(); if(evtn%50==0)cout << "Event Number "<Clear(); // Declare some variables CbmSttPoint* point = NULL; CbmSttPoint* point2 = NULL; // Loop over SttPoints Int_t nPoints = fPointArray->GetEntriesFast(); for (Int_t iPoint=0; iPointAt(iPoint); if (! point) continue; double InOut[6]; memset(InOut, 0, sizeof(InOut)); Int_t j=0; Bool_t flag=0; if(point->GetPostFlag()==1) { for(j=0;j<(nPoints-iPoint-1);j++) { if(((CbmSttPoint*) fPointArray->At(iPoint+j))->GetPreFlag()==1) { point2 = (CbmSttPoint*) fPointArray->At(iPoint+j); Bool_t ctrl=1; if(point->GetXOut()==-999||point->GetYOut()==-999||point->GetZOut()==-999||point2->GetXIn()==-999||point2->GetYIn()==-999||point2->GetZIn()==-999) { ctrl=0; } // point = point where the track is exiting the // mylar tube and entering the gas; // point2 = point where the track is exiting the // gas and entering the mylar tube if((point->GetNam()==point2->GetNam()) && (ctrl==1)) { InOut[0] = point->GetXOut(); InOut[1] = point->GetYOut(); InOut[2] = point->GetZOut(); InOut[3] = point2->GetXIn(); InOut[4] = point2->GetYIn(); InOut[5] = point2->GetZIn(); flag=1; } break; } } if(!flag) continue; // Geometry call --------------------------- TGeoVolume *vol = (TGeoVolume *) fVolumeArray->At(0); TGeoNode *tmpnod; TGeoMatrix *geoMat; double cosskew = 999; double arccos = 999; TVector3 rot1,rot2,rot3; Double_t cpos[3]; tmpnod = (TGeoNode *) vol->GetNode(point->GetNam()); geoMat = (TGeoMatrix *) tmpnod->GetMatrix(); // position of the centre of the tube cpos[0]=geoMat->GetTranslation()[0]; // cm! cpos[1]=geoMat->GetTranslation()[1]; cpos[2]=geoMat->GetTranslation()[2]; // rotation matrices of the tube rot1.SetXYZ(geoMat->GetRotationMatrix()[0], geoMat->GetRotationMatrix()[1],geoMat->GetRotationMatrix()[2]); rot2.SetXYZ(geoMat->GetRotationMatrix()[3], geoMat->GetRotationMatrix()[4],geoMat->GetRotationMatrix()[5]); rot3.SetXYZ(geoMat->GetRotationMatrix()[6], geoMat->GetRotationMatrix()[7],geoMat->GetRotationMatrix()[8]); // single straw tube simulation ----------------------- TStraw stt; //setting the single straw tube simulation constants // 3 options currently available: // TConst(tube radius (cm), gas pressure (bar), Ar%, CO2%) // stt.TConst(0.4, 1, 0.9, 0.1); stt.TConst(0.5, 1, 0.9, 0.1); // stt.TConst(0.5, 2, 0.8, 0.2); // wire position: // rotation applied to (0, 0, L/2) // L = 1500 mm = tube length TVector3 half_length(rot1.Z() * 75, rot2.Z() * 75, rot3.Z() * 75); TVector3 second_extr = cpos + half_length; Int_t cm=1; stt.PutWireXYZ(cpos[0]/cm,cpos[1]/cm,cpos[2]/cm,second_extr[0]/cm,second_extr[1]/cm,second_extr[2]/cm); // skew angle TVector3 diff1 = second_extr - cpos; TVector3 diff2(0,0,1); cosskew = (second_extr[2] - cpos[2])/sqrt((second_extr[0] - cpos[0])*(second_extr[0] - cpos[0]) + (second_extr[1] - cpos[1])*(second_extr[1] - cpos[1]) + (second_extr[2] - cpos[2])*(second_extr[2] - cpos[2])); arccos = acos(cosskew)*180/TMath::Pi(); // TO BE SET AS DEBUG // cout << "skew: " << arccos << " cosskew: " << cosskew << endl; //cout << "name: " << namelog << std::endl; //std::cout << "Translation x, y, z: " << pos.getX() << " " << pos.getY() << " " << pos.getZ() << std::endl; //std::cout << "Rotation xx, xy, xz: " << rot.xx() << " " << rot.xy() << " " << rot.xz() << std::endl; //std::cout << "Rotation yx, yy, yz: " << rot.yx() << " " << rot.yy() << " " << rot.yz() << std::endl; //std::cout << "Rotation zx, zy, zz: " << rot.zx() << " " << rot.zy() << " " << rot.zz() << std::endl; //std::cout << "cosskew: " << cosskew << std::endl; // get particle momentum TVector3 momentum(point->GetPxOut(),point->GetPyOut(),point->GetPzOut()); // GeV Double_t GeV=1.; // position in cm (already in cm); momentum in GeV (already in GeV); mass in GeV (already in GeV) // drift time calculation Double_t pulset = stt.PartToTime(point->GetMass()/GeV, momentum.Mag()/GeV, InOut); // simulated radius (cm) double radius = stt.TimnsToDiscm(pulset); // true radius (cm) double true_rad = stt.TrueDist(InOut); // dE calculation double depCharge = stt.PartToADC(); // dE/dx calculation //TVector3 diff3=point1->Get //double distance =diff3.Mag(); // //double dedx = 999; //if (distance != 0) dedx = depCharge/(1000000 * distance); // in arbitrary units TVector3 center(cpos[0],cpos[1],cpos[2]); TVector3 tubemin = cpos-half_length; TVector3 tubemax = cpos+half_length; // create digi AddDigi(point->GetTrackID(),point->GetEventID(), pulset, radius, true_rad, point->GetNam(), center,tubemax, tubemin); //--------------------------- end geometry call } }// Loop over MCPoints // Event summary // cout << "-I- CbmSttDigiProducer: " << nPoints << " SttPoints, " // << nPoints << "Digis created." << endl; } // ------------------------------------------------------------------------- // ----- Private method AddDigi -------------------------------------------- CbmSttDigi* CbmSttDigiProducer::AddDigi(Int_t trackID,Int_t eventID, Double_t p, Double_t rsim, Double_t rtrue,TString nam,TVector3 center,TVector3 tubemax,TVector3 tubemin){ // see CbmDigi for digit description TClonesArray& clref = *fDigiArray; Int_t size = clref.GetEntriesFast(); //cout << "-I- CbmSttDigiProducer: Adding Digi: track"<