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.108 2007/02/02 09:40:50 alibrary
22 * Includes required by ROOT head
24 * Revision 1.107 2007/02/01 10:34:47 hristov
25 * Removing warnings on Solaris x86
27 * Revision 1.106 2006/11/14 17:11:15 hristov
28 * 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
30 * Revision 1.105 2006/09/13 07:31:01 kharlov
31 * Effective C++ corrections (T.Pocheptsov)
33 * Revision 1.104 2005/05/28 14:19:05 schutz
34 * Compilation warnings fixed by T.P.
38 //_________________________________________________________________________
39 // Implementation version v1 of PHOS Manager class
42 // Layout EMC + CPV has name IHEP:
43 // Produces hits for CPV, cumulated hits
46 //*-- Author: Yves Schutz (SUBATECH)
49 // --- ROOT system ---
50 #include <TClonesArray.h>
51 #include <TParticle.h>
52 #include <TVirtualMC.h>
54 // --- Standard library ---
57 // --- AliRoot header files ---
58 #include "AliPHOSCPVDigit.h"
59 #include "AliPHOSGeometry.h"
60 #include "AliPHOSHit.h"
61 #include "AliPHOSv1.h"
67 //____________________________________________________________________________
68 AliPHOSv1::AliPHOSv1():
70 fIntrinsicPINEfficiency(0.),
71 fLightYieldAttenuation(0.),
72 fRecalibrationFactor(0.),
81 //____________________________________________________________________________
82 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
83 AliPHOSv0(name,title),
85 fIntrinsicPINEfficiency(0.),
86 fLightYieldAttenuation(0.),
87 fRecalibrationFactor(0.),
95 // - fHits (the "normal" one), which retains the hits associated with
96 // the current primary particle being tracked
97 // (this array is reset after each primary has been tracked).
102 // We do not want to save in TreeH the raw hits
103 // But save the cumulated hits instead (need to create the branch myself)
104 // It is put in the Digit Tree because the TreeH is filled after each primary
105 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
107 fHits= new TClonesArray("AliPHOSHit",1000) ;
108 gAlice->GetMCApp()->AddHitList(fHits) ;
112 fIshunt = 2 ; // All hits are associated with primary particles
114 //Photoelectron statistics:
115 // The light yield is a poissonian distribution of the number of
116 // photons created in the PbWo4 crystal, calculated using following formula
117 // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
118 // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
119 // LightYieldMean is parameter calculated to be over 47000 photons per GeV
120 // APDEfficiency is 0.02655
121 // k_0 is 0.0045 from Valery Antonenko
122 // The number of electrons created in the APD is
123 // NumberOfElectrons = APDGain * LightYield
124 // The APD Gain is 300
125 fLightYieldMean = 47000;
126 fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
127 fLightYieldAttenuation = 0.0045 ;
128 fRecalibrationFactor = 13.418/ fLightYieldMean ;
129 fElectronsPerGeV = 2.77e+8 ;
131 fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
132 fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
135 //____________________________________________________________________________
136 AliPHOSv1::~AliPHOSv1()
146 //____________________________________________________________________________
147 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
149 // Add a hit to the hit list.
150 // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
155 Bool_t deja = kFALSE ;
156 AliPHOSGeometry * geom = GetGeometry() ;
158 newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
160 for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
161 curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
162 if(curHit->GetPrimary() != primary) break ;
163 // We add hits with the same primary, while GEANT treats primaries succesively
164 if( *curHit == *newHit ) {
171 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
172 // get the block Id number
174 geom->AbsToRelNumbering(Id, relid) ;
182 //____________________________________________________________________________
183 void AliPHOSv1::FinishPrimary()
185 // called at the end of each track (primary) by AliRun
186 // hits are reset for each new track
187 // accumulate the total hit-multiplicity
191 //____________________________________________________________________________
192 void AliPHOSv1::FinishEvent()
194 // called at the end of each event by AliRun
195 // accumulate the hit-multiplicity and total energy per block
196 // if the values have been updated check it
198 AliDetector::FinishEvent();
200 //____________________________________________________________________________
201 void AliPHOSv1::StepManager(void)
203 // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
205 Int_t relid[4] ; // (box, layer, row, column) indices
206 Int_t absid ; // absolute cell ID number
207 Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
208 TLorentzVector pos ; // Lorentz vector of the track current position
211 TString name = GetGeometry()->GetName() ;
215 static Int_t idPCPQ = gMC->VolId("PCPQ");
216 if( gMC->CurrentVolID(copy) == idPCPQ &&
217 (gMC->IsTrackEntering() ) &&
218 gMC->TrackCharge() != 0) {
220 gMC -> TrackPosition(pos);
222 Float_t xyzm[3], xyzd[3] ;
224 for (i=0; i<3; i++) xyzm[i] = pos[i];
225 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
227 Float_t xyd[3]={0,0,0} ; //local position of the entering
232 // Current momentum of the hit's track in the local ref. system
233 TLorentzVector pmom ; //momentum of the particle initiated hit
234 gMC -> TrackMomentum(pmom);
235 Float_t pm[3], pd[3];
239 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
244 // Digitize the current CPV hit:
246 // 1. find pad response and
247 gMC->CurrentVolOffID(3,moduleNumber);
250 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
251 CPVDigitize(pmom,xyd,cpvDigits);
256 Int_t idigit,ndigits;
258 // 2. go through the current digit list and sum digits in pads
260 ndigits = cpvDigits->GetEntriesFast();
261 for (idigit=0; idigit<ndigits-1; idigit++) {
262 AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
263 Float_t x1 = cpvDigit1->GetXpad() ;
264 Float_t z1 = cpvDigit1->GetYpad() ;
265 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
266 AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
267 Float_t x2 = cpvDigit2->GetXpad() ;
268 Float_t z2 = cpvDigit2->GetYpad() ;
269 if (x1==x2 && z1==z2) {
270 Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
271 cpvDigit2->SetQpad(qsumpad) ;
272 cpvDigits->RemoveAt(idigit) ;
276 cpvDigits->Compress() ;
278 // 3. add digits to temporary hit list fTmpHits
280 ndigits = cpvDigits->GetEntriesFast();
281 for (idigit=0; idigit<ndigits; idigit++) {
282 AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
283 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
284 relid[1] =-1 ; // means CPV
285 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
286 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
288 // get the absolute Id number
289 GetGeometry()->RelToAbsNumbering(relid, absid) ;
291 // add current digit to the temporary hit list
293 xyzte[3] = gMC->TrackTime() ;
294 xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
296 Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
297 AddHit(fIshunt, primary, absid, xyzte);
299 if (cpvDigit->GetQpad() > 0.02) {
300 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
301 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
302 qsum += cpvDigit->GetQpad();
313 static Int_t idPXTL = gMC->VolId("PXTL");
314 if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
316 gMC->TrackPosition(pos) ;
321 Float_t global[3], local[3] ;
325 Float_t lostenergy = gMC->Edep();
327 //Put in the TreeK particle entering PHOS and all its parents
328 if ( gMC->IsTrackEntering() ){
330 gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
331 if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
332 Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
333 TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
334 Float_t vert[3],vertd[3] ;
338 gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
339 if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
340 //0.1 to get rid of numerical errors
341 part->SetBit(kKeepBit);
342 while ( parent != -1 ) {
343 part = gAlice->GetMCApp()->Particle(parent) ;
344 part->SetBit(kKeepBit);
345 parent = part->GetFirstMother() ;
350 if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
351 xyzte[3] = gMC->TrackTime() ;
353 gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
356 gMC->CurrentVolOffID(3, strip);
358 gMC->CurrentVolOffID(2, cell);
360 //Old formula for row is wrong. For example, I have strip 56 (28 for 2 x 8), row must be 1.
361 //But row == 1 + 56 - 56 % 56 == 57 (row == 1 + 28 - 28 % 28 == 29)
362 //Int_t row = 1 + GetGeometry()->GetEMCAGeometry()->GetNStripZ() - strip % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
363 Int_t row = GetGeometry()->GetEMCAGeometry()->GetNStripZ() - (strip - 1) % (GetGeometry()->GetEMCAGeometry()->GetNStripZ()) ;
364 Int_t col = (Int_t) TMath::Ceil((Double_t) strip/(GetGeometry()->GetEMCAGeometry()->GetNStripZ())) -1 ;
366 // Absid for 8x2-strips. Looks nice :)
367 absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
368 row * 2 + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsXInStrip() + (cell - 1) / 2)*GetGeometry()->GetNZ() - (cell & 1 ? 1 : 0);
370 gMC->Gmtod(global, local, 1) ;
372 //Calculates the light yield, the number of photons produced in the
374 Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
375 exp(-fLightYieldAttenuation *
376 (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
379 //Calculates de energy deposited in the crystal
380 xyzte[4] = fAPDFactor * lightYield ;
384 primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
385 TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
386 while ( !part->TestBit(kKeepBit) ) {
387 primary = part->GetFirstMother() ;
389 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
390 break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
391 //surface of the crystal. In this case it may have no primary at all.
392 //We can not easily separate this case from the case when this is part of the shower,
393 //developed in the neighboring crystal.
395 part = gAlice->GetMCApp()->Particle(primary) ;
399 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
403 // add current hit to the hit list
404 // Info("StepManager","%d %d", primary, tracknumber) ;
405 AddHit(fIshunt, primary, absid, xyzte);
407 } // there is deposited energy
408 } // we are inside a PHOS Xtal
412 //____________________________________________________________________________
413 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
415 // ------------------------------------------------------------------------
416 // Digitize one CPV hit:
417 // On input take exact 4-momentum p and position zxhit of the hit,
418 // find the pad response around this hit and
419 // put the amplitudes in the pads into array digits
421 // Author: Yuri Kharlov (after Serguei Sadovsky)
423 // ------------------------------------------------------------------------
425 const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
426 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
427 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
428 const Int_t kNgamz = 5; // Ionization size in Z
429 const Int_t kNgamx = 9; // Ionization size in Phi
430 const Float_t kNoise = 0.03; // charge noise in one pad
434 // Just a reminder on axes notation in the CPV module:
435 // axis Z goes along the beam
436 // axis X goes across the beam in the module plane
437 // axis Y is a normal to the module plane showing from the IP
439 Float_t hitX = zxhit[0];
440 Float_t hitZ =-zxhit[1];
443 Float_t pNorm = p.Py();
444 Float_t eloss = kdEdx;
446 //Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
448 Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
449 Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
450 gRandom->Rannor(rnor1,rnor2);
451 eloss *= (1 + kDetR*rnor1) *
452 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
453 Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
454 Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
455 Float_t zhit2 = zhit1 + dZY;
456 Float_t xhit2 = xhit1 + dXY;
458 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
459 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
463 if (iwht1==iwht2) { // incline 1-wire hit
465 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
466 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
468 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
469 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
472 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
474 Int_t iwht3 = (iwht1 + iwht2) / 2;
475 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
476 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
477 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
478 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
479 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
480 Float_t dxw1 = xhit1 - xwr13;
481 Float_t dxw2 = xhit2 - xwr23;
482 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
483 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
484 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
485 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
487 zxe[2][0] = eloss * egm1;
488 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
490 zxe[2][1] = eloss * egm2;
491 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
493 zxe[2][2] = eloss * egm3;
495 else { // incline 2-wire hit
497 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
498 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
499 Float_t xwr12 = (xwht1 + xwht2) / 2;
500 Float_t dxw1 = xhit1 - xwr12;
501 Float_t dxw2 = xhit2 - xwr12;
502 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
503 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
504 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
506 zxe[2][0] = eloss * egm1;
507 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
509 zxe[2][1] = eloss * egm2;
512 // Finite size of ionization region
514 Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
515 Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
516 Int_t nz3 = (kNgamz+1)/2;
517 Int_t nx3 = (kNgamx+1)/2;
518 cpvDigits->Expand(nIter*kNgamx*kNgamz);
519 TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
521 for (Int_t iter=0; iter<nIter; iter++) {
523 Float_t zhit = zxe[0][iter];
524 Float_t xhit = zxe[1][iter];
525 Float_t qhit = zxe[2][iter];
526 Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
527 Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
528 if ( zcell<=0 || xcell<=0 ||
529 zcell>=nCellZ || xcell>=nCellX) return;
530 Int_t izcell = (Int_t) zcell;
531 Int_t ixcell = (Int_t) xcell;
532 Float_t zc = zcell - izcell - 0.5;
533 Float_t xc = xcell - ixcell - 0.5;
534 for (Int_t iz=1; iz<=kNgamz; iz++) {
535 Int_t kzg = izcell + iz - nz3;
536 if (kzg<=0 || kzg>nCellZ) continue;
537 Float_t zg = (Float_t)(iz-nz3) - zc;
538 for (Int_t ix=1; ix<=kNgamx; ix++) {
539 Int_t kxg = ixcell + ix - nx3;
540 if (kxg<=0 || kxg>nCellX) continue;
541 Float_t xg = (Float_t)(ix-nx3) - xc;
543 // Now calculate pad response
544 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
545 qpad += kNoise*rnor2;
546 if (qpad<0) continue;
548 // Fill the array with pad response ID and amplitude
549 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
555 //____________________________________________________________________________
556 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
557 // ------------------------------------------------------------------------
558 // Calculate the amplitude in one CPV pad using the
559 // cumulative pad response function
560 // Author: Yuri Kharlov (after Serguei Sadovski)
562 // ------------------------------------------------------------------------
564 Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
565 Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
566 Double_t z = zhit * GetGeometry()->GetPadSizeZ();
567 Double_t x = xhit * GetGeometry()->GetPadSizePhi();
568 Double_t amplitude = qhit *
569 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
570 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
571 return (Float_t)amplitude;
574 //____________________________________________________________________________
575 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
576 // ------------------------------------------------------------------------
577 // Cumulative pad response function
578 // It includes several terms from the CF decomposition in electrostatics
579 // Note: this cumulative function is wrong since omits some terms
580 // but the cell amplitude obtained with it is correct because
581 // these omitting terms cancel
582 // Author: Yuri Kharlov (after Serguei Sadovski)
584 // ------------------------------------------------------------------------
586 const Double_t kA=1.0;
587 const Double_t kB=0.7;
589 Double_t r2 = x*x + y*y;
591 Double_t cumulPRF = 0;
592 for (Int_t i=0; i<=4; i++) {
593 Double_t b1 = (2*i + 1) * kB;
594 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
596 cumulPRF *= kA/(2*TMath::Pi());