Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-2004, 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 : /* $Id$ */
17 :
18 : //_________________________________________________________________________
19 : //*-- Implementation version v2 of EMCAL Manager class; SHASHLYK version
20 : //*-- An object of this class does not produce digits
21 : //*-- It is the one to use if you do want to produce outputs in TREEH
22 : //*--
23 : //*-- Author : Alexei Pavlinov (WSU)
24 : // : Adapted for DCAL by M.L. Wang CCNU Wuhan & Subatech Oct-23-2012
25 : // This Class not stores information on all particles prior to EMCAL entry - in order to facilitate analysis.
26 : // This is done by setting fIShunt =2, and flagging all parents of particles entering the EMCAL.
27 :
28 : #include <cassert>
29 : // --- ROOT system ---
30 : #include <TBrowser.h>
31 : #include <TClonesArray.h>
32 : #include <TH2.h>
33 : #include <TParticle.h>
34 : #include <TROOT.h>
35 : #include <TVirtualMC.h>
36 :
37 : // --- Standard library ---
38 :
39 : // --- AliRoot header files ---
40 : #include "AliEMCALv2.h"
41 : #include "AliEMCALHit.h"
42 : #include "AliEMCALGeometry.h"
43 : #include "AliRun.h"
44 : #include "AliHeader.h"
45 : #include "AliMC.h"
46 : #include "AliStack.h"
47 : #include "AliTrackReference.h"
48 : // for TRD1 case only; May 31,2006
49 :
50 42 : ClassImp(AliEMCALv2)
51 :
52 : //______________________________________________________________________
53 : AliEMCALv2::AliEMCALv2()
54 12 : : AliEMCALv1()
55 60 : {
56 : // ctor
57 24 : }
58 :
59 : //______________________________________________________________________
60 : AliEMCALv2::AliEMCALv2(const char *name, const char *title,
61 : const Bool_t checkGeoAndRun)
62 1 : : AliEMCALv1(name,title,checkGeoAndRun)
63 3 : {
64 : // Standard Creator.
65 :
66 : //fHits= new TClonesArray("AliEMCALHit",1000); //Already done in ctor of v1
67 :
68 1 : gAlice->GetMCApp()->AddHitList(fHits);
69 :
70 1 : fNhits = 0;
71 1 : fIshunt = 2; // All hits are associated with particles entering the calorimeter
72 1 : fTimeCut = 30e-09;
73 :
74 : //fGeometry = GetGeometry();
75 : // pointer already initialized in AliEMCALv0 ctor, grandmother of this class
76 : // not sure why it was also invoked here, comment and leave the comment.
77 :
78 2 : AliInfo("Very good, V2 version is used!");
79 2 : }
80 :
81 : //______________________________________________________________________
82 52 : AliEMCALv2::~AliEMCALv2(){
83 : // dtor
84 :
85 : //Already done in dtor of v1
86 : // if ( fHits ) {
87 : // fHits->Clear();
88 : // delete fHits;
89 : // fHits = 0;
90 : // }
91 52 : }
92 :
93 : //______________________________________________________________________
94 : void AliEMCALv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t iparent, Float_t ienergy,
95 : Int_t id, Float_t * hits,Float_t * p){
96 : // Add a hit to the hit list.
97 : // An EMCAL hit is the sum of all hits in a tower section
98 : // originating from the same entering particle
99 : static Int_t hitCounter;
100 : static AliEMCALHit *newHit, *curHit;
101 : static Bool_t deja ;
102 :
103 129686 : deja = kFALSE;
104 :
105 129686 : newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
106 3600188 : for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- ) {
107 1147287 : curHit = (AliEMCALHit*) (*fHits)[hitCounter];
108 : // We add hits with the same tracknumber, while GEANT treats
109 : // primaries succesively
110 1147287 : if(curHit->GetPrimary() != primary)
111 : break;
112 1147287 : if( *curHit == *newHit ) {
113 129102 : *curHit = *curHit + *newHit;
114 64551 : deja = kTRUE;
115 : // break; // 30-aug-04 by PAI
116 64551 : } // end if
117 : } // end for hitCounter
118 :
119 64843 : if ( !deja ) {
120 292 : new((*fHits)[fNhits]) AliEMCALHit(*newHit);
121 292 : fNhits++;
122 292 : }
123 : // printf(" fNhits %i \n", fNhits);
124 129686 : delete newHit;
125 64843 : }
126 :
127 : //______________________________________________________________________
128 : void AliEMCALv2::StepManager(void){
129 : // Accumulates hits as long as the track stays in a tower
130 :
131 : // position wrt MRS and energy deposited
132 : static Float_t xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
133 : static Float_t pmom[4]={0.,0.,0.,0.};
134 2432575 : static TLorentzVector pos; // Lorentz vector of the track current position.
135 1216289 : static TLorentzVector mom; // Lorentz vector of the track current momentum.
136 : static Float_t ienergy = 0; // part->Energy();
137 1216289 : static TString curVolName="";
138 : static int supModuleNumber=-1, moduleNumber=-1, yNumber=-1, xNumber=-1, absid=-1;
139 : static int keyGeom=1; //real TRD1 geometry
140 : static const char *vn = "SCMX"; // Apr 13, 2006 - only TRD1 case now
141 : static Float_t depositedEnergy=0.0;
142 :
143 1216286 : if(keyGeom == 0) {
144 0 : keyGeom = 2;
145 0 : if(TVirtualMC::GetMC()->VolId("PBMO")==0 || TVirtualMC::GetMC()->VolId("WSUC")==1) {
146 0 : vn = "SCMX"; // old TRD2(TRD1) or WSUC
147 0 : keyGeom = 1;
148 0 : }
149 0 : printf("AliEMCALv2::StepManager(): keyGeom %i : Sensetive volume %s \n",
150 0 : keyGeom, vn);
151 0 : if(TVirtualMC::GetMC()->VolId("WSUC")==1) printf(" WSUC - cosmic ray stand geometry \n");
152 : }
153 1216286 : Int_t tracknumber = gAlice->GetMCApp()->GetCurrentTrackNumber();
154 : Int_t parent=0;
155 : TParticle* part=0;
156 :
157 1216286 : curVolName = TVirtualMC::GetMC()->CurrentVolName();
158 2131997 : if(curVolName.Contains(vn) || curVolName.Contains("SCX")) { // We are in a scintillator layer; SCX for 3X3
159 :
160 368460 : if( ((depositedEnergy = TVirtualMC::GetMC()->Edep()) > 0.) && (TVirtualMC::GetMC()->TrackTime() < fTimeCut)){// Track is inside a scintillator and deposits some energy
161 : // Info("StepManager "," entry %i DE %f",++ientry, depositedEnergy); // for testing
162 64843 : if (fCurPrimary==-1)
163 23 : fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber);
164 :
165 129663 : if (fCurParent==-1 || tracknumber != fCurTrack) {
166 : // Check parentage
167 : parent=tracknumber;
168 :
169 7021 : if (fCurParent != -1) {
170 100480 : while (parent != fCurParent && parent != -1) {
171 : //TParticle *part=gAlice->GetMCApp()->Particle(parent);
172 46741 : part=gAlice->GetMCApp()->Particle(parent);
173 46741 : parent=part->GetFirstMother();
174 : }
175 : }
176 7021 : if (fCurParent==-1 || parent==-1) {
177 : //Int_t parent=tracknumber;
178 : //TParticle *part=gAlice->GetMCApp()->Particle(parent);
179 : parent=tracknumber;
180 76 : part=gAlice->GetMCApp()->Particle(parent);
181 366 : while (parent != -1 && fGeometry->IsInEMCALOrDCAL(part->Vx(),part->Vy(),part->Vz())) {
182 69 : parent=part->GetFirstMother();
183 69 : if (parent!=-1)
184 69 : part=gAlice->GetMCApp()->Particle(parent);
185 : }
186 76 : fCurParent=parent;
187 76 : if (fCurParent==-1)
188 0 : Error("StepManager","Cannot find parent");
189 : else {
190 : //TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
191 76 : part=gAlice->GetMCApp()->Particle(fCurParent);
192 76 : ienergy = part->Energy();
193 :
194 : //Add reference to parent in TR tree.
195 76 : AddTrackReference(tracknumber, AliTrackReference::kEMCAL);
196 :
197 : }
198 638 : while (parent != -1) {
199 281 : part=gAlice->GetMCApp()->Particle(parent);
200 281 : part->SetBit(kKeepBit);
201 281 : parent=part->GetFirstMother();
202 : }
203 : }
204 7021 : fCurTrack=tracknumber;
205 7021 : }
206 64843 : TVirtualMC::GetMC()->TrackPosition(pos);
207 64843 : xyzte[0] = pos[0];
208 64843 : xyzte[1] = pos[1];
209 64843 : xyzte[2] = pos[2];
210 64843 : xyzte[3] = TVirtualMC::GetMC()->TrackTime() ;
211 :
212 64843 : TVirtualMC::GetMC()->TrackMomentum(mom);
213 64843 : pmom[0] = mom[0];
214 64843 : pmom[1] = mom[1];
215 64843 : pmom[2] = mom[2];
216 64843 : pmom[3] = mom[3];
217 :
218 : // if(ientry%200 > 0) return; // testing
219 64843 : supModuleNumber = moduleNumber = yNumber = xNumber = absid = 0;
220 129686 : if(keyGeom >= 1) { // TRD1 case now
221 129686 : TVirtualMC::GetMC()->CurrentVolOffID(4, supModuleNumber);
222 64843 : TVirtualMC::GetMC()->CurrentVolOffID(3, moduleNumber);
223 64843 : TVirtualMC::GetMC()->CurrentVolOffID(1, yNumber);
224 64843 : TVirtualMC::GetMC()->CurrentVolOffID(0, xNumber); // really x number now
225 : Int_t CurrentSMType = 0;
226 128707 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"SMOD")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_Standard ;
227 979 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"SM10")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_Half ;
228 1716 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"SM3rd")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_3rd ;
229 484 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"DCSM")==0) CurrentSMType = AliEMCALGeometry::kDCAL_Standard ;
230 0 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"DCEXT")==0) CurrentSMType = AliEMCALGeometry::kDCAL_Ext ;
231 0 : else AliError("Unkown SM Type!!");
232 :
233 : Int_t preSM = 0;
234 150234 : while( fGeometry->GetSMType(preSM) != CurrentSMType ) preSM++;
235 64843 : supModuleNumber += preSM;
236 : // Nov 10,2006
237 64843 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(0),vn) != 0) { // 3X3 case
238 0 : if (strcmp(TVirtualMC::GetMC()->CurrentVolOffName(0),"SCX1")==0) xNumber=1;
239 0 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(0),"SCX2")==0) xNumber=2;
240 0 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(0),"SCX3")==0) xNumber=3;
241 0 : else Fatal("StepManager()", "Wrong name of sensitive volume in 3X3 case : %s ", TVirtualMC::GetMC()->CurrentVolOffName(0));
242 : }
243 64843 : } else {
244 0 : TVirtualMC::GetMC()->CurrentVolOffID(5, supModuleNumber);
245 0 : TVirtualMC::GetMC()->CurrentVolOffID(4, moduleNumber);
246 0 : TVirtualMC::GetMC()->CurrentVolOffID(1, yNumber);
247 0 : TVirtualMC::GetMC()->CurrentVolOffID(0, xNumber);
248 0 : if (strcmp(TVirtualMC::GetMC()->CurrentVolOffName(5),"SMOP")==0) supModuleNumber = 2*(supModuleNumber-1)+1;
249 0 : else if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(5),"SMON")==0) supModuleNumber = 2*(supModuleNumber-1)+2;
250 0 : else assert(0); // something wrong
251 : }
252 :
253 : // Due to problem with index ordering conventions the calcultation of absid is no more like this:
254 : //absid = fGeometry->GetAbsCellId(smNumber, moduleNumber-1, yNumber-1, xNumber-1);
255 :
256 : //Swap A side in Phi and C side in Eta due to wrong indexing.
257 64843 : Int_t iphi = -1;
258 64843 : Int_t ieta = -1;
259 64843 : Int_t smNumber = supModuleNumber-1;
260 : Int_t smType = 1;
261 64843 : fGeometry->GetCellPhiEtaIndexInSModule(smNumber,moduleNumber-1,yNumber-1,xNumber-1, iphi, ieta);
262 129686 : if (smNumber%2 == 0) {
263 67996 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"DCSM")==0) smType = 3; //DCal supermodule. previous design/idea
264 : else smType = 2;
265 2911 : ieta = ((fGeometry->GetCentersOfCellsEtaDir()).GetSize()* 2/smType -1)-ieta;// 47/31-ieta, revert the ordering on A side in order to keep convention.
266 2911 : } else {
267 61932 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"SM10")==0) smType = 2 ; //half supermodule. previous design/idea
268 62661 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"SM3rd")==0) smType = 3 ; //one third (installed in 2012) supermodule
269 61932 : if(strcmp(TVirtualMC::GetMC()->CurrentVolOffName(4),"DCEXT")==0) smType = 3 ; //one third (installed in 2012) supermodule
270 61932 : iphi= ((fGeometry->GetCentersOfCellsPhiDir()).GetSize()/smType-1)-iphi;// 23/7-iphi, revert the ordering on C side in order to keep convention.
271 : }
272 :
273 : //Once we know the indexes, calculate the absolute ID
274 64843 : absid = fGeometry->GetAbsCellIdFromCellIndexes(smNumber, iphi, ieta);
275 :
276 64843 : if (absid < 0) {
277 0 : printf(" supModuleNumber %i : moduleNumber %i : yNumber %i : xNumber %i \n",
278 0 : supModuleNumber, moduleNumber, yNumber, xNumber);
279 0 : Fatal("StepManager()", "Wrong id : %i ", absid) ;
280 0 : }
281 :
282 64843 : Float_t lightYield = depositedEnergy ;
283 : // Apply Birk's law (copied from G3BIRK)
284 :
285 64843 : if (TVirtualMC::GetMC()->TrackCharge()!=0) { // Check
286 : Float_t birkC1Mod = 0;
287 64370 : if (fBirkC0==1){ // Apply correction for higher charge states
288 64370 : if (TMath::Abs(TVirtualMC::GetMC()->TrackCharge())>=2) birkC1Mod = fBirkC1*7.2/12.6;
289 64370 : else birkC1Mod = fBirkC1;
290 : }
291 :
292 : Float_t dedxcm=0.;
293 128724 : if (TVirtualMC::GetMC()->TrackStep()>0) dedxcm=1000.*TVirtualMC::GetMC()->Edep()/TVirtualMC::GetMC()->TrackStep();
294 : else dedxcm=0;
295 64370 : lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm);
296 64370 : }
297 :
298 : // use sampling fraction to get original energy --HG
299 : // xyzte[4] = lightYield * fGeometry->GetSampling();
300 64843 : xyzte[4] = lightYield; // 15-dec-04
301 :
302 64843 : if (gDebug == -2)
303 0 : printf("#sm %2i #m %3i #x %1i #z %1i -> absid %i : xyzte[4] = %f\n",
304 0 : supModuleNumber,moduleNumber,yNumber,xNumber,absid, xyzte[4]);
305 :
306 64843 : AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid, xyzte, pmom);
307 64843 : } // there is deposited energy
308 : }
309 1216286 : }
|
