AliRoot Core  3dc7879 (3dc7879)
AliEMCALv2.cxx
Go to the documentation of this file.
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 // --- Standard library ---
17 #include <cassert>
18 
19 // --- ROOT system ---
20 #include <TBrowser.h>
21 #include <TClonesArray.h>
22 #include <TH2.h>
23 #include <TParticle.h>
24 #include <TROOT.h>
25 #include <TVirtualMC.h>
26 
27 // --- AliRoot header files ---
28 #include "AliEMCALv2.h"
29 #include "AliEMCALHit.h"
30 #include "AliEMCALGeometry.h"
31 #include "AliRun.h"
32 #include "AliHeader.h"
33 #include "AliMC.h"
34 #include "AliStack.h"
35 #include "AliTrackReference.h"
36 
38 ClassImp(AliEMCALv2) ;
40 
43 //______________________________________________________________________
45 {}
46 
53 //______________________________________________________________________
54 AliEMCALv2::AliEMCALv2(const char *name, const char *title,
55  const Bool_t checkGeoAndRun)
56 : AliEMCALv1(name,title,checkGeoAndRun)
57 {
58  //fHits= new TClonesArray("AliEMCALHit",1000); //Already done in ctor of v1
59 
61 
62  fNhits = 0;
63  fIshunt = 2; // All hits are associated with particles entering the calorimeter
64  fTimeCut = 30e-09;
65 
66  //fGeometry = GetGeometry();
67  // pointer already initialized in AliEMCALv0 ctor, grandmother of this class
68  // not sure why it was also invoked here, comment and leave the comment.
69 
70  AliInfo("Very good, V2 version is used!");
71 }
72 
76 //______________________________________________________________________
78 {}
79 
95 //______________________________________________________________________
96 void AliEMCALv2::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber,
97  Int_t iparent, Float_t ienergy,
98  Int_t id, Float_t * hits,Float_t * p)
99 {
100  static Int_t hitCounter;
101  static AliEMCALHit *newHit, *curHit;
102  static Bool_t deja ;
103 
104  deja = kFALSE;
105 
106  newHit = new AliEMCALHit(shunt, primary, tracknumber, iparent, ienergy, id, hits, p);
107 
108  for ( hitCounter = fNhits-1; hitCounter >= 0 && !deja; hitCounter-- )
109  {
110  curHit = (AliEMCALHit*) (*fHits)[hitCounter];
111 
112  // We add hits with the same tracknumber, while GEANT treats
113  // primaries succesively
114  if(curHit->GetPrimary() != primary)
115  break;
116 
117  if( *curHit == *newHit )
118  {
119  *curHit = *curHit + *newHit;
120  deja = kTRUE;
121  // break; // 30-aug-04 by PAI
122  } // end if
123  } // end for hitCounter
124 
125  if ( !deja )
126  {
127  new((*fHits)[fNhits]) AliEMCALHit(*newHit);
128  fNhits++;
129  }
130 
131  // printf(" fNhits %i \n", fNhits);
132  delete newHit;
133 }
134 
137 //______________________________________________________________________
139 {
140  // position wrt MRS and energy deposited
141  static Float_t xyzte[5]={0.,0.,0.,0.,0.};// position wrt MRS, time and energy deposited
142  static Float_t pmom[4]={0.,0.,0.,0.};
143  static TLorentzVector pos; // Lorentz vector of the track current position.
144  static TLorentzVector mom; // Lorentz vector of the track current momentum.
145  static Float_t ienergy = 0; // part->Energy();
146  static TString curVolName="";
147  static int supModuleNumber=-1, moduleNumber=-1, yNumber=-1, xNumber=-1, absid=-1;
148  static int keyGeom=1; //real TRD1 geometry
149  static const char *vn = "SCMX"; // Apr 13, 2006 - only TRD1 case now
150  static Float_t depositedEnergy=0.0;
151 
152  if(keyGeom == 0)
153  {
154  keyGeom = 2;
155  if(fMC->VolId("PBMO")==0 || fMC->VolId("WSUC")==1)
156  {
157  vn = "SCMX"; // old TRD2(TRD1) or WSUC
158  keyGeom = 1;
159  }
160 
161  printf("AliEMCALv2::StepManager(): keyGeom %i : Sensetive volume %s \n",
162  keyGeom, vn);
163 
164  if(fMC->VolId("WSUC")==1) printf(" WSUC - cosmic ray stand geometry \n");
165  }
166 
167  Int_t tracknumber = gAlice->GetMCApp()->GetCurrentTrackNumber();
168  Int_t parent=0;
169  TParticle* part=0;
170 
171  curVolName = fMC->CurrentVolName();
172  if(curVolName.Contains(vn) || curVolName.Contains("SCX"))
173  { // We are in a scintillator layer; SCX for 3X3
174 
175  if( ((depositedEnergy = fMC->Edep()) > 0.) && (fMC->TrackTime() < fTimeCut))
176  {// Track is inside a scintillator and deposits some energy
177  // Info("StepManager "," entry %i DE %f",++ientry, depositedEnergy); // for testing
178  if (fCurPrimary==-1)
179  fCurPrimary=gAlice->GetMCApp()->GetPrimary(tracknumber);
180 
181  if (fCurParent==-1 || tracknumber != fCurTrack)
182  {
183  // Check parentage
184  parent=tracknumber;
185 
186  if (fCurParent != -1)
187  {
188  while (parent != fCurParent && parent != -1)
189  {
190  //TParticle *part=gAlice->GetMCApp()->Particle(parent);
191  part=gAlice->GetMCApp()->Particle(parent);
192  parent=part->GetFirstMother();
193  }
194  }
195 
196  if (fCurParent==-1 || parent==-1)
197  {
198  //Int_t parent=tracknumber;
199  //TParticle *part=gAlice->GetMCApp()->Particle(parent);
200  parent=tracknumber;
201  part=gAlice->GetMCApp()->Particle(parent);
202  while (parent != -1 && fGeometry->IsInEMCALOrDCAL(part->Vx(),part->Vy(),part->Vz())) {
203  parent=part->GetFirstMother();
204  if (parent!=-1)
205  part=gAlice->GetMCApp()->Particle(parent);
206  }
207 
208  fCurParent=parent;
209  if (fCurParent==-1)
210  Error("StepManager","Cannot find parent");
211  else
212  {
213  //TParticle *part=gAlice->GetMCApp()->Particle(fCurParent);
215  ienergy = part->Energy();
216 
217  //Add reference to parent in TR tree.
219 
220  }
221  while (parent != -1)
222  {
223  part=gAlice->GetMCApp()->Particle(parent);
224  part->SetBit(kKeepBit);
225  parent=part->GetFirstMother();
226  }
227  }
228 
229  fCurTrack=tracknumber;
230  }
231 
232  fMC->TrackPosition(pos);
233  xyzte[0] = pos[0];
234  xyzte[1] = pos[1];
235  xyzte[2] = pos[2];
236  xyzte[3] = fMC->TrackTime() ;
237 
238  fMC->TrackMomentum(mom);
239  pmom[0] = mom[0];
240  pmom[1] = mom[1];
241  pmom[2] = mom[2];
242  pmom[3] = mom[3];
243 
244  // if(ientry%200 > 0) return; // testing
245  supModuleNumber = moduleNumber = yNumber = xNumber = absid = 0;
246  if(keyGeom >= 1)
247  { // TRD1 case now
248  fMC->CurrentVolOffID(4, supModuleNumber);
249  fMC->CurrentVolOffID(3, moduleNumber);
250  fMC->CurrentVolOffID(1, yNumber);
251  fMC->CurrentVolOffID(0, xNumber); // really x number now
252 
253  Int_t CurrentSMType = 0;
254  if (strcmp(fMC->CurrentVolOffName(4),"SMOD")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_Standard ;
255  else if(strcmp(fMC->CurrentVolOffName(4),"SM10")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_Half ;
256  else if(strcmp(fMC->CurrentVolOffName(4),"SM3rd")==0) CurrentSMType = AliEMCALGeometry::kEMCAL_3rd ;
257  else if(strcmp(fMC->CurrentVolOffName(4),"DCSM")==0) CurrentSMType = AliEMCALGeometry::kDCAL_Standard ;
258  else if(strcmp(fMC->CurrentVolOffName(4),"DCEXT")==0) CurrentSMType = AliEMCALGeometry::kDCAL_Ext ;
259  else AliError("Unkown SM Type!!");
260 
261  Int_t preSM = 0;
262  while( fGeometry->GetSMType(preSM) != CurrentSMType ) preSM++;
263  supModuleNumber += preSM;
264 
265  // Nov 10,2006
266  if(strcmp(fMC->CurrentVolOffName(0),vn) != 0) { // 3X3 case
267  if (strcmp(fMC->CurrentVolOffName(0),"SCX1")==0) xNumber=1;
268  else if(strcmp(fMC->CurrentVolOffName(0),"SCX2")==0) xNumber=2;
269  else if(strcmp(fMC->CurrentVolOffName(0),"SCX3")==0) xNumber=3;
270  else Fatal("StepManager()", "Wrong name of sensitive volume in 3X3 case : %s ", fMC->CurrentVolOffName(0));
271  }
272  }
273  else
274  {
275  fMC->CurrentVolOffID(5, supModuleNumber);
276  fMC->CurrentVolOffID(4, moduleNumber);
277  fMC->CurrentVolOffID(1, yNumber);
278  fMC->CurrentVolOffID(0, xNumber);
279 
280  if (strcmp(fMC->CurrentVolOffName(5),"SMOP")==0) supModuleNumber = 2*(supModuleNumber-1)+1;
281  else if(strcmp(fMC->CurrentVolOffName(5),"SMON")==0) supModuleNumber = 2*(supModuleNumber-1)+2;
282  else assert(0); // something wrong
283  }
284 
285  // Due to problem with index ordering conventions the calcultation of absid is no more like this:
286  // absid = fGeometry->GetAbsCellId(smNumber, moduleNumber-1, yNumber-1, xNumber-1);
287 
288  //
289  // Swap A side in Phi and C side in Eta due to wrong indexing.
290  //
291  Int_t iphi = -1;
292  Int_t ieta = -1;
293  Int_t smNumber = supModuleNumber-1;
294  Int_t smType = 1;
295  fGeometry->GetCellPhiEtaIndexInSModule(smNumber,moduleNumber-1,yNumber-1,xNumber-1, iphi, ieta);
296 
297  if (smNumber%2 == 0)
298  {
299  if(strcmp(fMC->CurrentVolOffName(4),"DCSM")==0) smType = 3; //DCal supermodule. previous design/idea
300  else smType = 2;
301  ieta = ((fGeometry->GetCentersOfCellsEtaDir()).GetSize()* 2/smType -1)-ieta;// 47/31-ieta, revert the ordering on A side in order to keep convention.
302  }
303  else
304  {
305  if(strcmp(fMC->CurrentVolOffName(4),"SM10")==0) smType = 2 ; //half supermodule. previous design/idea
306  if(strcmp(fMC->CurrentVolOffName(4),"SM3rd")==0) smType = 3 ; //one third (installed in 2012) supermodule
307  if(strcmp(fMC->CurrentVolOffName(4),"DCEXT")==0) smType = 3 ; //one third (installed in 2012) supermodule
308  iphi= ((fGeometry->GetCentersOfCellsPhiDir()).GetSize()/smType-1)-iphi;// 23/7-iphi, revert the ordering on C side in order to keep convention.
309  }
310 
311  // Once we know the indexes, calculate the absolute ID
312  absid = fGeometry->GetAbsCellIdFromCellIndexes(smNumber, iphi, ieta);
313 
314  if (absid < 0)
315  {
316  printf(" supModuleNumber %i : moduleNumber %i : yNumber %i : xNumber %i \n",
317  supModuleNumber, moduleNumber, yNumber, xNumber);
318  Fatal("StepManager()", "Wrong id : %i ", absid) ;
319  }
320 
321  Float_t lightYield = depositedEnergy ;
322  // Apply Birk's law (copied from G3BIRK)
323 
324  if (fMC->TrackCharge()!=0)
325  { // Check
326  Float_t birkC1Mod = 0;
327  if (fBirkC0==1)
328  { // Apply correction for higher charge states
329  if (TMath::Abs(fMC->TrackCharge())>=2) birkC1Mod = fBirkC1*7.2/12.6;
330  else birkC1Mod = fBirkC1;
331  }
332 
333  Float_t dedxcm=0.;
334  if (fMC->TrackStep()>0) dedxcm=1000.*fMC->Edep()/fMC->TrackStep();
335  else dedxcm=0;
336  lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm);
337  }
338 
339  // use sampling fraction to get original energy --HG
340  // xyzte[4] = lightYield * fGeometry->GetSampling();
341  xyzte[4] = lightYield; // 15-dec-04
342 
343  if (gDebug == -2)
344  printf("#sm %2i #m %3i #x %1i #z %1i -> absid %i : xyzte[4] = %f\n",
345  supModuleNumber,moduleNumber,yNumber,xNumber,absid, xyzte[4]);
346 
347  AddHit(fIshunt, fCurPrimary,tracknumber, fCurParent, ienergy, absid, xyzte, pmom);
348  } // there is deposited energy
349  }
350 }
TClonesArray * fHits
Counter for the hit iterator.
Definition: AliDetector.h:79
virtual void AddHit(Int_t shunt, Int_t primary, Int_t track, Int_t iparent, Float_t ienergy, Int_t id, Float_t *hits, Float_t *p)
Definition: AliEMCALv2.cxx:96
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
virtual void StepManager(void)
Accumulates hits as long as the track stays in a tower.
Definition: AliEMCALv2.cxx:138
TArrayD GetCentersOfCellsEtaDir() const
Int_t fCurTrack
Current track.
Definition: AliEMCALv1.h:63
virtual Int_t GetPrimary(Int_t track) const
Definition: AliMC.cxx:1888
Int_t GetAbsCellIdFromCellIndexes(Int_t nSupMod, Int_t iphi, Int_t ieta) const
EMCal simulation manager class v1.
Definition: AliEMCALv1.h:35
TVirtualMC * fMC
Definition: AliModule.h:155
Int_t fNhits
Definition: AliDetector.h:74
Float_t p[]
Definition: kNNTest.C:133
Int_t fCurParent
Current parent.
Definition: AliEMCALv1.h:62
void GetCellPhiEtaIndexInSModule(Int_t nSupMod, Int_t nModule, Int_t nIphi, Int_t nIeta, Int_t &iphi, Int_t &ieta) const
virtual ~AliEMCALv2(void)
Definition: AliEMCALv2.cxx:77
Int_t fCurPrimary
Current primary track.
Definition: AliEMCALv1.h:61
AliEMCALv2(void)
Default Constructor.
Definition: AliEMCALv2.cxx:44
EMCal hits object.
Definition: AliEMCALHit.h:24
Bool_t checkGeoAndRun
check or not the year to configure the detector
Definition: Config.C:97
#define AliInfo(message)
Definition: AliLog.h:484
Int_t GetSMType(Int_t nSupMod) const
AliRun * gAlice
Definition: AliRun.cxx:62
Int_t GetCurrentTrackNumber() const
Definition: AliMC.cxx:1846
Double_t fBirkC1
Constant 1 for Birk&#39;s Law implementation.
Definition: AliEMCAL.h:80
Int_t fBirkC0
Constant 0 for Birk&#39;s Law implementation.
Definition: AliEMCAL.h:79
Float_t fTimeCut
Cut to remove the background from the ALICE system.
Definition: AliEMCALv1.h:64
virtual AliTrackReference * AddTrackReference(Int_t label, Int_t id=-999)
Definition: AliModule.cxx:378
Int_t GetPrimary(void) const
Definition: AliEMCALHit.h:51
Int_t fIshunt
Definition: AliDetector.h:73
#define AliError(message)
Definition: AliLog.h:591
TArrayD GetCentersOfCellsPhiDir() const
Double_t fBirkC2
Constant 2 for Birk&#39;s Law implementation.
Definition: AliEMCAL.h:81
EMCal simulation manager class v2.
Definition: AliEMCALv2.h:30
AliMC * GetMCApp() const
Definition: AliRun.h:53
virtual void AddHitList(TCollection *hitList)
Definition: AliMC.h:83
Int_t IsInEMCALOrDCAL(Double_t x, Double_t y, Double_t z) const
TParticle * Particle(Int_t i) const
Definition: AliMC.cxx:1901
AliEMCALGeometry * fGeometry
! EMCal geometry access
Definition: AliEMCAL.h:85