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.106 2006/11/14 17:11:15 hristov
22 * 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
24 * Revision 1.105 2006/09/13 07:31:01 kharlov
25 * Effective C++ corrections (T.Pocheptsov)
27 * Revision 1.104 2005/05/28 14:19:05 schutz
28 * Compilation warnings fixed by T.P.
32 //_________________________________________________________________________
33 // Implementation version v1 of PHOS Manager class
36 // Layout EMC + CPV has name IHEP:
37 // Produces hits for CPV, cumulated hits
40 //*-- Author: Yves Schutz (SUBATECH)
43 // --- ROOT system ---
44 #include <TParticle.h>
45 #include <TVirtualMC.h>
47 // --- Standard library ---
50 // --- AliRoot header files ---
51 #include "AliPHOSCPVDigit.h"
52 #include "AliPHOSGeometry.h"
53 #include "AliPHOSHit.h"
54 #include "AliPHOSv1.h"
60 //____________________________________________________________________________
61 AliPHOSv1::AliPHOSv1():
63 fIntrinsicPINEfficiency(0.),
64 fLightYieldAttenuation(0.),
65 fRecalibrationFactor(0.),
74 //____________________________________________________________________________
75 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
76 AliPHOSv0(name,title),
78 fIntrinsicPINEfficiency(0.),
79 fLightYieldAttenuation(0.),
80 fRecalibrationFactor(0.),
88 // - fHits (the "normal" one), which retains the hits associated with
89 // the current primary particle being tracked
90 // (this array is reset after each primary has been tracked).
95 // We do not want to save in TreeH the raw hits
96 // But save the cumulated hits instead (need to create the branch myself)
97 // It is put in the Digit Tree because the TreeH is filled after each primary
98 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
100 fHits= new TClonesArray("AliPHOSHit",1000) ;
101 gAlice->GetMCApp()->AddHitList(fHits) ;
105 fIshunt = 2 ; // All hits are associated with primary particles
107 //Photoelectron statistics:
108 // The light yield is a poissonian distribution of the number of
109 // photons created in the PbWo4 crystal, calculated using following formula
110 // NumberOfPhotons = EnergyLost * LightYieldMean* APDEfficiency *
111 // exp (-LightYieldAttenuation * DistanceToPINdiodeFromTheHit);
112 // LightYieldMean is parameter calculated to be over 47000 photons per GeV
113 // APDEfficiency is 0.02655
114 // k_0 is 0.0045 from Valery Antonenko
115 // The number of electrons created in the APD is
116 // NumberOfElectrons = APDGain * LightYield
117 // The APD Gain is 300
118 fLightYieldMean = 47000;
119 fIntrinsicPINEfficiency = 0.02655 ; //APD= 0.1875/0.1271 * 0.018 (PIN)
120 fLightYieldAttenuation = 0.0045 ;
121 fRecalibrationFactor = 13.418/ fLightYieldMean ;
122 fElectronsPerGeV = 2.77e+8 ;
124 fLightFactor = fLightYieldMean * fIntrinsicPINEfficiency ;
125 fAPDFactor = (fRecalibrationFactor/100.) * fAPDGain ;
128 //____________________________________________________________________________
129 AliPHOSv1::~AliPHOSv1()
139 //____________________________________________________________________________
140 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t Id, Float_t * hits)
142 // Add a hit to the hit list.
143 // A PHOS hit is the sum of all hits in a single crystal from one primary and within some time gate
148 Bool_t deja = kFALSE ;
149 AliPHOSGeometry * geom = GetGeometry() ;
151 newHit = new AliPHOSHit(shunt, primary, Id, hits) ;
153 for ( hitCounter = fNhits-1 ; hitCounter >= 0 && !deja ; hitCounter-- ) {
154 curHit = dynamic_cast<AliPHOSHit*>((*fHits)[hitCounter]) ;
155 if(curHit->GetPrimary() != primary) break ;
156 // We add hits with the same primary, while GEANT treats primaries succesively
157 if( *curHit == *newHit ) {
164 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
165 // get the block Id number
167 geom->AbsToRelNumbering(Id, relid) ;
175 //____________________________________________________________________________
176 void AliPHOSv1::FinishPrimary()
178 // called at the end of each track (primary) by AliRun
179 // hits are reset for each new track
180 // accumulate the total hit-multiplicity
184 //____________________________________________________________________________
185 void AliPHOSv1::FinishEvent()
187 // called at the end of each event by AliRun
188 // accumulate the hit-multiplicity and total energy per block
189 // if the values have been updated check it
191 AliDetector::FinishEvent();
193 //____________________________________________________________________________
194 void AliPHOSv1::StepManager(void)
196 // Accumulates hits as long as the track stays in a single crystal or CPV gas Cell
198 Int_t relid[4] ; // (box, layer, row, column) indices
199 Int_t absid ; // absolute cell ID number
200 Float_t xyzte[5]={-1000.,-1000.,-1000.,0.,0.} ; // position wrt MRS, time and energy deposited
201 TLorentzVector pos ; // Lorentz vector of the track current position
204 TString name = GetGeometry()->GetName() ;
208 static Int_t idPCPQ = gMC->VolId("PCPQ");
209 if( gMC->CurrentVolID(copy) == idPCPQ &&
210 (gMC->IsTrackEntering() ) &&
211 gMC->TrackCharge() != 0) {
213 gMC -> TrackPosition(pos);
215 Float_t xyzm[3], xyzd[3] ;
217 for (i=0; i<3; i++) xyzm[i] = pos[i];
218 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
220 Float_t xyd[3]={0,0,0} ; //local position of the entering
225 // Current momentum of the hit's track in the local ref. system
226 TLorentzVector pmom ; //momentum of the particle initiated hit
227 gMC -> TrackMomentum(pmom);
228 Float_t pm[3], pd[3];
232 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
237 // Digitize the current CPV hit:
239 // 1. find pad response and
240 gMC->CurrentVolOffID(3,moduleNumber);
243 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
244 CPVDigitize(pmom,xyd,cpvDigits);
249 Int_t idigit,ndigits;
251 // 2. go through the current digit list and sum digits in pads
253 ndigits = cpvDigits->GetEntriesFast();
254 for (idigit=0; idigit<ndigits-1; idigit++) {
255 AliPHOSCPVDigit *cpvDigit1 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
256 Float_t x1 = cpvDigit1->GetXpad() ;
257 Float_t z1 = cpvDigit1->GetYpad() ;
258 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
259 AliPHOSCPVDigit *cpvDigit2 = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(jdigit));
260 Float_t x2 = cpvDigit2->GetXpad() ;
261 Float_t z2 = cpvDigit2->GetYpad() ;
262 if (x1==x2 && z1==z2) {
263 Float_t qsumpad = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
264 cpvDigit2->SetQpad(qsumpad) ;
265 cpvDigits->RemoveAt(idigit) ;
269 cpvDigits->Compress() ;
271 // 3. add digits to temporary hit list fTmpHits
273 ndigits = cpvDigits->GetEntriesFast();
274 for (idigit=0; idigit<ndigits; idigit++) {
275 AliPHOSCPVDigit *cpvDigit = dynamic_cast<AliPHOSCPVDigit*>(cpvDigits->UncheckedAt(idigit));
276 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
277 relid[1] =-1 ; // means CPV
278 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
279 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
281 // get the absolute Id number
282 GetGeometry()->RelToAbsNumbering(relid, absid) ;
284 // add current digit to the temporary hit list
286 xyzte[3] = gMC->TrackTime() ;
287 xyzte[4] = cpvDigit->GetQpad() ; // amplitude in a pad
289 Int_t primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
290 AddHit(fIshunt, primary, absid, xyzte);
292 if (cpvDigit->GetQpad() > 0.02) {
293 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
294 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
295 qsum += cpvDigit->GetQpad();
306 static Int_t idPXTL = gMC->VolId("PXTL");
307 if(gMC->CurrentVolID(copy) == idPXTL ) { // We are inside a PBWO crystal
309 gMC->TrackPosition(pos) ;
314 Float_t global[3], local[3] ;
318 Float_t lostenergy = gMC->Edep();
320 //Put in the TreeK particle entering PHOS and all its parents
321 if ( gMC->IsTrackEntering() ){
323 gMC -> Gmtod (xyzte, xyzd, 1); // transform coordinate from master to daughter system
324 if (xyzd[1] < -GetGeometry()->GetCrystalSize(1)/2.+0.1){ //Entered close to forward surface
325 Int_t parent = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
326 TParticle * part = gAlice->GetMCApp()->Particle(parent) ;
327 Float_t vert[3],vertd[3] ;
331 gMC -> Gmtod (vert, vertd, 1); // transform coordinate from master to daughter system
332 if(vertd[1]<-GetGeometry()->GetCrystalSize(1)/2.-0.1){ //Particle is created in foront of PHOS
333 //0.1 to get rid of numerical errors
334 part->SetBit(kKeepBit);
335 while ( parent != -1 ) {
336 part = gAlice->GetMCApp()->Particle(parent) ;
337 part->SetBit(kKeepBit);
338 parent = part->GetFirstMother() ;
343 if ( lostenergy != 0 ) { // Track is inside the crystal and deposits some energy
344 xyzte[3] = gMC->TrackTime() ;
346 gMC->CurrentVolOffID(10, moduleNumber) ; // get the PHOS module number ;
349 gMC->CurrentVolOffID(3, strip);
351 gMC->CurrentVolOffID(2, cell);
353 Int_t row = 1 + GetGeometry()->GetNZ() - strip % GetGeometry()->GetNZ() ;
354 Int_t col = (Int_t) TMath::Ceil((Double_t) strip/GetGeometry()->GetNZ()) -1 ;
356 absid = (moduleNumber-1)*GetGeometry()->GetNCristalsInModule() +
357 row + (col*GetGeometry()->GetEMCAGeometry()->GetNCellsInStrip() + cell-1)*GetGeometry()->GetNZ() ;
359 gMC->Gmtod(global, local, 1) ;
361 //Calculates the light yield, the number of photons produced in the
363 Float_t lightYield = gRandom->Poisson(fLightFactor * lostenergy *
364 exp(-fLightYieldAttenuation *
365 (local[1]+GetGeometry()->GetCrystalSize(1)/2.0 ))
368 //Calculates de energy deposited in the crystal
369 xyzte[4] = fAPDFactor * lightYield ;
373 primary = gAlice->GetMCApp()->GetCurrentTrackNumber() ;
374 TParticle * part = gAlice->GetMCApp()->Particle(primary) ;
375 while ( !part->TestBit(kKeepBit) ) {
376 primary = part->GetFirstMother() ;
378 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
379 break ; //there is a possibility that particle passed e.g. thermal isulator and hits a side
380 //surface of the crystal. In this case it may have no primary at all.
381 //We can not easily separate this case from the case when this is part of the shower,
382 //developed in the neighboring crystal.
384 part = gAlice->GetMCApp()->Particle(primary) ;
388 primary = gAlice->GetMCApp()->GetPrimary( gAlice->GetMCApp()->GetCurrentTrackNumber() );
392 // add current hit to the hit list
393 // Info("StepManager","%d %d", primary, tracknumber) ;
394 AddHit(fIshunt, primary, absid, xyzte);
396 } // there is deposited energy
397 } // we are inside a PHOS Xtal
401 //____________________________________________________________________________
402 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, TClonesArray *cpvDigits)
404 // ------------------------------------------------------------------------
405 // Digitize one CPV hit:
406 // On input take exact 4-momentum p and position zxhit of the hit,
407 // find the pad response around this hit and
408 // put the amplitudes in the pads into array digits
410 // Author: Yuri Kharlov (after Serguei Sadovsky)
412 // ------------------------------------------------------------------------
414 const Float_t kCelWr = GetGeometry()->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
415 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
416 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
417 const Int_t kNgamz = 5; // Ionization size in Z
418 const Int_t kNgamx = 9; // Ionization size in Phi
419 const Float_t kNoise = 0.03; // charge noise in one pad
423 // Just a reminder on axes notation in the CPV module:
424 // axis Z goes along the beam
425 // axis X goes across the beam in the module plane
426 // axis Y is a normal to the module plane showing from the IP
428 Float_t hitX = zxhit[0];
429 Float_t hitZ =-zxhit[1];
432 Float_t pNorm = p.Py();
433 Float_t eloss = kdEdx;
435 //Info("CPVDigitize", "YVK : %f %f | %f %f %d", hitX, hitZ, pX, pZ, pNorm) ;
437 Float_t dZY = pZ/pNorm * GetGeometry()->GetCPVGasThickness();
438 Float_t dXY = pX/pNorm * GetGeometry()->GetCPVGasThickness();
439 gRandom->Rannor(rnor1,rnor2);
440 eloss *= (1 + kDetR*rnor1) *
441 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(GetGeometry()->GetCPVGasThickness(),2)));
442 Float_t zhit1 = hitZ + GetGeometry()->GetCPVActiveSize(1)/2 - dZY/2;
443 Float_t xhit1 = hitX + GetGeometry()->GetCPVActiveSize(0)/2 - dXY/2;
444 Float_t zhit2 = zhit1 + dZY;
445 Float_t xhit2 = xhit1 + dXY;
447 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
448 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
452 if (iwht1==iwht2) { // incline 1-wire hit
454 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
455 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
457 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
458 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
461 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
463 Int_t iwht3 = (iwht1 + iwht2) / 2;
464 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
465 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
466 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
467 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
468 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
469 Float_t dxw1 = xhit1 - xwr13;
470 Float_t dxw2 = xhit2 - xwr23;
471 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
472 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
473 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
474 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
476 zxe[2][0] = eloss * egm1;
477 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
479 zxe[2][1] = eloss * egm2;
480 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
482 zxe[2][2] = eloss * egm3;
484 else { // incline 2-wire hit
486 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
487 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
488 Float_t xwr12 = (xwht1 + xwht2) / 2;
489 Float_t dxw1 = xhit1 - xwr12;
490 Float_t dxw2 = xhit2 - xwr12;
491 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
492 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
493 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
495 zxe[2][0] = eloss * egm1;
496 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
498 zxe[2][1] = eloss * egm2;
501 // Finite size of ionization region
503 Int_t nCellZ = GetGeometry()->GetNumberOfCPVPadsZ();
504 Int_t nCellX = GetGeometry()->GetNumberOfCPVPadsPhi();
505 Int_t nz3 = (kNgamz+1)/2;
506 Int_t nx3 = (kNgamx+1)/2;
507 cpvDigits->Expand(nIter*kNgamx*kNgamz);
508 TClonesArray &ldigits = *(static_cast<TClonesArray *>(cpvDigits));
510 for (Int_t iter=0; iter<nIter; iter++) {
512 Float_t zhit = zxe[0][iter];
513 Float_t xhit = zxe[1][iter];
514 Float_t qhit = zxe[2][iter];
515 Float_t zcell = zhit / GetGeometry()->GetPadSizeZ();
516 Float_t xcell = xhit / GetGeometry()->GetPadSizePhi();
517 if ( zcell<=0 || xcell<=0 ||
518 zcell>=nCellZ || xcell>=nCellX) return;
519 Int_t izcell = (Int_t) zcell;
520 Int_t ixcell = (Int_t) xcell;
521 Float_t zc = zcell - izcell - 0.5;
522 Float_t xc = xcell - ixcell - 0.5;
523 for (Int_t iz=1; iz<=kNgamz; iz++) {
524 Int_t kzg = izcell + iz - nz3;
525 if (kzg<=0 || kzg>nCellZ) continue;
526 Float_t zg = (Float_t)(iz-nz3) - zc;
527 for (Int_t ix=1; ix<=kNgamx; ix++) {
528 Int_t kxg = ixcell + ix - nx3;
529 if (kxg<=0 || kxg>nCellX) continue;
530 Float_t xg = (Float_t)(ix-nx3) - xc;
532 // Now calculate pad response
533 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
534 qpad += kNoise*rnor2;
535 if (qpad<0) continue;
537 // Fill the array with pad response ID and amplitude
538 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
544 //____________________________________________________________________________
545 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
546 // ------------------------------------------------------------------------
547 // Calculate the amplitude in one CPV pad using the
548 // cumulative pad response function
549 // Author: Yuri Kharlov (after Serguei Sadovski)
551 // ------------------------------------------------------------------------
553 Double_t dz = GetGeometry()->GetPadSizeZ() / 2;
554 Double_t dx = GetGeometry()->GetPadSizePhi() / 2;
555 Double_t z = zhit * GetGeometry()->GetPadSizeZ();
556 Double_t x = xhit * GetGeometry()->GetPadSizePhi();
557 Double_t amplitude = qhit *
558 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
559 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
560 return (Float_t)amplitude;
563 //____________________________________________________________________________
564 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
565 // ------------------------------------------------------------------------
566 // Cumulative pad response function
567 // It includes several terms from the CF decomposition in electrostatics
568 // Note: this cumulative function is wrong since omits some terms
569 // but the cell amplitude obtained with it is correct because
570 // these omitting terms cancel
571 // Author: Yuri Kharlov (after Serguei Sadovski)
573 // ------------------------------------------------------------------------
575 const Double_t kA=1.0;
576 const Double_t kB=0.7;
578 Double_t r2 = x*x + y*y;
580 Double_t cumulPRF = 0;
581 for (Int_t i=0; i<=4; i++) {
582 Double_t b1 = (2*i + 1) * kB;
583 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
585 cumulPRF *= kA/(2*TMath::Pi());