trk/include/PamVMCTrkSD.h_old

#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 <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;
  }

  
  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);
  }

   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;
    switch(f) {
    case ENTERING:
      CleanHit();
       FillVolID();
      InitHit();         
    default: 
      UpdateHit(g);
     
      if (g->Edep()) PreDigit();
      
      break;      
    }

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

  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);
  }

  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){

        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;
      }

  }


 virtual void Digitize(){

   gpdspe_();

   fGPSSPE = &gpsspe_; 

   fgps->SetData(fGPSSPE);


 }


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

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

  ClassDef(PamVMCTrkSD,1)

};


#endif // PAMVMCTRKSD_H
1	pamelats	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 <TDatabasePDG.h>
37			#include <TGeoManager.h>
38			#include <TSystem.h>
39			#include <TFile.h>
40			#include <TH1F.h>
41
42
43			class pGPSSPEData;
44
45
46			class PamVMCTrkSD: public PamVMCDetectorSD{
47
48			private:
49			cGPCSPE * fGPCSPE;
50			cGPSSPE * fGPSSPE;
51			pGPSSPEHits * fgps;
52
53			public:
54			PamVMCTrkSD():PamVMCDetectorSD("PamVMCTrkHit","TSPA",1000)
55			{
56			fGPCSPE = 0;
57			fGPSSPE = 0;
58			fgps = new pGPSSPEHits();
59			Init();
60			};
61
62			virtual ~PamVMCTrkSD(){ delete fgps; }
63
64			//Loading SPE hit resolution data from resxy_new.root
65			//which should present in the same directory then application library
66
67			void Init(){
68			TString G4WORKDIR=gSystem->Getenv("G4WORKDIR");
69			TString PLATFORM=gSystem->Getenv("PLATFORM");
70			TFile * resspe =
71			new TFile(G4WORKDIR+"/lib/tgt_"+PLATFORM+"/resxy_new.root");
72			if(resspe->IsOpen()){
73			TH1F * hsx;
74			TH1F * hsy;
75			TString num;
76			float cwx[11][1000];
77			float cwy[11][1000];
78			for(Int_t j=0; j<11; j++){
79			num=(::Form("%4i",1000+2*j));
80			resspe->GetObject("h"+num+";1",hsx);
81			num=(::Form("%4i",2000+2*j));
82			resspe->GetObject("h"+num+";1",hsy);
83			if ((hsx)&&(hsy)){
84			for(Int_t i = 0; i<1000; i++){
85			cwx[j][i]=hsx->GetBinContent(i);
86			cwy[j][i]=hsy->GetBinContent(i);
87			}
88			}else{
89			cout<<"!!!WARNING: SPE hit resolution data file corrupted!!!"<<endl
90			<<"!!!Application terminated..."<<endl;
91			// exit(-1);
92			}
93			}
94
95			gpinit_((cwx), (cwy));
96
97			delete hsx;
98			delete hsy;
99			} else cout<<"!!!WARNING: SPE hit resolution data file resxy_new.root"
100			<<" is not present in: " <<G4WORKDIR+"/lib/tgt_"+PLATFORM
101			<<" or corrupted!!!"<<endl;
102
103			delete resspe;
104			}
105
106
107			virtual void PreEvent(){ gpbegin_(); }
108
109			void Register(){
110
111			PamRootManager::Instance()->
112			Register(fdname.Data(),"TClonesArray", &fHitColl);
113			TString a=fHitColl?"exist":"not exist";
114			PamRootManager::Instance()->
115			Register("GPSSPE","pGPSSPEHits",&fgps);
116			}
117
118			virtual void ClearHitColl(){
119			fnohit=0;
120			(*fHitColl).Clear("C");
121			fgps->Clear();
122			}
123
124
125			virtual void Compress(){fgps->Compress(); (*fHitColl).Compress();}
126
127			virtual void SaveTrkHit(const char * dname){
128			PamVMCDetectorHit * t=(PamVMCDetectorHit*) CreateTrkHit(dname);
129			t=fhit;
130			}
131
132			virtual PamVMCTrkHit * CreateTrkHit(const char * dname)
133			{
134			fdetID=(PamVMCGeoID*)fdmap(dname);
135			return CreateTrkHit();
136			}
137
138			virtual PamVMCTrkHit * CreateTrkHit()
139			{
140			PamVMCTrkHit hit= (PamVMCTrkHit ) fHitColl->New(fnohit++);
141			InitHit();
142			SetTrkHit(hit);
143			return hit;
144			}
145
146			void FillHit(fin f,TVirtualMC *g){
147			ffi=f;
148			switch(f) {
149			case ENTERING:
150			CleanHit();
151			FillVolID();
152			InitHit();
153			default:
154			UpdateHit(g);
155
156			if (g->Edep()) PreDigit();
157
158			break;
159			}
160
161			switch(f){
162			case EXITING:
163			// Save hit if energy release is greater than zero
164			if(fhit->GetEREL()){
165			SaveTrkHit(fdname.Data());
166			}
167			break;
168			default:
169			break;
170			}
171			}
172
173
174
175
176			void SetTrkHit(PamVMCTrkHit *t){
177			float xin = (float)fhit->GetXIN();
178			float yin = (float)fhit->GetYIN();
179			float zin = (float)fhit->GetZIN();
180			float xout = (float)fhit->GetXOUT();
181			float yout = (float)fhit->GetYOUT();
182			float zout = (float)fhit->GetZOUT();
183
184			gpdcspe_(&xin, &yin, &zin, &xout, &yout, &zout);
185
186			fGPCSPE=&gpcspe_;
187
188			// cout<<"GPCSPE"<<fGPCSPE->xave<<" "<<fGPCSPE->yave<<" "<<
189			// fGPCSPE->zave<<" " <<fGPCSPE->nxmult<<" "<<fGPCSPE->nymult<<endl;
190
191
192
193			fGPCSPE->SetHit(t);
194			}
195
196			virtual void PreDigit(){
197
198			Double_t x,y,z,px,py,pz,e,p;
199
200			gMC->TrackMomentum(px,py,pz,e);
201
202
203			p=sqrt(pxpx+pypy+pz*pz);
204
205			if(p){
206
207			gGeoManager->cd(gMC->CurrentVolPath());
208
209			int iact=ffi;
210
211			float *trpar = new float[6];
212
213			gMC->TrackPosition(x,y,z);
214
215			trpar[0]=(float)x;
216			trpar[1]=(float)y;
217			trpar[2]=(float)z;
218
219			trpar[3]=(float)px/p;
220			trpar[4]=(float)py/p;
221			trpar[5]=(float)pz/p;
222
223			int *numvol = new int[2];
224			numvol[0]=GetPlaneID();
225			if (numvol[0]==6) numvol[0]=0;
226			numvol[1]=GetPadID();
227
228			float deloss=(float)gMC->Edep();
229
230			float step=(float)gMC->TrackStep();
231
232			int itypar=TDatabasePDG::Instance()->ConvertPdgToGeant3(gMC->TrackPid());
233
234			gpudiffusion_(&iact, trpar, numvol, &deloss, &step, &itypar);
235
236			delete numvol;
237			delete trpar;
238			}
239
240			}
241
242
243
244
245			virtual void Digitize(){
246
247			gpdspe_();
248
249			fGPSSPE = &gpsspe_;
250
251			fgps->SetData(fGPSSPE);
252
253
254			}
255
256
257
258			Int_t GetPlaneID(){
259			Int_t pl=Int_t(fhit->GetPOS()/6.+1);
260			if (pl==7) pl=6;
261			return pl;
262			}
263
264			Int_t GetPadID(){
265			return (fhit->GetPOS()-(GetPlaneID()-1)*6);
266			}
267
268			ClassDef(PamVMCTrkSD,1)
269
270			};
271
272
273
274			#endif // PAMVMCTRKSD_H