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.107 2007/02/01 10:34:47 hristov
22 * Removing warnings on Solaris x86
24 * Revision 1.106 2006/11/14 17:11:15 hristov
25 * 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
27 * Revision 1.105 2006/09/13 07:31:01 kharlov
28 * Effective C++ corrections (T.Pocheptsov)
30 * Revision 1.104 2005/05/28 14:19:05 schutz
31 * Compilation warnings fixed by T.P.
35 //_________________________________________________________________________
36 // Implementation version v1 of PHOS Manager class
39 // Layout EMC + CPV has name IHEP:
40 // Produces hits for CPV, cumulated hits
43 //*-- Author: Yves Schutz (SUBATECH)
46 // --- ROOT system ---
47 #include <TClonesArray.h>
48 #include <TParticle.h>
49 #include <TVirtualMC.h>
51 // --- Standard library ---
54 // --- AliRoot header files ---
55 #include "AliPHOSCPVDigit.h"
56 #include "AliPHOSGeometry.h"
57 #include "AliPHOSHit.h"
58 #include "AliPHOSv1.h"
64 //____________________________________________________________________________
65 AliPHOSv1::AliPHOSv1():
67 fIntrinsicPINEfficiency(0.),
68 fLightYieldAttenuation(0.),
69 fRecalibrationFactor(0.),
78 //____________________________________________________________________________
79 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
80 AliPHOSv0(name,title),
82 fIntrinsicPINEfficiency(0.),
83 fLightYieldAttenuation(0.),
84 fRecalibrationFactor(0.),
92 // - fHits (the "normal" one), which retains the hits associated with
93 // the current primary particle being tracked
94 // (this array is reset after each primary has been tracked).
99 // We do not want to save in TreeH the raw hits
100 // But save the cumulated hits instead (need to create the branch myself)
101 // It is put in the Digit Tree because the TreeH is filled after each primary
102 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
104 fHits= new TClonesArray("AliPHOSHit",1000) ;
105 gAlice->GetMCApp()->AddHitList(fHits) ;
109 fIshunt = 2 ; // All hits are associated with primary particles
111 //Photoelectron statistics:
112 // The light yield is a poissonian distribution of the number of
113 // photons created in the PbWo4 crystal, calculated using following formula
114 // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
115 // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
116 // LightYieldMean is parameter calculated to be over 47000 photons per GeV
117 // APDEfficiency is 0.02655
118 // k_0 is 0.0045 from Valery Antonenko
119 // The number of electrons created in the APD is
120 // NumberOfElectrons = APDGain * LightYield
121 // The APD Gain is 300
122 fLightYieldMean = 47000;
123 fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
124 fLightYieldAttenuation = 0.0045 ;
125 fRecalibrationFactor = 13.418/ fLightYieldMean ;
126 fElectronsPerGeV = 2.77e+8 ;
128 fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
129 fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
132 //____________________________________________________________________________
133 AliPHOSv1::~AliPHOSv1()
143 //____________________________________________________________________________
144 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
146 // Add a hit to the hit list.
147 // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
152 Bool_t deja = kFALSE ;
153 AliPHOSGeometry * geom = GetGeometry() ;
155 newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
157 for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
158 curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
159 if(curHit->GetPrimary() != primary) break ;
160 // We add hits with the same primary, while GEANT treats primaries succesively
161 if( *curHit == *newHit ) {
168 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
169 // get the block Id number
171 geom->AbsToRelNumbering(Id, relid) ;
179 //____________________________________________________________________________
180 void AliPHOSv1::FinishPrimary()
182 // called at the end of each track (primary) by AliRun
183 // hits are reset for each new track
184 // accumulate the total hit-multiplicity
188 //____________________________________________________________________________
189 void AliPHOSv1::FinishEvent()
191 // called at the end of each event by AliRun
192 // accumulate the hit-multiplicity and total energy per block
193 // if the values have been updated check it
195 AliDetector::FinishEvent();
197 //____________________________________________________________________________
198 void AliPHOSv1::StepManager(void)
200 // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
202 Int_t relid[4] ; // (box, layer, row, column) indices
203 Int_t absid ; // absolute cell ID number
204 Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
205 TLorentzVector pos ; // Lorentz vector of the track current position
208 TString name = GetGeometry()->GetName() ;
212 static Int_t idPCPQ = gMC->VolId("PCPQ");
213 if( gMC->CurrentVolID(copy) == idPCPQ &&
214 (gMC->IsTrackEntering() ) &&
215 gMC->TrackCharge() != 0) {
217 gMC -> TrackPosition(pos);
219 Float_t xyzm[3], xyzd[3] ;
221 for (i=0; i<3; i++) xyzm[i] = pos[i];
222 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
224 Float_t xyd[3]={0,0,0} ; //local position of the entering
229 // Current momentum of the hit's track in the local ref. system
230 TLorentzVector pmom ; //momentum of the particle initiated hit
231 gMC -> TrackMomentum(pmom);
232 Float_t pm[3], pd[3];
236 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
241 // Digitize the current CPV hit:
243 // 1. find pad response and
244 gMC->CurrentVolOffID(3,moduleNumber);
247 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
248 CPVDigitize(pmom,xyd,cpvDigits);
253 Int_t idigit,ndigits;
255 // 2. go through the current digit list and sum digits in pads
257 ndigits = cpvDigits->GetEntriesFast();
258 for (idigit=0; idigit<ndigits-1; idigit++) {
259 AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
260 Float_t x1 = cpvDigit1->GetXpad() ;
261 Float_t z1 = cpvDigit1->GetYpad() ;
262 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
263 AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
264 Float_t x2 = cpvDigit2->GetXpad() ;
265 Float_t z2 = cpvDigit2->GetYpad() ;
266 if (x1==x2 && z1==z2) {
267 Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
268 cpvDigit2->SetQpad(qsumpad) ;
269 cpvDigits->RemoveAt(idigit) ;
273 cpvDigits->Compress() ;
275 // 3. add digits to temporary hit list fTmpHits
277 ndigits = cpvDigits->GetEntriesFast();
278 for (idigit=0; idigit<ndigits; idigit++) {
279 AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
280 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
281 relid[1] =-1 ; // means CPV
282 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
283 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
285 // get the absolute Id number
286 GetGeometry()->RelToAbsNumbering(relid, absid) ;
288 // add current digit to the temporary hit list
290 xyzte[3] = gMC->TrackTime() ;
291 xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
293 Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
294 AddHit(fIshunt, primary, absid, xyzte);
296 if (cpvDigit->GetQpad() > 0.02) {
297 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
298 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
299 qsum += cpvDigit->GetQpad();
310 static Int_t idPXTL = gMC->VolId("PXTL");
311 if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
313 gMC->TrackPosition(pos) ;
318 Float_t global[3], local[3] ;
322 Float_t lostenergy = gMC->Edep();
324 //Put in the TreeK particle entering PHOS and all its parents
325 if ( gMC->IsTrackEntering() ){
327 gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
328 if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
329 Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
330 TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
331 Float_t vert[3],vertd[3] ;
335 gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
336 if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
337 //0.1 to get rid of numerical errors
338 part->SetBit(kKeepBit);
339 while ( parent != -1 ) {
340 part = gAlice->GetMCApp()->Particle(parent) ;
341 part->SetBit(kKeepBit);
342 parent = part->GetFirstMother() ;
347 if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
348 xyzte[3] = gMC->TrackTime() ;
350 gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
353 gMC->CurrentVolOffID(3, strip);
355 gMC->CurrentVolOffID(2, cell);
357 Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
358 Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
360 absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
361 row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
363 gMC->Gmtod(global, local, 1) ;
365 //Calculates the light yield, the number of photons produced in the
367 Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
368 exp(-fLightYieldAttenuation *
369 (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
372 //Calculates de energy deposited in the crystal
373 xyzte[4] = fAPDFactor * lightYield ;
377 primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
378 TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
379 while ( !part->TestBit(kKeepBit) ) {
380 primary = part->GetFirstMother() ;
382 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
383 break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
384 //surface of the crystal. In this case it may have no primary at all.
385 //We can not easily separate this case from the case when this is part of the shower,
386 //developed in the neighboring crystal.
388 part = gAlice->GetMCApp()->Particle(primary) ;
392 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
396 // add current hit to the hit list
397 // Info("StepManager","%d %d", primary, tracknumber) ;
398 AddHit(fIshunt, primary, absid, xyzte);
400 } // there is deposited energy
401 } // we are inside a PHOS Xtal
405 //____________________________________________________________________________
406 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
408 // ------------------------------------------------------------------------
409 // Digitize one CPV hit:
410 // On input take exact 4-momentum p and position zxhit of the hit,
411 // find the pad response around this hit and
412 // put the amplitudes in the pads into array digits
414 // Author: Yuri Kharlov (after Serguei Sadovsky)
416 // ------------------------------------------------------------------------
418 const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
419 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
420 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
421 const Int_t kNgamz = 5; // Ionization size in Z
422 const Int_t kNgamx = 9; // Ionization size in Phi
423 const Float_t kNoise = 0.03; // charge noise in one pad
427 // Just a reminder on axes notation in the CPV module:
428 // axis Z goes along the beam
429 // axis X goes across the beam in the module plane
430 // axis Y is a normal to the module plane showing from the IP
432 Float_t hitX = zxhit[0];
433 Float_t hitZ =-zxhit[1];
436 Float_t pNorm = p.Py();
437 Float_t eloss = kdEdx;
439 //Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
441 Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
442 Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
443 gRandom->Rannor(rnor1,rnor2);
444 eloss *= (1 + kDetR*rnor1) *
445 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
446 Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
447 Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
448 Float_t zhit2 = zhit1 + dZY;
449 Float_t xhit2 = xhit1 + dXY;
451 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
452 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
456 if (iwht1==iwht2) { // incline 1-wire hit
458 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
459 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
461 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
462 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
465 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
467 Int_t iwht3 = (iwht1 + iwht2) / 2;
468 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
469 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
470 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
471 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
472 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
473 Float_t dxw1 = xhit1 - xwr13;
474 Float_t dxw2 = xhit2 - xwr23;
475 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
476 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
477 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
478 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
480 zxe[2][0] = eloss * egm1;
481 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
483 zxe[2][1] = eloss * egm2;
484 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
486 zxe[2][2] = eloss * egm3;
488 else { // incline 2-wire hit
490 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
491 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
492 Float_t xwr12 = (xwht1 + xwht2) / 2;
493 Float_t dxw1 = xhit1 - xwr12;
494 Float_t dxw2 = xhit2 - xwr12;
495 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
496 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
497 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
499 zxe[2][0] = eloss * egm1;
500 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
502 zxe[2][1] = eloss * egm2;
505 // Finite size of ionization region
507 Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
508 Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
509 Int_t nz3 = (kNgamz+1)/2;
510 Int_t nx3 = (kNgamx+1)/2;
511 cpvDigits->Expand(nIter*kNgamx*kNgamz);
512 TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
514 for (Int_t iter=0; iter<nIter; iter++) {
516 Float_t zhit = zxe[0][iter];
517 Float_t xhit = zxe[1][iter];
518 Float_t qhit = zxe[2][iter];
519 Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
520 Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
521 if ( zcell<=0 || xcell<=0 ||
522 zcell>=nCellZ || xcell>=nCellX) return;
523 Int_t izcell = (Int_t) zcell;
524 Int_t ixcell = (Int_t) xcell;
525 Float_t zc = zcell - izcell - 0.5;
526 Float_t xc = xcell - ixcell - 0.5;
527 for (Int_t iz=1; iz<=kNgamz; iz++) {
528 Int_t kzg = izcell + iz - nz3;
529 if (kzg<=0 || kzg>nCellZ) continue;
530 Float_t zg = (Float_t)(iz-nz3) - zc;
531 for (Int_t ix=1; ix<=kNgamx; ix++) {
532 Int_t kxg = ixcell + ix - nx3;
533 if (kxg<=0 || kxg>nCellX) continue;
534 Float_t xg = (Float_t)(ix-nx3) - xc;
536 // Now calculate pad response
537 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
538 qpad += kNoise*rnor2;
539 if (qpad<0) continue;
541 // Fill the array with pad response ID and amplitude
542 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
548 //____________________________________________________________________________
549 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
550 // ------------------------------------------------------------------------
551 // Calculate the amplitude in one CPV pad using the
552 // cumulative pad response function
553 // Author: Yuri Kharlov (after Serguei Sadovski)
555 // ------------------------------------------------------------------------
557 Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
558 Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
559 Double_t z = zhit * GetGeometry()->GetPadSizeZ();
560 Double_t x = xhit * GetGeometry()->GetPadSizePhi();
561 Double_t amplitude = qhit *
562 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
563 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
564 return (Float_t)amplitude;
567 //____________________________________________________________________________
568 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
569 // ------------------------------------------------------------------------
570 // Cumulative pad response function
571 // It includes several terms from the CF decomposition in electrostatics
572 // Note: this cumulative function is wrong since omits some terms
573 // but the cell amplitude obtained with it is correct because
574 // these omitting terms cancel
575 // Author: Yuri Kharlov (after Serguei Sadovski)
577 // ------------------------------------------------------------------------
579 const Double_t kA=1.0;
580 const Double_t kB=0.7;
582 Double_t r2 = x*x + y*y;
584 Double_t cumulPRF = 0;
585 for (Int_t i=0; i<=4; i++) {
586 Double_t b1 = (2*i + 1) * kB;
587 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
589 cumulPRF *= kA/(2*TMath::Pi());