MpdSftHitProducer.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429
  1. // --------------------------------------------------------------------------
  2. // ----- Class MpdSftHitProducer ------
  3. // ----- Created by E. Cordier 14/09/05 ------
  4. // ----- Modified by D. Gonzalez-Diaz 07/09/06 ------
  5. // ----- Modified by D. Gonzalez-Diaz 02/02/07 ------
  6. // --------------------------------------------------------------------------
  7. #include <iostream>
  8. #include "TSystem.h"
  9. #include "TClonesArray.h"
  10. #include "TRandom.h"
  11. #include "TString.h"
  12. #include "TVector3.h"
  13. #include "TMath.h"
  14. #include "MpdSftHitProducer.h"
  15. #include "FairMCApplication.h"
  16. #include "FairDetector.h"
  17. #include "FairRootManager.h"
  18. #include "MpdSftPoint.h"
  19. #include "MpdSftHit.h"
  20. //#include "FairMCTrack.h"
  21. using namespace std;
  22. // ---- Default constructor -------------------------------------------
  23. MpdSftHitProducer::MpdSftHitProducer()
  24. {
  25. fVerbose = 1;
  26. fSigmaT = 0.;
  27. fSigmaXY = 0.;
  28. fSigmaY = 0.;
  29. fSigmaZ = 0.;
  30. fNHits = -1;
  31. }
  32. // ---- Constructor ----------------------------------------------------
  33. MpdSftHitProducer::MpdSftHitProducer(const char *name, Int_t verbose)
  34. :FairTask(name,verbose)
  35. {
  36. fVerbose = verbose;
  37. fSigmaT = 0.;
  38. fSigmaXY = 0.;
  39. fSigmaY = 0.;
  40. fSigmaZ = 0.;
  41. fNHits = -1;
  42. }
  43. // ---- Destructor ----------------------------------------------------
  44. MpdSftHitProducer::~MpdSftHitProducer()
  45. {
  46. // FairRootManager *fManager =FairRootManager::Instance();
  47. // fManager->Write();
  48. }
  49. // ---- Init ----------------------------------------------------------
  50. InitStatus MpdSftHitProducer::Init()
  51. {
  52. FairRootManager *fManager = FairRootManager::Instance();
  53. // fTofPoints = (TClonesArray *) fManager->ActivateBranch("SftPoint");
  54. fTofPoints = (TClonesArray *) fManager->GetObject("SftPoint"); // EL
  55. //fMCTracks = (TClonesArray *) fManager->ActivateBranch("MCTrack");
  56. FILE *par;
  57. //Reading the parameter file. In the future this must be done in a different way.
  58. char header='#', cell_type='#';
  59. int region, module, cell;
  60. Int_t nregions=10, nmodules=500, ncells=500;
  61. Float_t X_tmp, Y_tmp, Dx_tmp, Dy_tmp;
  62. // Initialize the matrixes [make this index visible in all the macro]. FIXME
  63. for(int i=0;i<nregions;i++){
  64. for(int j=0;j<nmodules;j++){
  65. for(int k=0;k<ncells;k++){
  66. X[i][j][k] = -1;
  67. Y[i][j][k] = -1;
  68. Dx[i][j][k]= -1;
  69. Dy[i][j][k]= -1;
  70. }
  71. }
  72. }
  73. TString tofGeoFile = gSystem->Getenv("VMCWORKDIR");
  74. tofGeoFile += "/parameters/tof/tof_standard.geom.par";
  75. par=fopen(tofGeoFile,"r");
  76. if(par==NULL) {
  77. printf("\n ERROR WHILE OPENING THE PARAMETER FILE IN SFT HIT PRODUCER!");
  78. return kFATAL;
  79. }
  80. //Skip the header. In the future the header structure must be defined. FIXME
  81. while (fscanf(par,"%c",&header)>=0){
  82. if((int)(header-'0')==0) break;
  83. }
  84. //Read the first line
  85. region=0;
  86. fscanf(par,"%d%d%s%f%f%f%f", &module, &cell, &cell_type, &X_tmp, &Y_tmp, &Dx_tmp, &Dy_tmp);
  87. X[region][module][cell] = X_tmp;
  88. Y[region][module][cell] = Y_tmp;
  89. Dx[region][module][cell] = Dx_tmp;
  90. Dy[region][module][cell] = Dy_tmp;
  91. type[region][module][cell] = cell_type;
  92. //Read all the lines
  93. while(fscanf(par,"%d%d%d%s%f%f%f%f", &region, &module, &cell, &cell_type, &X_tmp, &Y_tmp, &Dx_tmp, &Dy_tmp)>=0){
  94. X[region][module][cell] = X_tmp;
  95. Y[region][module][cell] = Y_tmp;
  96. Dx[region][module][cell] = Dx_tmp;
  97. Dy[region][module][cell] = Dy_tmp;
  98. type[region][module][cell] = cell_type;
  99. }
  100. fclose(par);
  101. fHitCollection = new TClonesArray("MpdSftHit");
  102. fManager->Register("SftHit","Sft",fHitCollection, kTRUE);
  103. cout << "-I- MpdSftHitProducer: Intialization successfull" << endl;
  104. return kSUCCESS;
  105. }
  106. // ---- Exec ----------------------------------------------------------
  107. void MpdSftHitProducer::Exec(Option_t * option)
  108. {
  109. fHitCollection->Clear();
  110. fNHits = -1; //Must start in -1
  111. MpdSftPoint *pt;
  112. //FairMCTrack *mc;
  113. Int_t nTofPoint = fTofPoints->GetEntries();
  114. //Int_t nMCTracks = fMCTracks ->GetEntries();
  115. Int_t tof_tracks = 0, tof_tracks_vert = 0, tof_tracks_local = 0;
  116. cout << "-I- MpdSftHitProducer : " << nTofPoint
  117. << " points in Tof for this event" << endl;
  118. //Some numbers on TOF distributions
  119. // for(Int_t p=0;p<nMCTracks;p++) {
  120. // mc = (FairMCTrack*) fMCTracks->At(p);
  121. // if(mc->GetNPoints(kTOF)>0) tof_tracks++;
  122. // if(mc->GetNPoints(kTOF)>0 && mc->GetStartZ()>996) tof_tracks_local++;
  123. // if(mc->GetNPoints(kTOF)>0 && mc->GetMotherID()==-1) tof_tracks_vert++;
  124. // }
  125. cout << "-I- MpdSftHitProducer : " << tof_tracks
  126. << " tracks in Tof " << endl;
  127. cout << "-I- MpdSftHitProducer : " << tof_tracks_vert
  128. << " tracks in Tof from vertex" << endl;
  129. cout << "-I- MpdSftHitProducer : " << tof_tracks_local
  130. << " tracks in Tof produced locally (Z>996 cm)" << endl;
  131. cout << "-I- MpdSftHitProducer : " << tof_tracks-tof_tracks_local
  132. << " tracks in Tof able to produce a hit" << endl;
  133. TVector3 pos;
  134. Int_t nregions=10, nmodules=500, ncells=500, ngaps=8;
  135. //FIXME: these parameters must be provided externally
  136. Double_t xHit, yHit, zHit, tHit;
  137. Double_t xHitErr, yHitErr, zHitErr;
  138. Double_t tl_new, tr_new;
  139. Double_t Z=1000, Dz=2.04; //FIXME: Introduce also Dz and Z as (constant) parameters
  140. Float_t sigma_T=0.08, sigma_Y=0.7, sigma_t_gap, t_o, T_smearing = 0, sigma_el,
  141. vprop = 20, Pgap = 0.75;
  142. //time[ns], position[cm], velocity[cm/ns]
  143. //FIXME: these parameters must be provided externally
  144. Int_t trackID, region, module, cell, flag, ref;
  145. //Here check for the validity of the parameters
  146. if(fSigmaY>1) cout<<"IRREALISTIC SFT POSITION RESOLUTION!! (HitProducer may crash)"<<endl;
  147. if((fSigmaT<0.07 && fSigmaT>0)||fSigmaT>0.2) cout<<"IRREALISTIC SFT RESOLUTION!! (HitProducer may crash)"<<endl;
  148. //Parameterizations. Now they depend on the geometry/algorithm. FIXME
  149. if(fSigmaY!=0) sigma_el = sqrt(2.)*fSigmaY/vprop*1.3;
  150. else sigma_el = sqrt(2.)*sigma_Y/vprop*1.3;
  151. if(fSigmaT!=0) {
  152. sigma_t_gap = sqrt(pow(fSigmaT/(0.5*pow(ngaps,-0.361) + 0.40),2)-1.4/1.3*pow(sigma_el,2));
  153. t_o = 1.93*fSigmaT;
  154. }
  155. else {
  156. sigma_t_gap = sqrt(pow(sigma_T/(0.5*pow(ngaps,-0.361) + 0.40),2)-1.4/1.3*pow(sigma_el,2));
  157. t_o = 1.93*sigma_T;
  158. }
  159. //Initialization of cell times
  160. for(Int_t i=0;i<nregions;i++){
  161. for(Int_t j=0;j<nmodules;j++){
  162. for(Int_t k=0;k<ncells;k++){
  163. tl[i][j][k]= 1e+6;
  164. tr[i][j][k]= 1e+6;
  165. }
  166. }
  167. }
  168. // Loop over the TOF points
  169. for (Int_t j=0; j < nTofPoint; j++ ) {
  170. // Probability that the avalanche is detected
  171. if(gRandom->Uniform(1)>Pgap) continue;
  172. // Get a pointer to the TOF point
  173. pt = (MpdSftPoint*) fTofPoints->At(j);
  174. if(pt == NULL) {
  175. cout<<"Be careful: hole in the MpdEtofPoint TClonesArray!"<<endl;
  176. continue;
  177. }
  178. // Reject particles produced in the last 4 cm. Better job must be
  179. // done here. For example: it could better to go up to the parent
  180. // particle and get its trackID, then the secondary is
  181. // processed. FIXME.
  182. //mc = (FairMCTrack*) fMCTracks-> At(pt->GetTrackID());
  183. // if((mc->GetStartZ())>996) continue;
  184. //Get relevant information from the point
  185. trackID = pt->GetTrackID();
  186. cell = pt->GetCell();
  187. module = pt->GetModule();
  188. region = pt->GetRegion();
  189. pt->Position(pos);
  190. if(fVerbose > 2) {
  191. pt->Print(""); //FIXME
  192. cout << endl;
  193. }
  194. T_smearing = gRandom->Gaus(t_o, sigma_t_gap);
  195. tl_new = pt->GetTime() + T_smearing - pos.Y()/vprop + gRandom->Gaus(0,sigma_el);
  196. tr_new = pt->GetTime() + T_smearing + pos.Y()/vprop + gRandom->Gaus(0,sigma_el);
  197. flag = 1;
  198. //Increase the counter for the TofHit TClonesArray if the first time a hit is attached to this cell
  199. if(tl[region][module][cell]==1e+6 && tr[region][module][cell]==1e+6) fNHits++;
  200. else {
  201. if(trackID_left[region][module][cell]!=trackID || trackID_right[region][module][cell]!=trackID) {
  202. flag = 2;
  203. //Worsening of the time response in the presence of 2 (or more) hits is not possible to implement
  204. //within the present algorith (taking the fastest hit). FIXME
  205. }
  206. }
  207. //Take the fastest time from all the points/gaps in this cell
  208. if(tl_new<tl[region][module][cell]) {
  209. tl[region][module][cell] = tl_new;
  210. trackID_left[region][module][cell] = trackID;
  211. point_left[region][module][cell] = j;
  212. }
  213. if(tr_new<tr[region][module][cell]) {
  214. tr[region][module][cell] = tr_new;
  215. trackID_right[region][module][cell] = trackID;
  216. point_right[region][module][cell] = j;
  217. }
  218. //X and Y depend on the orientation of the cell. FIXME
  219. //Parameters of the Hit. Different for pad or strip read-out.
  220. if (type[region][module][cell]=="s") {
  221. yHit = (tr[region][module][cell] - tl[region][module][cell])/2*vprop;
  222. yHitErr = sigma_el/sqrt(2.)*vprop/1.3; //This must be better done.
  223. //1.3 comes from a relatin lines above. FIXME
  224. tHit = (tl[region][module][cell] + tr[region][module][cell])/2;
  225. }
  226. else if (type[region][module][cell]=="p") {
  227. yHit = Dy[region][module][cell]/2/10 + Y[region][module][cell]/10;
  228. yHitErr = Dy[region][module][cell]/10/sqrt(12.);
  229. tHit = tl[region][module][cell] + pos.Y()/vprop;
  230. }
  231. xHit = Dx[region][module][cell]/2/10 + X[region][module][cell]/10;
  232. xHitErr = Dx[region][module][cell]/10/sqrt(12.);
  233. zHit = Dz/2 + Z;
  234. zHitErr = Dz/sqrt(12.);
  235. //Reference to the point that contributes to the left side.
  236. ref = point_left[region][module][cell];
  237. //Reference to the point that contributes to the right side.
  238. if(type[region][module][cell]=="s" &&
  239. (trackID_left[region][module][cell]!=trackID_right[region][module][cell])) {
  240. flag = 10 + point_right[region][module][cell];
  241. //More than one TRACK contributed to this Hit: store the reference as flag + an offset of 10.
  242. }
  243. TVector3 hitPos(xHit, yHit, zHit);
  244. TVector3 hitPosErr(xHitErr, yHitErr, zHitErr);
  245. AddHit(pt->GetDetectorID(), hitPos, hitPosErr,
  246. j, tHit, flag);
  247. }
  248. cout << "-I- MpdSftHitProducer : " << fNHits
  249. << " hits in Tof created" << endl;
  250. }
  251. // ---- Add Hit to HitCollection --------------------------------------
  252. void MpdSftHitProducer::AddHit(Int_t detID, TVector3 &posHit, TVector3 &posHitErr,
  253. Int_t ref, Double_t tHit, Int_t flag)
  254. {
  255. new((*fHitCollection)[fNHits]) MpdSftHit(detID, posHit, posHitErr,
  256. ref, tHit, flag);
  257. if(fVerbose > 1) {
  258. MpdSftHit* tofHit = (MpdSftHit*) fHitCollection->At(fNHits);
  259. tofHit->Print();
  260. cout << endl;
  261. }
  262. }
  263. // ---- Finish --------------------------------------------------------
  264. void MpdSftHitProducer::Finish()
  265. {
  266. }
  267. // ---- SetSigmaT -----------------------------------------------------
  268. void MpdSftHitProducer::SetSigmaT(Double_t sigma)
  269. {
  270. fSigmaT = sigma;
  271. }
  272. // ---- SetSigmaXY -----------------------------------------------------
  273. void MpdSftHitProducer::SetSigmaXY(Double_t sigma)
  274. {
  275. fSigmaXY = sigma;
  276. }
  277. // ---- SetSigmaY -----------------------------------------------------
  278. void MpdSftHitProducer::SetSigmaY(Double_t sigma)
  279. {
  280. fSigmaY = sigma;
  281. }
  282. // ---- SetSigmaZ -----------------------------------------------------
  283. void MpdSftHitProducer::SetSigmaZ(Double_t sigma)
  284. {
  285. fSigmaZ = sigma;
  286. }
  287. // ---- GetSigmaXY -----------------------------------------------------
  288. Double_t MpdSftHitProducer::GetSigmaT()
  289. {
  290. return fSigmaT;
  291. }
  292. // ---- GetSigmaXY -----------------------------------------------------
  293. Double_t MpdSftHitProducer::GetSigmaXY()
  294. {
  295. return fSigmaXY;
  296. }
  297. // ---- GetSigmaY -----------------------------------------------------
  298. Double_t MpdSftHitProducer::GetSigmaY()
  299. {
  300. return fSigmaY;
  301. }
  302. // ---- GetSigmaZ -----------------------------------------------------
  303. Double_t MpdSftHitProducer::GetSigmaZ()
  304. {
  305. return fSigmaZ;
  306. }
  307. ClassImp(MpdSftHitProducer)