Line data Source code
1 : // **************************************************************************
2 : // * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : // * *
4 : // * Author: The ALICE Off-line Project. *
5 : // * Contributors are mentioned in the code where appropriate. *
6 : // * *
7 : // * Permission to use, copy, modify and distribute this software and its *
8 : // * documentation strictly for non-commercial purposes is hereby granted *
9 : // * without fee, provided that the above copyright notice appears in all *
10 : // * copies and that both the copyright notice and this permission notice *
11 : // * appear in the supporting documentation. The authors make no claims *
12 : // * about the suitability of this software for any purpose. It is *
13 : // * provided "as is" without express or implied warranty. *
14 : // **************************************************************************
15 :
16 :
17 : #include "AliHMPIDvG4.h" // class header
18 : #include "AliHMPIDParam.h" // StepManager()
19 : #include "AliHMPIDHit.h" // Hits2SDigs(),StepManager()
20 : #include "AliHMPIDDigit.h" // Digits2Raw(), Raw2SDigits()
21 : #include "AliHMPIDRawStream.h" // Digits2Raw(), Raw2SDigits()
22 : #include "AliRawReader.h" // Raw2SDigits()
23 : #include "AliTrackReference.h" //
24 : #include <TVirtualMC.h> // StepManager() for TVirtualMC::GetMC()
25 : #include <TPDGCode.h> // StepHistory()
26 : #include <AliStack.h> // StepManager(),Hits2SDigits()78.6
27 : #include <AliLoader.h> // Hits2SDigits()
28 : #include <AliRunLoader.h> // Hits2SDigits()
29 : #include <AliMC.h> // StepManager()
30 : #include <AliRun.h> // CreateMaterials()
31 : #include <AliMagF.h> // CreateMaterials()
32 : #include "AliGeomManager.h" // AddAlignableVolumes()
33 : #include <AliCDBEntry.h> // CreateMaterials()
34 : #include <AliCDBManager.h> // CreateMaterials()
35 : #include <TF1.h> // DefineOpticalProperties()
36 : #include <TF2.h> // DefineOpticalProperties()
37 : #include <TGeoCompositeShape.h> // CradleBaseVolume()
38 : #include <TGeoGlobalMagField.h> //
39 : #include <TGeoPhysicalNode.h> // AddAlignableVolumes()
40 : #include <TGeoXtru.h> // CradleBaseVolume()
41 : #include <TLorentzVector.h> // IsLostByFresnel()
42 : #include <TString.h> // StepManager()
43 : #include <TTree.h> //
44 :
45 12 : ClassImp(AliHMPIDvG4)
46 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 : void AliHMPIDvG4::GenFee(Float_t qtot)
48 : {
49 : // Generate FeedBack photons for the current particle. To be invoked from StepManager().
50 : // eloss=0 means photon so only pulse height distribution is to be analysed.
51 0 : TLorentzVector x4;
52 0 : TVirtualMC::GetMC()->TrackPosition(x4);
53 0 : Int_t iNphotons=TVirtualMC::GetMC()->GetRandom()->Poisson(0.02*qtot); //# of feedback photons is proportional to the charge of hit
54 0 : AliDebug(1,Form("N photons=%i",iNphotons));
55 0 : if (iNphotons > fMaxFeed) return;
56 : Int_t j;
57 0 : Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop,dir[3], phi,pol[3], mom[4];
58 : //Generate photons
59 0 : for(Int_t i=0;i<iNphotons;i++){//feedbacks loop
60 0 : Double_t ranf[2];
61 0 : TVirtualMC::GetMC()->GetRandom()->RndmArray(2,ranf); //Sample direction
62 0 : cthf=ranf[0]*2-1.0;
63 0 : if(cthf<0) continue;
64 0 : sthf = TMath::Sqrt((1. - cthf) * (1. + cthf));
65 0 : phif = ranf[1] * 2 * TMath::Pi();
66 :
67 0 : if(Double_t randomNumber=TVirtualMC::GetMC()->GetRandom()->Rndm()<=0.57)
68 0 : enfp = 7.5e-9;
69 0 : else if(randomNumber<=0.7)
70 0 : enfp = 6.4e-9;
71 : else
72 : enfp = 7.9e-9;
73 :
74 :
75 0 : dir[0] = sthf * TMath::Sin(phif); dir[1] = cthf; dir[2] = sthf * TMath::Cos(phif);
76 0 : TVirtualMC::GetMC()->Gdtom(dir, mom, 2);
77 0 : mom[0]*=enfp; mom[1]*=enfp; mom[2]*=enfp;
78 0 : mom[3] = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
79 :
80 : // Polarisation
81 0 : e1[0]= 0; e1[1]=-dir[2]; e1[2]= dir[1];
82 0 : e2[0]=-dir[1]; e2[1]= dir[0]; e2[2]= 0;
83 0 : e3[0]= dir[1]; e3[1]= 0; e3[2]=-dir[0];
84 :
85 : vmod=0;
86 0 : for(j=0;j<3;j++) vmod+=e1[j]*e1[j];
87 0 : if (!vmod) for(j=0;j<3;j++) {
88 0 : uswop=e1[j];
89 0 : e1[j]=e3[j];
90 0 : e3[j]=uswop;
91 : }
92 : vmod=0;
93 0 : for(j=0;j<3;j++) vmod+=e2[j]*e2[j];
94 0 : if (!vmod) for(j=0;j<3;j++) {
95 0 : uswop=e2[j];
96 0 : e2[j]=e3[j];
97 0 : e3[j]=uswop;
98 : }
99 :
100 0 : vmod=0; for(j=0;j<3;j++) vmod+=e1[j]*e1[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e1[j]*=vmod;
101 0 : vmod=0; for(j=0;j<3;j++) vmod+=e2[j]*e2[j]; vmod=TMath::Sqrt(1/vmod); for(j=0;j<3;j++) e2[j]*=vmod;
102 :
103 0 : phi = TVirtualMC::GetMC()->GetRandom()->Rndm()* 2 * TMath::Pi();
104 0 : for(j=0;j<3;j++) pol[j]=e1[j]*TMath::Sin(phi)+e2[j]*TMath::Cos(phi);
105 0 : TVirtualMC::GetMC()->Gdtom(pol, pol, 2);
106 0 : Int_t outputNtracksStored;
107 0 : gAlice->GetMCApp()->PushTrack(1, //transport
108 0 : gAlice->GetMCApp()->GetCurrentTrackNumber(),//parent track
109 : 50000051, //PID
110 0 : mom[0],mom[1],mom[2],mom[3], //track momentum
111 0 : x4.X(),x4.Y(),x4.Z(),x4.T(), //track origin
112 0 : pol[0],pol[1],pol[2], //polarization
113 : kPFeedBackPhoton, //process ID
114 : outputNtracksStored, //on return how many new photons stored on stack
115 : 1.0); //weight
116 0 : }//feedbacks loop
117 0 : AliDebug(1,"Stop.");
118 0 : }//GenerateFeedbacks()
119 :
120 : //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
121 :
122 : void AliHMPIDvG4::StepManager()
123 : {
124 : // Full Step Manager.
125 : // Arguments: none
126 : // Returns: none
127 : // StepHistory(); return; //uncomment to print tracks history
128 : // StepCount(); return; //uncomment to count photons
129 :
130 0 : TString volname = TVirtualMC::GetMC()->CurrentVolName();
131 :
132 : //Treat photons
133 0 : if((TVirtualMC::GetMC()->TrackPid()==50000050||TVirtualMC::GetMC()->TrackPid()==50000051)&&volname.Contains("Hpad")){ //photon (Ckov or feedback) hits on module PC (Hpad)
134 0 : if(TVirtualMC::GetMC()->Edep()>0){ //photon survided QE test i.e. produces electron
135 0 : if(IsLostByFresnel()){ TVirtualMC::GetMC()->StopTrack(); return;} //photon lost due to fersnel reflection on PC
136 0 : Int_t tid= TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber(); //take TID
137 0 : Int_t pid= TVirtualMC::GetMC()->TrackPid(); //take PID
138 0 : Float_t etot= TVirtualMC::GetMC()->Etot(); //total hpoton energy, [GeV]
139 0 : Double_t x[3]; TVirtualMC::GetMC()->TrackPosition(x[0],x[1],x[2]); //take MARS position at entrance to PC
140 0 : Float_t hitTime= (Float_t)TVirtualMC::GetMC()->TrackTime(); //hit formation time
141 0 : TString tmpname = volname; tmpname.Remove(0,4); Int_t idch = tmpname.Atoi(); //retrieve the chamber number
142 0 : Float_t xl,yl; AliHMPIDParam::Instance()->Mars2Lors(idch,x,xl,yl); //take LORS position
143 0 : new((*fHits)[fNhits++])AliHMPIDHit(idch,etot,pid,tid,xl,yl,hitTime,x); //HIT for photon, position at P, etot will be set to Q
144 0 : if(fDoFeed) {
145 0 : Int_t nfeed = etot * 0.02;
146 0 : if (nfeed > fMaxFeed) {
147 0 : printf("Nfeed: %5d eloss: %13.3f pid: %5d tid: %5d \n", nfeed, etot, pid, tid);
148 0 : StepHistory();
149 : }
150 0 : GenFee(etot); //generate feedback photons etot is modified in hit ctor to Q of hit
151 0 : }
152 0 : }//photon hit PC and DE >0
153 : }//photon hit PC
154 :
155 :
156 : //Treat charged particles
157 : static Float_t eloss; //need to store mip parameters between different steps
158 : static Double_t in[3];
159 :
160 0 : if(TVirtualMC::GetMC()->IsTrackEntering() && TVirtualMC::GetMC()->TrackCharge() && volname.Contains("Hpad")) //Trackref stored when entering in the pad volume
161 0 : AddTrackReference(TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber(), AliTrackReference::kHMPID); //for acceptance calculations
162 :
163 :
164 0 : if(TVirtualMC::GetMC()->TrackCharge() && volname.Contains("Hcel")){ //charged particle in amplification gap (Hcel)
165 0 : if(TVirtualMC::GetMC()->IsTrackEntering()||TVirtualMC::GetMC()->IsNewTrack()) { //entering or newly created
166 0 : eloss=0; //reset Eloss collector
167 0 : TVirtualMC::GetMC()->TrackPosition(in[0],in[1],in[2]); //take position at the entrance
168 0 : }else if(TVirtualMC::GetMC()->IsTrackExiting()||TVirtualMC::GetMC()->IsTrackStop()||TVirtualMC::GetMC()->IsTrackDisappeared()){ //exiting or disappeared
169 0 : eloss +=TVirtualMC::GetMC()->Edep(); //take into account last step Eloss
170 0 : Int_t tid= TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber(); //take TID
171 0 : Int_t pid= TVirtualMC::GetMC()->TrackPid(); //take PID
172 0 : Double_t out[3]; TVirtualMC::GetMC()->TrackPosition(out[0],out[1],out[2]); //take MARS position at exit
173 0 : Float_t hitTime= (Float_t)TVirtualMC::GetMC()->TrackTime(); //hit formation time
174 0 : out[0]=0.5*(out[0]+in[0]); //
175 0 : out[1]=0.5*(out[1]+in[1]); //take hit position at the anod plane
176 0 : out[2]=0.5*(out[2]+in[2]);
177 0 : TString tmpname = volname; tmpname.Remove(0,4); Int_t idch = tmpname.Atoi(); //retrieve the chamber number
178 0 : Float_t xl,yl;AliHMPIDParam::Instance()->Mars2Lors(idch,out,xl,yl); //take LORS position
179 0 : if(eloss>0) {
180 0 : new((*fHits)[fNhits++])AliHMPIDHit(idch,eloss,pid,tid,xl,yl,hitTime,out); //HIT for MIP, position near anod plane, eloss will be set to Q
181 0 : if(fDoFeed){
182 0 : Int_t nfeed = 0.02 * eloss;
183 0 : if (nfeed > fMaxFeed) {
184 0 : printf("Nfeed: %5d eloss: %13.3f pid: %5d tid: %5d \n", nfeed, eloss, pid, tid);
185 0 : StepHistory();
186 : }
187 0 : GenFee(eloss); //generate feedback photons
188 0 : }
189 0 : eloss=0;
190 0 : }
191 0 : }else //just going inside
192 0 : eloss += TVirtualMC::GetMC()->Edep(); //collect this step eloss
193 : }//MIP in GAP
194 :
195 0 : }//StepManager()
|
