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 //_________________________________________________________________________
19 // Implementation version v1 of PHOS Manager class
21 // Layout EMC + PPSD has name GPS2:
22 // Produces cumulated hits (no hits) and digits
24 // Layout EMC + CPV has name IHEP:
25 // Produces hits for CPV, cumulated hits and digits
26 //*-- Author: Yves Schutz (SUBATECH)
29 // --- ROOT system ---
37 // --- Standard library ---
42 #include <strstream.h>
44 // --- AliRoot header files ---
46 #include "AliPHOSv1.h"
47 #include "AliPHOSHit.h"
48 #include "AliPHOSDigit.h"
49 #include "AliPHOSReconstructioner.h"
56 //____________________________________________________________________________
57 AliPHOSv1::AliPHOSv1()
63 // Create an empty array of AliPHOSCPVModule to satisfy
64 // AliPHOSv1::Streamer when reading root file
66 if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
67 Error("AliPHOSv1","Can not create array of CPV modules");
73 //____________________________________________________________________________
74 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
77 // ctor : title is used to identify the layout
78 // GPS2 = 5 modules (EMC + PPSD)
79 // IHEP = 5 modules (EMC + CPV )
80 // We use 2 arrays of hits :
82 // - fHits (the "normal" one), which retains the hits associated with
83 // the current primary particle being tracked
84 // (this array is reset after each primary has been tracked).
86 // - fTmpHits, which retains all the hits of the current event. It
87 // is used for the digitization part.
89 fPinElectronicNoise = 0.010 ;
90 fDigitThreshold = 0.1 ; // 1 GeV
92 // We do not want to save in TreeH the raw hits
93 // But save the cumulated hits instead (need to create the branch myself)
94 // It is put in the Digit Tree because the TreeH is filled after each primary
95 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
97 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
99 fNTmpHits = fNhits = 0 ;
101 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
103 fIshunt = 1 ; // All hits are associated with primary particles
105 // Create array of CPV modules for the IHEP's version of CPV
107 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
108 // Create array of AliPHOSCPVmodule of the size of PHOS modules number
110 if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
111 Error("AliPHOSv1","Can not create array of CPV modules");
114 TClonesArray &lcpvmodule = *fCPVModules;
115 for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
118 // Create an empty array of AliPHOSCPVModule to satisfy
119 // AliPHOSv1::Streamer when writing root file
121 fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
126 //____________________________________________________________________________
127 AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
128 AliPHOSv0(name,title)
130 // ctor : title is used to identify the layout
131 // GPS2 = 5 modules (EMC + PPSD)
132 // We use 2 arrays of hits :
134 // - fHits (the "normal" one), which retains the hits associated with
135 // the current primary particle being tracked
136 // (this array is reset after each primary has been tracked).
138 // - fTmpHits, which retains all the hits of the current event. It
139 // is used for the digitization part.
141 fPinElectronicNoise = 0.010 ;
143 // We do not want to save in TreeH the raw hits
144 //fHits = new TClonesArray("AliPHOSHit",100) ;
146 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
147 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
149 fNTmpHits = fNhits = 0 ;
151 fIshunt = 1 ; // All hits are associated with primary particles
153 // gets an instance of the geometry parameters class
154 fGeom = AliPHOSGeometry::GetInstance(title, "") ;
156 if (fGeom->IsInitialized() )
157 cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
159 cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
161 // Defining the PHOS Reconstructioner
163 fReconstructioner = Reconstructioner ;
167 //____________________________________________________________________________
168 AliPHOSv1::~AliPHOSv1()
178 // if ( fEmcRecPoints ) {
179 // fEmcRecPoints->Delete() ;
180 // delete fEmcRecPoints ;
181 // fEmcRecPoints = 0 ;
184 // if ( fPpsdRecPoints ) {
185 // fPpsdRecPoints->Delete() ;
186 // delete fPpsdRecPoints ;
187 // fPpsdRecPoints = 0 ;
190 // if ( fTrackSegments ) {
191 // fTrackSegments->Delete() ;
192 // delete fTrackSegments ;
193 // fTrackSegments = 0 ;
198 //____________________________________________________________________________
199 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
201 // Add a hit to the hit list.
202 // A PHOS hit is the sum of all hits in a single crystal
203 // or in a single PPSD gas cell
206 TClonesArray <mphits = *fTmpHits ;
209 Bool_t deja = kFALSE ;
211 // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
213 newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
215 // We do not want to save in TreeH the raw hits
216 // TClonesArray &lhits = *fHits;
218 for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
219 curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
220 if( *curHit == *newHit ) {
221 *curHit = *curHit + *newHit ;
227 new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
231 // We do not want to save in TreeH the raw hits
232 // new(lhits[fNhits]) AliPHOSHit(*newHit) ;
235 // Please note that the fTmpHits array must survive up to the
236 // end of the events, so it does not appear e.g. in ResetHits() (
237 // which is called at the end of each primary).
243 //___________________________________________________________________________
244 Int_t AliPHOSv1::Digitize(Float_t Energy)
246 // Applies the energy calibration
248 Float_t fB = 100000000. ;
250 Int_t chan = Int_t(fA + Energy*fB ) ;
254 //___________________________________________________________________________
255 void AliPHOSv1::FinishEvent()
257 // Makes the digits from the sum of summed hit in a single crystal or PPSD gas cell
258 // Adds to the energy the electronic noise
259 // Keeps digits with energy above fDigitThreshold
261 // Save the cumulated hits instead of raw hits (need to create the branch myself)
262 // It is put in the Digit Tree because the TreeH is filled after each primary
263 // and the TreeD at the end of the event.
269 TClonesArray &lDigits = *fDigits ;
271 AliPHOSDigit * newdigit ;
272 AliPHOSDigit * curdigit ;
273 Bool_t deja = kFALSE ;
275 for ( i = 0 ; i < fNTmpHits ; i++ ) {
276 hit = (AliPHOSHit*)fTmpHits->At(i) ;
278 // Assign primary number only if contribution is significant
279 if( hit->GetEnergy() > fDigitThreshold)
280 newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
282 newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
284 for ( j = 0 ; j < fNdigits ; j++) {
285 curdigit = (AliPHOSDigit*) lDigits[j] ;
286 if ( *curdigit == *newdigit) {
287 *curdigit = *curdigit + *newdigit ;
292 new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ;
299 // Noise induced by the PIN diode of the PbWO crystals
301 Float_t energyandnoise ;
302 for ( i = 0 ; i < fNdigits ; i++ ) {
303 newdigit = (AliPHOSDigit * ) fDigits->At(i) ;
305 fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
307 if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals)
308 energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
310 if (energyandnoise < 0 )
313 if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
314 fDigits->RemoveAt(i) ;
318 fDigits->Compress() ;
320 fNdigits = fDigits->GetEntries() ;
321 for (i = 0 ; i < fNdigits ; i++) {
322 newdigit = (AliPHOSDigit *) fDigits->At(i) ;
323 newdigit->SetIndexInList(i) ;
325 // fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
326 // printf("FinishEvent(): relid=(%d,%d,%d,%d) Amp=%d\n",
327 // relid[0],relid[1],relid[2],relid[3], newdigit->GetAmp());
332 //___________________________________________________________________________
333 void AliPHOSv1::MakeBranch(Option_t* opt)
335 // Create new branche in the current Root Tree in the digit Tree
336 AliDetector::MakeBranch(opt) ;
339 sprintf(branchname,"%s",GetName());
340 char *cdD = strstr(opt,"D");
341 if (fDigits && gAlice->TreeD() && cdD) {
342 gAlice->TreeD()->Branch(branchname, &fDigits, fBufferSize);
345 // Create new branche PHOSCH in the current Root Tree in the digit Tree for accumulated Hits
346 if ( ! (gAlice->IsLegoRun()) ) { // only when not in lego plot mode
347 if ( fTmpHits && gAlice->TreeD() && cdD) {
348 char branchname[10] ;
349 sprintf(branchname, "%sCH", GetName()) ;
350 gAlice->TreeD()->Branch(branchname, &fTmpHits, fBufferSize) ;
354 // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
355 // Yuri Kharlov, 28 September 2000.
357 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
358 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).MakeBranch(i+1);
363 //_____________________________________________________________________________
364 void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
366 // 1. Reinitializes the existing RecPoint, TrackSegment, and RecParticles Lists and
367 // 2. Creates TreeR with a branch for each list
368 // 3. Steers the reconstruction processes
369 // 4. Saves the 3 lists in TreeR
370 // 5. Write the Tree to File
372 fReconstructioner = Reconstructioner ;
374 char branchname[10] ;
378 // gAlice->MakeTree("R") ;
379 Int_t splitlevel = 0 ;
381 fEmcRecPoints->Delete() ;
383 if ( fEmcRecPoints && gAlice->TreeR() ) {
384 sprintf(branchname,"%sEmcRP",GetName()) ;
385 gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
388 fPpsdRecPoints->Delete() ;
390 if ( fPpsdRecPoints && gAlice->TreeR() ) {
391 sprintf(branchname,"%sPpsdRP",GetName()) ;
392 gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
395 fTrackSegments->Delete() ;
397 if ( fTrackSegments && gAlice->TreeR() ) {
398 sprintf(branchname,"%sTS",GetName()) ;
399 gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
402 fRecParticles->Delete() ;
404 if (strcmp(fGeom->GetName(),"GPS2") == 0) {
405 if ( fRecParticles && gAlice->TreeR() ) {
406 sprintf(branchname,"%sRP",GetName()) ;
407 gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
412 if (strcmp(fGeom->GetName(),"GPS2") == 0)
413 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
414 else if (strcmp(fGeom->GetName(),"IHEP") == 0)
415 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints);
417 // 4. Expand or Shrink the arrays to the proper size
421 size = fEmcRecPoints->GetEntries() ;
422 fEmcRecPoints->Expand(size) ;
424 size = fPpsdRecPoints->GetEntries() ;
425 fPpsdRecPoints->Expand(size) ;
427 size = fTrackSegments->GetEntries() ;
428 fTrackSegments->Expand(size) ;
430 size = fRecParticles->GetEntries() ;
431 fRecParticles->Expand(size) ;
433 gAlice->TreeR()->Fill() ;
436 gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
438 // Deleting reconstructed objects
439 ResetReconstruction();
443 //____________________________________________________________________________
444 void AliPHOSv1::ResetHits()
446 // Reset hit tree for CPV in IHEP geometry
447 // Yuri Kharlov, 28 September 2000
449 AliDetector::ResetHits();
450 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
451 for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
454 //____________________________________________________________________________
455 void AliPHOSv1::ResetDigits()
457 // May sound strange, but cumulative hits are store in digits Tree
458 AliDetector::ResetDigits();
464 //____________________________________________________________________________
465 void AliPHOSv1::ResetReconstruction()
467 // Deleting reconstructed objects
469 if ( fEmcRecPoints ) fEmcRecPoints->Delete();
470 if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
471 if ( fTrackSegments ) fTrackSegments->Delete();
472 if ( fRecParticles ) fRecParticles->Delete();
476 //____________________________________________________________________________
477 void AliPHOSv1::SetTreeAddress()
480 AliPHOS::SetTreeAddress();
482 // //Branch address for TreeR: RecPpsdRecPoint
483 // TTree *treeR = gAlice->TreeR();
484 // if ( treeR && fPpsdRecPoints ) {
485 // branch = treeR->GetBranch("PHOSPpsdRP");
486 // if (branch) branch->SetAddress(&fPpsdRecPoints) ;
489 // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
490 // Yuri Kharlov, 28 September 2000.
492 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
493 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).SetTreeAddress(i+1);
498 //____________________________________________________________________________
500 void AliPHOSv1::StepManager(void)
502 // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
504 Int_t relid[4] ; // (box, layer, row, column) indices
505 Int_t absid ; // absolute cell ID number
506 Float_t xyze[4] ; // position wrt MRS and energy deposited
507 TLorentzVector pos ; // Lorentz vector of the track current position
510 Int_t tracknumber = gAlice->CurrentTrack() ;
511 Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
512 TString name = fGeom->GetName() ;
514 if ( name == "GPS2" ) { // ======> CPV is a GPS' PPSD
516 if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
518 gMC->TrackPosition(pos) ;
522 xyze[3] = gMC->Edep() ;
524 if ( xyze[3] != 0 ) { // there is deposited energy
525 gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number
526 gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
527 // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
528 // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
529 gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
530 gMC->CurrentVolID(relid[3]) ; // get the column number
532 // get the absolute Id number
534 fGeom->RelToAbsNumbering(relid, absid) ;
536 // add current hit to the hit list
537 AddHit(fIshunt, primary, tracknumber, absid, xyze);
539 } // there is deposited energy
540 } // We are inside the gas of the CPV
541 } // GPS2 configuration
543 else if ( name == "IHEP" ) { // ======> CPV is a IHEP's one
545 // Yuri Kharlov, 28 September 2000
547 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
548 gMC->IsTrackEntering() &&
549 gMC->TrackCharge() != 0) {
551 // Charged track has just entered to the CPV sensitive plane
553 AliPHOSv1 &phos = *(AliPHOSv1*)gAlice->GetModule("PHOS");
556 gMC->CurrentVolOffID(3,moduleNumber);
559 // Current position of the hit in the CPV module ref. system
561 gMC -> TrackPosition(pos);
562 Float_t xyzm[3], xyzd[3], xyd[2];
564 for (i=0; i<3; i++) xyzm[i] = pos[i];
565 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
569 // Current momentum of the hit's track in the CPV module ref. system
572 gMC -> TrackMomentum(pmom);
573 Float_t pm[3], pd[3];
574 for (i=0; i<3; i++) pm[i] = pmom[i];
575 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
580 // Current particle type of the hit's track
582 Int_t ipart = gMC->TrackPid();
584 // Add the current particle in the list of the CPV hits.
586 phos.GetCPVModule(moduleNumber).AddHit(pmom,xyd,ipart);
588 if (fDebugLevel == 1) {
589 printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
590 moduleNumber+1,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
591 printf( " xy = (%8.4f, %8.4f) cm, ipart = %d\n",
592 xyd[0],xyd[1],ipart);
595 // Digitize the current CPV hit:
597 // 1. find pad response and
599 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
600 CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
605 Int_t idigit,ndigits;
607 // 2. go through the current digit list and sum digits in pads
609 ndigits = cpvDigits->GetEntriesFast();
610 for (idigit=0; idigit<ndigits-1; idigit++) {
611 AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
612 Float_t x1 = cpvDigit1->GetXpad() ;
613 Float_t z1 = cpvDigit1->GetYpad() ;
614 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
615 AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
616 Float_t x2 = cpvDigit2->GetXpad() ;
617 Float_t z2 = cpvDigit2->GetYpad() ;
618 if (x1==x2 && z1==z2) {
619 Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
620 cpvDigit2->SetQpad(qsum) ;
621 cpvDigits->RemoveAt(idigit) ;
625 cpvDigits->Compress() ;
627 // 3. add digits to temporary hit list fTmpHits
629 ndigits = cpvDigits->GetEntriesFast();
630 for (idigit=0; idigit<ndigits; idigit++) {
631 AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
632 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
633 relid[1] =-1 ; // means CPV
634 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
635 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
637 // get the absolute Id number
638 fGeom->RelToAbsNumbering(relid, absid) ;
640 // add current digit to the temporary hit list
644 xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
645 primary = -1; // No need in primary for CPV
646 AddHit(fIshunt, primary, tracknumber, absid, xyze);
648 if (cpvDigit->GetQpad() > 0.02) {
649 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
650 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
651 qsum += cpvDigit->GetQpad();
656 } // end of IHEP configuration
658 if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
659 gMC->TrackPosition(pos) ;
663 xyze[3] = gMC->Edep() ;
665 if ( xyze[3] != 0 ) {
666 gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
667 relid[1] = 0 ; // means PBW04
668 gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
669 gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
671 // get the absolute Id number
673 fGeom->RelToAbsNumbering(relid, absid) ;
675 // add current hit to the hit list
677 AddHit(fIshunt, primary,tracknumber, absid, xyze);
679 } // there is deposited energy
680 } // we are inside a PHOS Xtal
683 //____________________________________________________________________________
684 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
686 // ------------------------------------------------------------------------
687 // Digitize one CPV hit:
688 // On input take exact 4-momentum p and position zxhit of the hit,
689 // find the pad response around this hit and
690 // put the amplitudes in the pads into array digits
692 // Author: Yuri Kharlov (after Serguei Sadovsky)
694 // ------------------------------------------------------------------------
696 const Float_t kCelWr = fGeom->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
697 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
698 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
699 const Int_t kNgamz = 5; // Ionization size in Z
700 const Int_t kNgamx = 9; // Ionization size in Phi
701 const Float_t kNoise = 0.03; // charge noise in one pad
705 // Just a reminder on axes notation in the CPV module:
706 // axis Z goes along the beam
707 // axis X goes across the beam in the module plane
708 // axis Y is a normal to the module plane showing from the IP
710 Float_t hitX = zxhit[0];
711 Float_t hitZ =-zxhit[1];
714 Float_t pNorm = p.Py();
715 Float_t eloss = kdEdx;
717 Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
718 Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
719 gRandom->Rannor(rnor1,rnor2);
720 eloss *= (1 + kDetR*rnor1) *
721 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
722 Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
723 Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
724 Float_t zhit2 = zhit1 + dZY;
725 Float_t xhit2 = xhit1 + dXY;
727 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
728 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
732 if (iwht1==iwht2) { // incline 1-wire hit
734 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
735 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
737 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
738 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
741 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
743 Int_t iwht3 = (iwht1 + iwht2) / 2;
744 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
745 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
746 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
747 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
748 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
749 Float_t dxw1 = xhit1 - xwr13;
750 Float_t dxw2 = xhit2 - xwr23;
751 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
752 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
753 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
754 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
756 zxe[2][0] = eloss * egm1;
757 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
759 zxe[2][1] = eloss * egm2;
760 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
762 zxe[2][2] = eloss * egm3;
764 else { // incline 2-wire hit
766 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
767 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
768 Float_t xwr12 = (xwht1 + xwht2) / 2;
769 Float_t dxw1 = xhit1 - xwr12;
770 Float_t dxw2 = xhit2 - xwr12;
771 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
772 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
773 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
775 zxe[2][0] = eloss * egm1;
776 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
778 zxe[2][1] = eloss * egm2;
781 // Finite size of ionization region
783 Int_t nCellZ = fGeom->GetNumberOfPadsZ();
784 Int_t nCellX = fGeom->GetNumberOfPadsPhi();
785 Int_t nz3 = (kNgamz+1)/2;
786 Int_t nx3 = (kNgamx+1)/2;
787 cpvDigits->Expand(nIter*kNgamx*kNgamz);
788 TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
790 for (Int_t iter=0; iter<nIter; iter++) {
792 Float_t zhit = zxe[0][iter];
793 Float_t xhit = zxe[1][iter];
794 Float_t qhit = zxe[2][iter];
795 Float_t zcell = zhit / fGeom->GetPadSizeZ();
796 Float_t xcell = xhit / fGeom->GetPadSizePhi();
797 if ( zcell<=0 || xcell<=0 ||
798 zcell>=nCellZ || xcell>=nCellX) return;
799 Int_t izcell = (Int_t) zcell;
800 Int_t ixcell = (Int_t) xcell;
801 Float_t zc = zcell - izcell - 0.5;
802 Float_t xc = xcell - ixcell - 0.5;
803 for (Int_t iz=1; iz<=kNgamz; iz++) {
804 Int_t kzg = izcell + iz - nz3;
805 if (kzg<=0 || kzg>nCellZ) continue;
806 Float_t zg = (Float_t)(iz-nz3) - zc;
807 for (Int_t ix=1; ix<=kNgamx; ix++) {
808 Int_t kxg = ixcell + ix - nx3;
809 if (kxg<=0 || kxg>nCellX) continue;
810 Float_t xg = (Float_t)(ix-nx3) - xc;
812 // Now calculate pad response
813 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
814 qpad += kNoise*rnor2;
815 if (qpad<0) continue;
817 // Fill the array with pad response ID and amplitude
818 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
824 //____________________________________________________________________________
825 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
826 // ------------------------------------------------------------------------
827 // Calculate the amplitude in one CPV pad using the
828 // cumulative pad response function
829 // Author: Yuri Kharlov (after Serguei Sadovski)
831 // ------------------------------------------------------------------------
833 Double_t dz = fGeom->GetPadSizeZ() / 2;
834 Double_t dx = fGeom->GetPadSizePhi() / 2;
835 Double_t z = zhit * fGeom->GetPadSizeZ();
836 Double_t x = xhit * fGeom->GetPadSizePhi();
837 Double_t amplitude = qhit *
838 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
839 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
840 return (Float_t)amplitude;
843 //____________________________________________________________________________
844 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
845 // ------------------------------------------------------------------------
846 // Cumulative pad response function
847 // It includes several terms from the CF decomposition in electrostatics
848 // Note: this cumulative function is wrong since omits some terms
849 // but the cell amplitude obtained with it is correct because
850 // these omitting terms cancel
851 // Author: Yuri Kharlov (after Serguei Sadovski)
853 // ------------------------------------------------------------------------
855 const Double_t kA=1.0;
856 const Double_t kB=0.7;
858 Double_t r2 = x*x + y*y;
860 Double_t cumulPRF = 0;
861 for (Int_t i=0; i<=4; i++) {
862 Double_t b1 = (2*i + 1) * kB;
863 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
865 cumulPRF *= kA/(2*TMath::Pi());