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 PAM_VMC=gSystem->Getenv("PAM_VMC");
    TString PLATFORM=gSystem->Getenv("PLATFORM");
    TFile * resspe = 
      new TFile(PAM_VMC+"/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: " <<PAM_VMC+"/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, Bool_t is_prim){
    ffi=f;
    //g->SetMaxStep(0.001); // 10 um
    switch(f) {
    case ENTERING:
      CleanHit();
       FillVolID();
      InitHit();
      //cout<<">>>ENTERING in PLANE: "<<GetPlaneID()<<endl;
    default: 
      UpdateHit(g, is_prim);
      //TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(gMC->TrackPid());
      //Double_t mass = particlePDG->Mass();
      //Double_t x[3];
      //Double_t B[3];
      //Double_t * xp = &x[0];
      //Double_t *Bp = &B[0];
      //g->TrackPosition(x[0],x[1],x[2]);
      //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;
      //PamVMCFieldMgr::Instance()->Field(xp,Bp);
      //cout<<"X,Y,Z: "<<x[0]<<" "<<x[1]<<" "<<x[2]<<" "<<"Bx, By, Bz: "<<B[0]<<" "<<B[1]<<" "<<B[2]<<endl;
      //fhit->Print();
      //cout<<"Edep"<<g->Edep()*1e6<<" 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
      //cout<<"PDG:"<<fhit->GetPDG()<<"TOF="<<fhit->GetTOF()<<endl;
      
      if(fhit->GetEREL() ){
        //cout<<"Saving.."<<endl;
        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


#ifndef PAMVMCTPANSD_H
#define PAMVMCTPANSD_H

class PamVMCTPANSD: public PamVMCDetectorSD{

public:
 PamVMCTPANSD():PamVMCDetectorSD("PamVMCDetectorHit","TPAN",1000)
    {
    }; 
  virtual void FillHit(fin f,TVirtualMC *g, Bool_t is_prim){
    ffi=f;
    switch(f) {
    case ENTERING:
      CleanHit();
      FillVolID();
      InitHit();         
    default: 
      UpdateHit(g, is_prim);
      break;      
    }
    
    switch(f){
    case EXITING:
      // Save hit if PATH, not energy release is greater than zero
      if(fhit->GetPATH()){
        SaveHit(fdname.Data());
      }
      break;
    default:
      break;
    }
  }      
  ClassDef(PamVMCTPANSD,1)
    
    };

#endif // PAMVMCTPAN_H


#ifndef PAMVMCTRCNSD_H
#define PAMVMCTRCNSD_H

class PamVMCTRCNSD: public PamVMCDetectorSD{

public:
 PamVMCTRCNSD():PamVMCDetectorSD("PamVMCDetectorHit","TRCN",1000)
    {
      
    };
  ClassDef(PamVMCTRCNSD,1)
    
    };

#endif // PAMVMCTRCN_H


#ifndef PAMVMCTRSLSD_H
#define PAMVMCTRSLSD_H

class PamVMCTRSLSD: public PamVMCDetectorSD{

public:
 PamVMCTRSLSD():PamVMCDetectorSD("PamVMCDetectorHit","TRSL",1000)
    {
      
    };
  ClassDef(PamVMCTRSLSD,1)
    
    };

#endif // PAMVMCTRSL_H
1	nikolas	1.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	pam-rm2	1.5	TString PAM_VMC=gSystem->Getenv("PAM_VMC");
70	nikolas	1.1	TString PLATFORM=gSystem->Getenv("PLATFORM");
71			TFile * resspe =
72	pam-rm2	1.5	new TFile(PAM_VMC+"/lib/tgt_"+PLATFORM+"/resxy_new.root");
73	nikolas	1.1	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	pam-rm2	1.5	<<" is not present in: " <<PAM_VMC+"/lib/tgt_"+PLATFORM
102	nikolas	1.1	<<" 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	pam-rm2	1.5	void FillHit(fin f,TVirtualMC *g, Bool_t is_prim){
152	nikolas	1.1	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	pam-rm2	1.5	UpdateHit(g, is_prim);
162	nikolas	1.1	//TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(gMC->TrackPid());
163			//Double_t mass = particlePDG->Mass();
164	pam-rm2	1.5	//Double_t x[3];
165			//Double_t B[3];
166			//Double_t * xp = &x[0];
167			//Double_t *Bp = &B[0];
168			//g->TrackPosition(x[0],x[1],x[2]);
169	nikolas	1.1	//g->TrackMomentum(x,y,z,e);
170			//Double_t P0 = TMath::Sqrt(xx+yy+z*z);
171			//Double_t energy = TMath::Sqrt(P0P0+massmass)-mass;
172			//cout<<"stepping... Stepsize"<<gMC->TrackStep()10000<<"um "<<" particle PDG is:"<<gMC->TrackPid()<<" Kin.energy:"<<energy1000000<<" keV"<<endl;
173	pam-rm2	1.5	//PamVMCFieldMgr::Instance()->Field(xp,Bp);
174			//cout<<"X,Y,Z: "<<x[0]<<" "<<x[1]<<" "<<x[2]<<" "<<"Bx, By, Bz: "<<B[0]<<" "<<B[1]<<" "<<B[2]<<endl;
175			//fhit->Print();
176			//cout<<"Edep"<<g->Edep()*1e6<<" kev"<<endl;
177	nikolas	1.1	if (g->Edep()) PreDigit();
178
179			break;
180			}
181
182			switch(f){
183			case EXITING:
184			//cout<<"<<<EXITING out"<<gMC->TrackPid()<<endl;
185			// Save hit if energy release is greater than zero
186	pam-rm2	1.5	//cout<<"PDG:"<<fhit->GetPDG()<<"TOF="<<fhit->GetTOF()<<endl;
187
188			if(fhit->GetEREL() ){
189			//cout<<"Saving.."<<endl;
190	nikolas	1.1	SaveTrkHit(fdname.Data());
191			}
192			break;
193			default:
194			break;
195			}
196			}
197
198
199
200			/*This method for simulation spatial tracker resolution, use data from
201			resxy_new.root file. Calls on track exiting when we know XOUT,YOUT,ZOUT */
202			void SetTrkHit(PamVMCTrkHit *t){
203			float xin = (float)fhit->GetXIN();
204			float yin = (float)fhit->GetYIN();
205			float zin = (float)fhit->GetZIN();
206			float xout = (float)fhit->GetXOUT();
207			float yout = (float)fhit->GetYOUT();
208			float zout = (float)fhit->GetZOUT();
209
210			gpdcspe_(&xin, &yin, &zin, &xout, &yout, &zout);
211
212			fGPCSPE=&gpcspe_;
213
214			// cout<<"GPCSPE"<<fGPCSPE->xave<<" "<<fGPCSPE->yave<<" "<<
215			// fGPCSPE->zave<<" " <<fGPCSPE->nxmult<<" "<<fGPCSPE->nymult<<endl;
216
217			fGPCSPE->SetHit(t);
218			}
219
220
221			/*This method designed for calculaton energy release on trakcer strips
222			is calls for each step*/
223
224			virtual void PreDigit(){
225
226			Double_t x,y,z,px,py,pz,e,p;
227
228			gMC->TrackMomentum(px,py,pz,e);
229
230
231			p=sqrt(pxpx+pypy+pz*pz);
232
233			if(p){ //sometimes FLUKA vives u p=0
234
235			gGeoManager->cd(gMC->CurrentVolPath());
236
237			int iact=ffi;
238
239			float *trpar = new float[6];
240
241			gMC->TrackPosition(x,y,z);
242
243			trpar[0]=(float)x;
244			trpar[1]=(float)y;
245			trpar[2]=(float)z;
246
247			trpar[3]=(float)px/p;
248			trpar[4]=(float)py/p;
249			trpar[5]=(float)pz/p;
250
251			int *numvol = new int[2];
252			numvol[0]=GetPlaneID();
253			//if (numvol[0]==6) numvol[0]=0;
254			numvol[1]=GetPadID();
255
256			float deloss=(float)gMC->Edep();
257
258			float step=(float)gMC->TrackStep();
259
260			int itypar=TDatabasePDG::Instance()->ConvertPdgToGeant3(gMC->TrackPid());
261
262			gpudiffusion_(&iact, trpar, numvol, &deloss, &step, &itypar);
263
264			delete numvol;
265			delete trpar;
266			}
267
268			}
269
270
271			/*This method calculate ADC amplitudes from strips. It use data arrays
272			PROXTANTI, PROYTANTI, GLOBSTRIPX, GLOBSTRIPY in common. Calls after
273			each event from application */
274
275			virtual void Digitize(){
276
277			gpdspe_();
278
279			fGPSSPE = &gpsspe_;
280
281			fgps->SetData(fGPSSPE);
282			}
283
284
285
286			Int_t GetPlaneID(){
287			Int_t pl=Int_t((fhit->GetPOS()-1)/6.)+1;
288			return pl;
289			}
290
291			Int_t GetPadID(){
292			return (fhit->GetPOS()-(GetPlaneID()-1)*6);
293			}
294
295			ClassDef(PamVMCTrkSD,1)
296
297			};
298
299	pam-rm2	1.5	#endif // PAMVMCTRKSD_H
300
301
302
303			#ifndef PAMVMCTPANSD_H
304			#define PAMVMCTPANSD_H
305
306			class PamVMCTPANSD: public PamVMCDetectorSD{
307
308			public:
309			PamVMCTPANSD():PamVMCDetectorSD("PamVMCDetectorHit","TPAN",1000)
310			{
311			};
312			virtual void FillHit(fin f,TVirtualMC *g, Bool_t is_prim){
313			ffi=f;
314			switch(f) {
315			case ENTERING:
316			CleanHit();
317			FillVolID();
318			InitHit();
319			default:
320			UpdateHit(g, is_prim);
321			break;
322			}
323
324			switch(f){
325			case EXITING:
326			// Save hit if PATH, not energy release is greater than zero
327			if(fhit->GetPATH()){
328			SaveHit(fdname.Data());
329			}
330			break;
331			default:
332			break;
333			}
334			}
335			ClassDef(PamVMCTPANSD,1)
336
337			};
338
339			#endif // PAMVMCTPAN_H
340
341
342			#ifndef PAMVMCTRCNSD_H
343			#define PAMVMCTRCNSD_H
344	nikolas	1.1
345	pam-rm2	1.5	class PamVMCTRCNSD: public PamVMCDetectorSD{
346	nikolas	1.1
347	pam-rm2	1.5	public:
348			PamVMCTRCNSD():PamVMCDetectorSD("PamVMCDetectorHit","TRCN",1000)
349			{
350
351			};
352			ClassDef(PamVMCTRCNSD,1)
353
354			};
355
356			#endif // PAMVMCTRCN_H
357
358
359			#ifndef PAMVMCTRSLSD_H
360			#define PAMVMCTRSLSD_H
361
362			class PamVMCTRSLSD: public PamVMCDetectorSD{
363
364			public:
365			PamVMCTRSLSD():PamVMCDetectorSD("PamVMCDetectorHit","TRSL",1000)
366			{
367
368			};
369			ClassDef(PamVMCTRSLSD,1)
370
371			};
372
373			#endif // PAMVMCTRSL_H