trk/include/PamVMCTrkSD.h

#ifndef PAMVMCTRKSD_H
#define PAMVMCTRKSD_H

#ifdef __cplusplus
extern "C" {
#endif

#define GPINIT gpinit_
#define GPUDIFFUSION gpudiffusion_
#define GPBEGIN gpbegin_
#define GPDSPE gpdspe_
#define GPDCSPE gpdcspe_

  extern void GPINIT(float*, float*);

  extern void GPBEGIN();

  extern void GPUDIFFUSION(int* ,float* ,int*, float* ,float* ,int*);

  extern void GPDSPE();

  extern void GPDCSPE(float*, float*, float*, float*, float*, float* );

#ifdef __cplusplus
}
#endif


#include <iostream>

#include "PamVMCDetectorSD.h"
#include "PamVMCTrkHit.h"
#include "PamVMCTrkF77GpsSpe.h"
#include "PamVMCTrkF77GpcSpe.h"
#include "PamVMCTrkDig.h"
#include <TDatabasePDG.h>
#include <TGeoManager.h>
#include <TSystem.h>
#include <TFile.h>
#include <TH1F.h>
 

class pGPSSPEData;


class PamVMCTrkSD: public PamVMCDetectorSD{
  
 private:
  cGPCSPE * fGPCSPE;
  cGPSSPE * fGPSSPE;
  pGPSSPEHits * fgps;

 public:
  PamVMCTrkSD():PamVMCDetectorSD("PamVMCTrkHit","TSPA",1000)
    {
      fGPCSPE = 0;
      fGPSSPE = 0;
      fgps = new pGPSSPEHits();
      Init();
    };
  
  virtual ~PamVMCTrkSD(){ delete fgps;} 

 //Loading SPE hit resolution data from resxy_new.root
 //which should present in the same directory then application library
 
  void Init(){ 
    TString G4WORKDIR=gSystem->Getenv("G4WORKDIR");
    TString PLATFORM=gSystem->Getenv("PLATFORM");
    TFile * resspe = 
      new TFile(G4WORKDIR+"/lib/tgt_"+PLATFORM+"/resxy_new.root");
    if(resspe->IsOpen()){
      TH1F * hsx;
      TH1F * hsy;
      TString num;
      float cwx[11][1000];
      float cwy[11][1000];
      for(Int_t j=0; j<11; j++){
        num=(::Form("%4i",1000+2*j));
        resspe->GetObject("h"+num+";1",hsx);
        num=(::Form("%4i",2000+2*j));
        resspe->GetObject("h"+num+";1",hsy);
        if ((hsx)&&(hsy)){
          for(Int_t i = 0; i<1000; i++){
            cwx[j][i]=hsx->GetBinContent(i);
            cwy[j][i]=hsy->GetBinContent(i);
          }
        }else{
          cout<<"!!!WARNING: SPE hit resolution data file corrupted!!!"<<endl
              <<"!!!Application terminated..."<<endl;
          // exit(-1);
        }
      }

       gpinit_((*cwx), (*cwy));
       
       delete hsx;
       delete hsy;
    } else cout<<"!!!WARNING: SPE hit resolution data file resxy_new.root"
             <<" is not present in: " <<G4WORKDIR+"/lib/tgt_"+PLATFORM
             <<" or corrupted!!!"<<endl;

    delete resspe;
  }

  /*This method for zeroing PROXTANTI, PROYTANTI, GLOBSTRIPX, GLOBSTRIPY
    arrays in fortran common. Calls before each event */
  virtual void PreEvent(){ gpbegin_(); }

  void Register(){
    
    PamRootManager::Instance()->
      Register(fdname.Data(),"TClonesArray", &fHitColl);
    TString a=fHitColl?"exist":"not exist"; 
    PamRootManager::Instance()->
      Register("GPSSPE","pGPSSPEHits",&fgps);
    
    if (fdig) ((PamVMCTrkDig*)fdig)->RegisterTrkData(fgps);

  }

   virtual void ClearHitColl(){
     fnohit=0;
     (*fHitColl).Clear("C");
     fgps->Clear();
   }


 virtual void Compress(){fgps->Compress(); (*fHitColl).Compress();}

 virtual void SaveTrkHit(const char * dname){
   PamVMCDetectorHit * t=(PamVMCDetectorHit*) CreateTrkHit(dname);
   *t=*fhit;
 }

 virtual PamVMCTrkHit * CreateTrkHit(const char * dname) 
  { 
    fdetID=(PamVMCGeoID*)fdmap(dname); 
    return CreateTrkHit();
 }

 virtual  PamVMCTrkHit * CreateTrkHit() 
  { 
    PamVMCTrkHit *hit=  (PamVMCTrkHit *) fHitColl->New(fnohit++);
    InitHit();
    SetTrkHit(hit);
    return hit;
  }

  void FillHit(fin f,TVirtualMC *g){
    ffi=f;
    //g->SetMaxStep(0.001); // 10 um
    switch(f) {
    case ENTERING:
      CleanHit();
       FillVolID();
      InitHit();
      //cout<<">>>ENTERING in PLANE: "<<GetPlaneID()<<endl;
    default: 
      UpdateHit(g);
      //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(gMC->TrackPid());
      //Double_t mass = particlePDG->Mass();
      //Double_t x,y,z,e;
      //g->TrackMomentum(x,y,z,e);
      //Double_t P0 = TMath::Sqrt(x*x+y*y+z*z);
      //Double_t energy =  TMath::Sqrt(P0*P0+mass*mass)-mass;
      //cout<<"stepping... Stepsize"<<gMC->TrackStep()*10000<<"um "<<" particle PDG is:"<<gMC->TrackPid()<<" Kin.energy:"<<energy*1000000<<" keV"<<endl;
     
      if (g->Edep()) PreDigit(); 
      
      break;      
    }

    switch(f){
    case EXITING:
      //cout<<"<<<EXITING out"<<gMC->TrackPid()<<endl;
    // Save hit if energy release is greater than zero
      if(fhit->GetEREL()){
        SaveTrkHit(fdname.Data());
      }
      break;
    default:
      break;
    }
  }
  

  /*This method for simulation spatial tracker resolution, use data from
    resxy_new.root file. Calls on track exiting when we know XOUT,YOUT,ZOUT */
  void SetTrkHit(PamVMCTrkHit *t){
    float xin = (float)fhit->GetXIN();
    float yin = (float)fhit->GetYIN();
    float zin = (float)fhit->GetZIN();
    float xout = (float)fhit->GetXOUT();
    float yout = (float)fhit->GetYOUT();
    float zout = (float)fhit->GetZOUT();

    gpdcspe_(&xin, &yin, &zin, &xout, &yout, &zout);

    fGPCSPE=&gpcspe_;

    //      cout<<"GPCSPE"<<fGPCSPE->xave<<" "<<fGPCSPE->yave<<" "<<
    //  fGPCSPE->zave<<" " <<fGPCSPE->nxmult<<" "<<fGPCSPE->nymult<<endl;

        fGPCSPE->SetHit(t);
  }


  /*This method designed for calculaton energy release on trakcer strips
    is calls for each step*/

  virtual void PreDigit(){
     
      Double_t x,y,z,px,py,pz,e,p;
      
      gMC->TrackMomentum(px,py,pz,e);
      
     
      p=sqrt(px*px+py*py+pz*pz);
     
      if(p){ //sometimes FLUKA vives u p=0

        gGeoManager->cd(gMC->CurrentVolPath());
      
        int iact=ffi;

        float *trpar = new float[6]; 

        gMC->TrackPosition(x,y,z);

        trpar[0]=(float)x;
        trpar[1]=(float)y;
        trpar[2]=(float)z;

        trpar[3]=(float)px/p;
        trpar[4]=(float)py/p;
        trpar[5]=(float)pz/p;
     
        int *numvol = new int[2];
        numvol[0]=GetPlaneID();
        //if (numvol[0]==6) numvol[0]=0;
        numvol[1]=GetPadID();

        float deloss=(float)gMC->Edep();
        
        float step=(float)gMC->TrackStep();
        
        int itypar=TDatabasePDG::Instance()->ConvertPdgToGeant3(gMC->TrackPid());
        
        gpudiffusion_(&iact, trpar, numvol, &deloss, &step, &itypar);
        
        delete numvol;
        delete trpar;
      }

  }


  /*This method calculate ADC amplitudes from strips. It use data arrays
    PROXTANTI, PROYTANTI, GLOBSTRIPX, GLOBSTRIPY in common. Calls after
    each event from application */

 virtual void Digitize(){

   gpdspe_();

   fGPSSPE = &gpsspe_; 

   fgps->SetData(fGPSSPE);
 }


Int_t GetPlaneID(){
  Int_t pl=Int_t((fhit->GetPOS()-1)/6.)+1;
  return pl;
}

Int_t GetPadID(){
    return (fhit->GetPOS()-(GetPlaneID()-1)*6);
}

  ClassDef(PamVMCTrkSD,1)

};


#endif // PAMVMCTRKSD_H
1	#ifndef PAMVMCTRKSD_H
2	#define PAMVMCTRKSD_H
3
4	#ifdef __cplusplus
5	extern "C" {
6	#endif
7
8	#define GPINIT gpinit_
9	#define GPUDIFFUSION gpudiffusion_
10	#define GPBEGIN gpbegin_
11	#define GPDSPE gpdspe_
12	#define GPDCSPE gpdcspe_
13
14	extern void GPINIT(float, float);
15
16	extern void GPBEGIN();
17
18	extern void GPUDIFFUSION(int* ,float* ,int, float ,float* ,int*);
19
20	extern void GPDSPE();
21
22	extern void GPDCSPE(float, float, float, float, float, float );
23
24	#ifdef __cplusplus
25	}
26	#endif
27
28
29
30	#include <iostream>
31
32	#include "PamVMCDetectorSD.h"
33	#include "PamVMCTrkHit.h"
34	#include "PamVMCTrkF77GpsSpe.h"
35	#include "PamVMCTrkF77GpcSpe.h"
36	#include "PamVMCTrkDig.h"
37	#include <TDatabasePDG.h>
38	#include <TGeoManager.h>
39	#include <TSystem.h>
40	#include <TFile.h>
41	#include <TH1F.h>
42
43
44	class pGPSSPEData;
45
46
47	class PamVMCTrkSD: public PamVMCDetectorSD{
48
49	private:
50	cGPCSPE * fGPCSPE;
51	cGPSSPE * fGPSSPE;
52	pGPSSPEHits * fgps;
53
54	public:
55	PamVMCTrkSD():PamVMCDetectorSD("PamVMCTrkHit","TSPA",1000)
56	{
57	fGPCSPE = 0;
58	fGPSSPE = 0;
59	fgps = new pGPSSPEHits();
60	Init();
61	};
62
63	virtual ~PamVMCTrkSD(){ delete fgps;}
64
65	//Loading SPE hit resolution data from resxy_new.root
66	//which should present in the same directory then application library
67
68	void Init(){
69	TString G4WORKDIR=gSystem->Getenv("G4WORKDIR");
70	TString PLATFORM=gSystem->Getenv("PLATFORM");
71	TFile * resspe =
72	new TFile(G4WORKDIR+"/lib/tgt_"+PLATFORM+"/resxy_new.root");
73	if(resspe->IsOpen()){
74	TH1F * hsx;
75	TH1F * hsy;
76	TString num;
77	float cwx[11][1000];
78	float cwy[11][1000];
79	for(Int_t j=0; j<11; j++){
80	num=(::Form("%4i",1000+2*j));
81	resspe->GetObject("h"+num+";1",hsx);
82	num=(::Form("%4i",2000+2*j));
83	resspe->GetObject("h"+num+";1",hsy);
84	if ((hsx)&&(hsy)){
85	for(Int_t i = 0; i<1000; i++){
86	cwx[j][i]=hsx->GetBinContent(i);
87	cwy[j][i]=hsy->GetBinContent(i);
88	}
89	}else{
90	cout<<"!!!WARNING: SPE hit resolution data file corrupted!!!"<<endl
91	<<"!!!Application terminated..."<<endl;
92	// exit(-1);
93	}
94	}
95
96	gpinit_((cwx), (cwy));
97
98	delete hsx;
99	delete hsy;
100	} else cout<<"!!!WARNING: SPE hit resolution data file resxy_new.root"
101	<<" is not present in: " <<G4WORKDIR+"/lib/tgt_"+PLATFORM
102	<<" or corrupted!!!"<<endl;
103
104	delete resspe;
105	}
106
107	/*This method for zeroing PROXTANTI, PROYTANTI, GLOBSTRIPX, GLOBSTRIPY
108	arrays in fortran common. Calls before each event */
109	virtual void PreEvent(){ gpbegin_(); }
110
111	void Register(){
112
113	PamRootManager::Instance()->
114	Register(fdname.Data(),"TClonesArray", &fHitColl);
115	TString a=fHitColl?"exist":"not exist";
116	PamRootManager::Instance()->
117	Register("GPSSPE","pGPSSPEHits",&fgps);
118
119	if (fdig) ((PamVMCTrkDig*)fdig)->RegisterTrkData(fgps);
120
121	}
122
123	virtual void ClearHitColl(){
124	fnohit=0;
125	(*fHitColl).Clear("C");
126	fgps->Clear();
127	}
128
129
130	virtual void Compress(){fgps->Compress(); (*fHitColl).Compress();}
131
132	virtual void SaveTrkHit(const char * dname){
133	PamVMCDetectorHit * t=(PamVMCDetectorHit*) CreateTrkHit(dname);
134	t=fhit;
135	}
136
137	virtual PamVMCTrkHit * CreateTrkHit(const char * dname)
138	{
139	fdetID=(PamVMCGeoID*)fdmap(dname);
140	return CreateTrkHit();
141	}
142
143	virtual PamVMCTrkHit * CreateTrkHit()
144	{
145	PamVMCTrkHit hit= (PamVMCTrkHit ) fHitColl->New(fnohit++);
146	InitHit();
147	SetTrkHit(hit);
148	return hit;
149	}
150
151	void FillHit(fin f,TVirtualMC *g){
152	ffi=f;
153	//g->SetMaxStep(0.001); // 10 um
154	switch(f) {
155	case ENTERING:
156	CleanHit();
157	FillVolID();
158	InitHit();
159	//cout<<">>>ENTERING in PLANE: "<<GetPlaneID()<<endl;
160	default:
161	UpdateHit(g);
162	//TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(gMC->TrackPid());
163	//Double_t mass = particlePDG->Mass();
164	//Double_t x,y,z,e;
165	//g->TrackMomentum(x,y,z,e);
166	//Double_t P0 = TMath::Sqrt(xx+yy+z*z);
167	//Double_t energy = TMath::Sqrt(P0P0+massmass)-mass;
168	//cout<<"stepping... Stepsize"<<gMC->TrackStep()10000<<"um "<<" particle PDG is:"<<gMC->TrackPid()<<" Kin.energy:"<<energy1000000<<" keV"<<endl;
169
170	if (g->Edep()) PreDigit();
171
172	break;
173	}
174
175	switch(f){
176	case EXITING:
177	//cout<<"<<<EXITING out"<<gMC->TrackPid()<<endl;
178	// Save hit if energy release is greater than zero
179	if(fhit->GetEREL()){
180	SaveTrkHit(fdname.Data());
181	}
182	break;
183	default:
184	break;
185	}
186	}
187
188
189
190	/*This method for simulation spatial tracker resolution, use data from
191	resxy_new.root file. Calls on track exiting when we know XOUT,YOUT,ZOUT */
192	void SetTrkHit(PamVMCTrkHit *t){
193	float xin = (float)fhit->GetXIN();
194	float yin = (float)fhit->GetYIN();
195	float zin = (float)fhit->GetZIN();
196	float xout = (float)fhit->GetXOUT();
197	float yout = (float)fhit->GetYOUT();
198	float zout = (float)fhit->GetZOUT();
199
200	gpdcspe_(&xin, &yin, &zin, &xout, &yout, &zout);
201
202	fGPCSPE=&gpcspe_;
203
204	// cout<<"GPCSPE"<<fGPCSPE->xave<<" "<<fGPCSPE->yave<<" "<<
205	// fGPCSPE->zave<<" " <<fGPCSPE->nxmult<<" "<<fGPCSPE->nymult<<endl;
206
207	fGPCSPE->SetHit(t);
208	}
209
210
211	/*This method designed for calculaton energy release on trakcer strips
212	is calls for each step*/
213
214	virtual void PreDigit(){
215
216	Double_t x,y,z,px,py,pz,e,p;
217
218	gMC->TrackMomentum(px,py,pz,e);
219
220
221	p=sqrt(pxpx+pypy+pz*pz);
222
223	if(p){ //sometimes FLUKA vives u p=0
224
225	gGeoManager->cd(gMC->CurrentVolPath());
226
227	int iact=ffi;
228
229	float *trpar = new float[6];
230
231	gMC->TrackPosition(x,y,z);
232
233	trpar[0]=(float)x;
234	trpar[1]=(float)y;
235	trpar[2]=(float)z;
236
237	trpar[3]=(float)px/p;
238	trpar[4]=(float)py/p;
239	trpar[5]=(float)pz/p;
240
241	int *numvol = new int[2];
242	numvol[0]=GetPlaneID();
243	//if (numvol[0]==6) numvol[0]=0;
244	numvol[1]=GetPadID();
245
246	float deloss=(float)gMC->Edep();
247
248	float step=(float)gMC->TrackStep();
249
250	int itypar=TDatabasePDG::Instance()->ConvertPdgToGeant3(gMC->TrackPid());
251
252	gpudiffusion_(&iact, trpar, numvol, &deloss, &step, &itypar);
253
254	delete numvol;
255	delete trpar;
256	}
257
258	}
259
260
261	/*This method calculate ADC amplitudes from strips. It use data arrays
262	PROXTANTI, PROYTANTI, GLOBSTRIPX, GLOBSTRIPY in common. Calls after
263	each event from application */
264
265	virtual void Digitize(){
266
267	gpdspe_();
268
269	fGPSSPE = &gpsspe_;
270
271	fgps->SetData(fGPSSPE);
272	}
273
274
275
276	Int_t GetPlaneID(){
277	Int_t pl=Int_t((fhit->GetPOS()-1)/6.)+1;
278	return pl;
279	}
280
281	Int_t GetPadID(){
282	return (fhit->GetPOS()-(GetPlaneID()-1)*6);
283	}
284
285	ClassDef(PamVMCTrkSD,1)
286
287	};
288
289
290
291	#endif // PAMVMCTRKSD_H