1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.111 2007/07/24 09:41:19 morsch
22 * AliStack included for kKeepBit.
24 * Revision 1.110 2007/03/10 08:58:52 kharlov
25 * Protection for noCPV geometry
27 * Revision 1.109 2007/03/01 11:37:37 kharlov
28 * Strip units changed from 8x1 to 8x2 (T.Pocheptsov)
30 * Revision 1.108 2007/02/02 09:40:50 alibrary
31 * Includes required by ROOT head
33 * Revision 1.107 2007/02/01 10:34:47 hristov
34 * Removing warnings on Solaris x86
36 * Revision 1.106 2006/11/14 17:11:15 hristov
37 * Removing inheritances from TAttLine, TAttMarker and AliRndm in AliModule. The copy constructor and assignment operators are moved to the private part of the class and not implemented. The corresponding changes are propagated to the detectors
39 * Revision 1.105 2006/09/13 07:31:01 kharlov
40 * Effective C++ corrections (T.Pocheptsov)
42 * Revision 1.104 2005/05/28 14:19:05 schutz
43 * Compilation warnings fixed by T.P.
47 //_________________________________________________________________________
48 // Implementation version v1 of PHOS Manager class
51 // Layout EMC + CPV has name IHEP:
52 // Produces hits for CPV, cumulated hits
55 //*-- Author: Yves Schutz (SUBATECH)
58 // --- ROOT system ---
59 #include <TClonesArray.h>
60 #include <TParticle.h>
61 #include <TVirtualMC.h>
63 // --- Standard library ---
66 // --- AliRoot header files ---
67 #include "AliPHOSCPVDigit.h"
68 #include "AliPHOSGeometry.h"
69 #include "AliPHOSHit.h"
70 #include "AliPHOSv1.h"
77 //____________________________________________________________________________
78 AliPHOSv1::AliPHOSv1():
80 fIntrinsicPINEfficiency(0.),
81 fLightYieldAttenuation(0.),
82 fRecalibrationFactor(0.),
91 //____________________________________________________________________________
92 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
93 AliPHOSv0(name,title),
95 fIntrinsicPINEfficiency(0.),
96 fLightYieldAttenuation(0.),
97 fRecalibrationFactor(0.),
105 // - fHits (the "normal" one), which retains the hits associated with
106 // the current primary particle being tracked
107 // (this array is reset after each primary has been tracked).
112 // We do not want to save in TreeH the raw hits
113 // But save the cumulated hits instead (need to create the branch myself)
114 // It is put in the Digit Tree because the TreeH is filled after each primary
115 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
117 fHits= new TClonesArray("AliPHOSHit",1000) ;
118 gAlice->GetMCApp()->AddHitList(fHits) ;
122 fIshunt = 2 ; // All hits are associated with primary particles
124 //Photoelectron statistics:
125 // The light yield is a poissonian distribution of the number of
126 // photons created in the PbWo4 crystal, calculated using following formula
127 // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
128 // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
129 // LightYieldMean is parameter calculated to be over 47000 photons per GeV
130 // APDEfficiency is 0.02655
131 // k_0 is 0.0045 from Valery Antonenko
132 // The number of electrons created in the APD is
133 // NumberOfElectrons = APDGain * LightYield
134 // The APD Gain is 300
135 fLightYieldMean = 47000;
136 fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
137 fLightYieldAttenuation = 0.0045 ;
138 fRecalibrationFactor = 13.418/ fLightYieldMean ;
139 fElectronsPerGeV = 2.77e+8 ;
141 fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
142 fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
145 //____________________________________________________________________________
146 AliPHOSv1::~AliPHOSv1()
156 //____________________________________________________________________________
157 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
159 // Add a hit to the hit list.
160 // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
165 Bool_t deja = kFALSE ;
166 AliPHOSGeometry * geom = GetGeometry() ;
168 newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
170 for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
171 curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
172 if(curHit->GetPrimary() != primary) break ;
173 // We add hits with the same primary, while GEANT treats primaries succesively
174 if( *curHit == *newHit ) {
181 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
182 // get the block Id number
184 geom->AbsToRelNumbering(Id, relid) ;
192 //____________________________________________________________________________
193 void AliPHOSv1::FinishPrimary()
195 // called at the end of each track (primary) by AliRun
196 // hits are reset for each new track
197 // accumulate the total hit-multiplicity
201 //____________________________________________________________________________
202 void AliPHOSv1::FinishEvent()
204 // called at the end of each event by AliRun
205 // accumulate the hit-multiplicity and total energy per block
206 // if the values have been updated check it
208 AliDetector::FinishEvent();
210 //____________________________________________________________________________
211 void AliPHOSv1::StepManager(void)
213 // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
215 Int_t relid[4] ; // (box, layer, row, column) indices
216 Int_t absid ; // absolute cell ID number
217 Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
218 TLorentzVector pos ; // Lorentz vector of the track current position
223 static Int_t idPCPQ = -1;
224 if (strstr(fTitle.Data(),"noCPV") == 0)
225 idPCPQ = gMC->VolId("PCPQ");
227 if( gMC->CurrentVolID(copy) == idPCPQ &&
228 (gMC->IsTrackEntering() ) &&
229 gMC->TrackCharge() != 0) {
231 gMC -> TrackPosition(pos);
233 Float_t xyzm[3], xyzd[3] ;
235 for (i=0; i<3; i++) xyzm[i] = pos[i];
236 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
239 Float_t xyd[3]={0,0,0} ; //local position of the entering
244 // Current momentum of the hit's track in the local ref. system
245 TLorentzVector pmom ; //momentum of the particle initiated hit
246 gMC -> TrackMomentum(pmom);
247 Float_t pm[3], pd[3];
251 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
256 // Digitize the current CPV hit:
258 // 1. find pad response and
259 gMC->CurrentVolOffID(3,moduleNumber);
262 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
263 CPVDigitize(pmom,xyd,cpvDigits);
268 Int_t idigit,ndigits;
270 // 2. go through the current digit list and sum digits in pads
272 ndigits = cpvDigits->GetEntriesFast();
273 for (idigit=0; idigit<ndigits-1; idigit++) {
274 AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
275 Float_t x1 = cpvDigit1->GetXpad() ;
276 Float_t z1 = cpvDigit1->GetYpad() ;
277 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
278 AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
279 Float_t x2 = cpvDigit2->GetXpad() ;
280 Float_t z2 = cpvDigit2->GetYpad() ;
281 if (x1==x2 && z1==z2) {
282 Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
283 cpvDigit2->SetQpad(qsumpad) ;
284 cpvDigits->RemoveAt(idigit) ;
288 cpvDigits->Compress() ;
290 // 3. add digits to temporary hit list fTmpHits
292 ndigits = cpvDigits->GetEntriesFast();
293 for (idigit=0; idigit<ndigits; idigit++) {
294 AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
295 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
296 relid[1] =-1 ; // means CPV
297 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
298 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
300 // get the absolute Id number
301 GetGeometry()->RelToAbsNumbering(relid, absid) ;
303 // add current digit to the temporary hit list
305 xyzte[3] = gMC->TrackTime() ;
306 xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
308 Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
309 AddHit(fIshunt, primary, absid, xyzte);
311 if (cpvDigit->GetQpad() > 0.02) {
312 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
313 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
314 qsum += cpvDigit->GetQpad();
325 static Int_t idPXTL = gMC->VolId("PXTL");
326 if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
328 gMC->TrackPosition(pos) ;
333 Float_t lostenergy = gMC->Edep();
335 //Put in the TreeK particle entering PHOS and all its parents
336 if ( gMC->IsTrackEntering() ){
338 gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
339 if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
340 Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
341 TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
342 Float_t vert[3],vertd[3] ;
346 gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
347 if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
348 //0.1 to get rid of numerical errors
349 part->SetBit(kKeepBit);
350 while ( parent != -1 ) {
351 part = gAlice->GetMCApp()->Particle(parent) ;
352 part->SetBit(kKeepBit);
353 parent = part->GetFirstMother() ;
358 if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
359 xyzte[3] = gMC->TrackTime() ;
361 gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
364 gMC->CurrentVolOffID(3, strip);
366 gMC->CurrentVolOffID(2, cell);
368 //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
369 //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
370 //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
371 Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
372 Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
374 // Absid for 8x2-strips. Looks nice :)
375 absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
376 row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
379 //Calculates the light yield, the number of photons produced in the
381 //There is no dependence of reponce on distance from energy deposition to APD
382 Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy) ;
384 //Calculates de energy deposited in the crystal
385 xyzte[4] = fAPDFactor * lightYield ;
389 primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
390 TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
391 while ( !part->TestBit(kKeepBit) ) {
392 primary = part->GetFirstMother() ;
394 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
395 break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
396 //surface of the crystal. In this case it may have no primary at all.
397 //We can not easily separate this case from the case when this is part of the shower,
398 //developed in the neighboring crystal.
400 part = gAlice->GetMCApp()->Particle(primary) ;
404 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
408 // add current hit to the hit list
409 // Info("StepManager","%d %d", primary, tracknumber) ;
410 AddHit(fIshunt, primary, absid, xyzte);
412 } // there is deposited energy
413 } // we are inside a PHOS Xtal
417 //____________________________________________________________________________
418 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
420 // ------------------------------------------------------------------------
421 // Digitize one CPV hit:
422 // On input take exact 4-momentum p and position zxhit of the hit,
423 // find the pad response around this hit and
424 // put the amplitudes in the pads into array digits
426 // Author: Yuri Kharlov (after Serguei Sadovsky)
428 // ------------------------------------------------------------------------
430 const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
431 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
432 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
433 const Int_t kNgamz = 5; // Ionization size in Z
434 const Int_t kNgamx = 9; // Ionization size in Phi
435 const Float_t kNoise = 0.03; // charge noise in one pad
439 // Just a reminder on axes notation in the CPV module:
440 // axis Z goes along the beam
441 // axis X goes across the beam in the module plane
442 // axis Y is a normal to the module plane showing from the IP
444 Float_t hitX = zxhit[0];
445 Float_t hitZ =-zxhit[1];
448 Float_t pNorm = p.Py();
449 Float_t eloss = kdEdx;
451 //Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
453 Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
454 Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
455 gRandom->Rannor(rnor1,rnor2);
456 eloss *= (1 + kDetR*rnor1) *
457 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
458 Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
459 Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
460 Float_t zhit2 = zhit1 + dZY;
461 Float_t xhit2 = xhit1 + dXY;
463 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
464 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
468 if (iwht1==iwht2) { // incline 1-wire hit
470 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
471 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
473 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
474 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
477 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
479 Int_t iwht3 = (iwht1 + iwht2) / 2;
480 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
481 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
482 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
483 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
484 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
485 Float_t dxw1 = xhit1 - xwr13;
486 Float_t dxw2 = xhit2 - xwr23;
487 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
488 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
489 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
490 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
492 zxe[2][0] = eloss * egm1;
493 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
495 zxe[2][1] = eloss * egm2;
496 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
498 zxe[2][2] = eloss * egm3;
500 else { // incline 2-wire hit
502 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
503 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
504 Float_t xwr12 = (xwht1 + xwht2) / 2;
505 Float_t dxw1 = xhit1 - xwr12;
506 Float_t dxw2 = xhit2 - xwr12;
507 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
508 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
509 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
511 zxe[2][0] = eloss * egm1;
512 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
514 zxe[2][1] = eloss * egm2;
517 // Finite size of ionization region
519 Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
520 Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
521 Int_t nz3 = (kNgamz+1)/2;
522 Int_t nx3 = (kNgamx+1)/2;
523 cpvDigits->Expand(nIter*kNgamx*kNgamz);
524 TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
526 for (Int_t iter=0; iter<nIter; iter++) {
528 Float_t zhit = zxe[0][iter];
529 Float_t xhit = zxe[1][iter];
530 Float_t qhit = zxe[2][iter];
531 Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
532 Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
533 if ( zcell<=0 || xcell<=0 ||
534 zcell>=nCellZ || xcell>=nCellX) return;
535 Int_t izcell = (Int_t) zcell;
536 Int_t ixcell = (Int_t) xcell;
537 Float_t zc = zcell - izcell - 0.5;
538 Float_t xc = xcell - ixcell - 0.5;
539 for (Int_t iz=1; iz<=kNgamz; iz++) {
540 Int_t kzg = izcell + iz - nz3;
541 if (kzg<=0 || kzg>nCellZ) continue;
542 Float_t zg = (Float_t)(iz-nz3) - zc;
543 for (Int_t ix=1; ix<=kNgamx; ix++) {
544 Int_t kxg = ixcell + ix - nx3;
545 if (kxg<=0 || kxg>nCellX) continue;
546 Float_t xg = (Float_t)(ix-nx3) - xc;
548 // Now calculate pad response
549 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
550 qpad += kNoise*rnor2;
551 if (qpad<0) continue;
553 // Fill the array with pad response ID and amplitude
554 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
560 //____________________________________________________________________________
561 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
562 // ------------------------------------------------------------------------
563 // Calculate the amplitude in one CPV pad using the
564 // cumulative pad response function
565 // Author: Yuri Kharlov (after Serguei Sadovski)
567 // ------------------------------------------------------------------------
569 Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
570 Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
571 Double_t z = zhit * GetGeometry()->GetPadSizeZ();
572 Double_t x = xhit * GetGeometry()->GetPadSizePhi();
573 Double_t amplitude = qhit *
574 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
575 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
576 return (Float_t)amplitude;
579 //____________________________________________________________________________
580 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
581 // ------------------------------------------------------------------------
582 // Cumulative pad response function
583 // It includes several terms from the CF decomposition in electrostatics
584 // Note: this cumulative function is wrong since omits some terms
585 // but the cell amplitude obtained with it is correct because
586 // these omitting terms cancel
587 // Author: Yuri Kharlov (after Serguei Sadovski)
589 // ------------------------------------------------------------------------
591 const Double_t kA=1.0;
592 const Double_t kB=0.7;
594 Double_t r2 = x*x + y*y;
596 Double_t cumulPRF = 0;
597 for (Int_t i=0; i<=4; i++) {
598 Double_t b1 = (2*i + 1) * kB;
599 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
601 cumulPRF *= kA/(2*TMath::Pi());