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
20 // Layout EMC + PPSD has name GPS2
21 // Produces cumulated hits (no hits) and digits
22 //*-- Author: Yves Schutz (SUBATECH)
25 // --- ROOT system ---
33 // --- Standard library ---
38 #include <strstream.h>
40 // --- AliRoot header files ---
42 #include "AliPHOSv1.h"
43 #include "AliPHOSHit.h"
44 #include "AliPHOSDigit.h"
45 #include "AliPHOSReconstructioner.h"
52 //____________________________________________________________________________
53 AliPHOSv1::AliPHOSv1()
59 // Create an empty array of CPVModule to satisfy
60 // AliPHOSv1::Streamer when reading root file
62 if ( NULL==(fCPVModules=new TClonesArray("CPVModule",0)) ) {
63 Error("AliPHOSv1","Can not create array of CPV modules");
69 //____________________________________________________________________________
70 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
73 // ctor : title is used to identify the layout
74 // GPS2 = 5 modules (EMC + PPSD)
75 // IHEP = 5 modules (EMC + CPV )
76 // We use 2 arrays of hits :
78 // - fHits (the "normal" one), which retains the hits associated with
79 // the current primary particle being tracked
80 // (this array is reset after each primary has been tracked).
82 // - fTmpHits, which retains all the hits of the current event. It
83 // is used for the digitization part.
85 fPinElectronicNoise = 0.010 ;
86 fDigitThreshold = 0.1 ; // 1 GeV
88 // We do not want to save in TreeH the raw hits
89 // But save the cumulated hits instead (need to create the branch myself)
90 // It is put in the Digit Tree because the TreeH is filled after each primary
91 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
93 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
95 fNTmpHits = fNhits = 0 ;
97 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
99 fIshunt = 1 ; // All hits are associated with primary particles
101 // Create array of CPV modules for the IHEP's version of CPV
103 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
104 // Create array of CPVmodule of the size of PHOS modules number
106 if ( NULL==(fCPVModules=new TClonesArray("CPVModule",fGeom->GetNModules())) ) {
107 Error("AliPHOSv1","Can not create array of CPV modules");
110 TClonesArray &lcpvmodule = *fCPVModules;
111 for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lcpvmodule[i]) CPVModule();
114 // Create an empty array of CPVModule to satisfy
115 // AliPHOSv1::Streamer when writing root file
117 fCPVModules=new TClonesArray("CPVModule",0);
122 //____________________________________________________________________________
123 AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
124 AliPHOSv0(name,title)
126 // ctor : title is used to identify the layout
127 // GPS2 = 5 modules (EMC + PPSD)
128 // We use 2 arrays of hits :
130 // - fHits (the "normal" one), which retains the hits associated with
131 // the current primary particle being tracked
132 // (this array is reset after each primary has been tracked).
134 // - fTmpHits, which retains all the hits of the current event. It
135 // is used for the digitization part.
137 fPinElectronicNoise = 0.010 ;
139 // We do not want to save in TreeH the raw hits
140 //fHits = new TClonesArray("AliPHOSHit",100) ;
142 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
143 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
145 fNTmpHits = fNhits = 0 ;
147 fIshunt = 1 ; // All hits are associated with primary particles
149 // gets an instance of the geometry parameters class
150 fGeom = AliPHOSGeometry::GetInstance(title, "") ;
152 if (fGeom->IsInitialized() )
153 cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
155 cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
157 // Defining the PHOS Reconstructioner
159 fReconstructioner = Reconstructioner ;
163 //____________________________________________________________________________
164 AliPHOSv1::~AliPHOSv1()
174 // if ( fEmcRecPoints ) {
175 // fEmcRecPoints->Delete() ;
176 // delete fEmcRecPoints ;
177 // fEmcRecPoints = 0 ;
180 // if ( fPpsdRecPoints ) {
181 // fPpsdRecPoints->Delete() ;
182 // delete fPpsdRecPoints ;
183 // fPpsdRecPoints = 0 ;
186 // if ( fTrackSegments ) {
187 // fTrackSegments->Delete() ;
188 // delete fTrackSegments ;
189 // fTrackSegments = 0 ;
194 //____________________________________________________________________________
195 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
197 // Add a hit to the hit list.
198 // A PHOS hit is the sum of all hits in a single crystal
199 // or in a single PPSD gas cell
202 TClonesArray <mphits = *fTmpHits ;
205 Bool_t deja = kFALSE ;
207 // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
209 newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
211 // We do not want to save in TreeH the raw hits
212 // TClonesArray &lhits = *fHits;
214 for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
215 curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
216 if( *curHit == *newHit ) {
217 *curHit = *curHit + *newHit ;
223 new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
227 // We do not want to save in TreeH the raw hits
228 // new(lhits[fNhits]) AliPHOSHit(*newHit) ;
231 // Please note that the fTmpHits array must survive up to the
232 // end of the events, so it does not appear e.g. in ResetHits() (
233 // which is called at the end of each primary).
239 //___________________________________________________________________________
240 Int_t AliPHOSv1::Digitize(Float_t Energy)
242 // Applies the energy calibration
244 Float_t fB = 100000000. ;
246 Int_t chan = Int_t(fA + Energy*fB ) ;
250 //___________________________________________________________________________
251 void AliPHOSv1::FinishEvent()
253 // Makes the digits from the sum of summed hit in a single crystal or PPSD gas cell
254 // Adds to the energy the electronic noise
255 // Keeps digits with energy above fDigitThreshold
257 // Save the cumulated hits instead of raw hits (need to create the branch myself)
258 // It is put in the Digit Tree because the TreeH is filled after each primary
259 // and the TreeD at the end of the event.
265 TClonesArray &lDigits = *fDigits ;
267 AliPHOSDigit * newdigit ;
268 AliPHOSDigit * curdigit ;
269 Bool_t deja = kFALSE ;
271 for ( i = 0 ; i < fNTmpHits ; i++ ) {
272 hit = (AliPHOSHit*)fTmpHits->At(i) ;
274 // Assign primary number only if contribution is significant
275 if( hit->GetEnergy() > fDigitThreshold)
276 newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
278 newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
280 for ( j = 0 ; j < fNdigits ; j++) {
281 curdigit = (AliPHOSDigit*) lDigits[j] ;
282 if ( *curdigit == *newdigit) {
283 *curdigit = *curdigit + *newdigit ;
288 new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ;
295 // Noise induced by the PIN diode of the PbWO crystals
297 Float_t energyandnoise ;
298 for ( i = 0 ; i < fNdigits ; i++ ) {
299 newdigit = (AliPHOSDigit * ) fDigits->At(i) ;
301 fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
303 if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals)
304 energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
306 if (energyandnoise < 0 )
309 if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
310 fDigits->RemoveAt(i) ;
314 fDigits->Compress() ;
316 fNdigits = fDigits->GetEntries() ;
317 for (i = 0 ; i < fNdigits ; i++) {
318 newdigit = (AliPHOSDigit *) fDigits->At(i) ;
319 newdigit->SetIndexInList(i) ;
321 // fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
322 // printf("FinishEvent(): relid=(%d,%d,%d,%d) Amp=%d\n",
323 // relid[0],relid[1],relid[2],relid[3], newdigit->GetAmp());
328 //___________________________________________________________________________
329 void AliPHOSv1::MakeBranch(Option_t* opt)
331 // Create new branche in the current Root Tree in the digit Tree
332 AliDetector::MakeBranch(opt) ;
335 sprintf(branchname,"%s",GetName());
336 char *cdD = strstr(opt,"D");
337 if (fDigits && gAlice->TreeD() && cdD) {
338 gAlice->TreeD()->Branch(branchname, &fDigits, fBufferSize);
341 // Create new branche PHOSCH in the current Root Tree in the digit Tree for accumulated Hits
342 if ( ! (gAlice->IsLegoRun()) ) { // only when not in lego plot mode
343 if ( fTmpHits && gAlice->TreeD() && cdD) {
344 char branchname[10] ;
345 sprintf(branchname, "%sCH", GetName()) ;
346 gAlice->TreeD()->Branch(branchname, &fTmpHits, fBufferSize) ;
350 // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
351 // Yuri Kharlov, 28 September 2000.
353 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
354 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).MakeBranch(i+1);
359 //_____________________________________________________________________________
360 void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
362 // 1. Reinitializes the existing RecPoint, TrackSegment, and RecParticles Lists and
363 // 2. Creates TreeR with a branch for each list
364 // 3. Steers the reconstruction processes
365 // 4. Saves the 3 lists in TreeR
366 // 5. Write the Tree to File
368 fReconstructioner = Reconstructioner ;
370 char branchname[10] ;
374 // gAlice->MakeTree("R") ;
375 Int_t splitlevel = 0 ;
377 fEmcRecPoints->Delete() ;
379 if ( fEmcRecPoints && gAlice->TreeR() ) {
380 sprintf(branchname,"%sEmcRP",GetName()) ;
382 gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
385 fPpsdRecPoints->Delete() ;
388 if ( fPpsdRecPoints && gAlice->TreeR() ) {
389 sprintf(branchname,"%sPpsdRP",GetName()) ;
391 gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
394 fTrackSegments->Delete() ;
396 if ( fTrackSegments && gAlice->TreeR() ) {
397 sprintf(branchname,"%sTS",GetName()) ;
398 gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
401 fRecParticles->Delete() ;
403 if (strcmp(fGeom->GetName(),"GPS2") == 0) {
404 if ( fRecParticles && gAlice->TreeR() ) {
405 sprintf(branchname,"%sRP",GetName()) ;
406 gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
411 if (strcmp(fGeom->GetName(),"GPS2") == 0)
412 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
413 else if (strcmp(fGeom->GetName(),"IHEP") == 0)
414 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints);
416 // 4. Expand or Shrink the arrays to the proper size
420 size = fEmcRecPoints->GetEntries() ;
421 fEmcRecPoints->Expand(size) ;
423 size = fPpsdRecPoints->GetEntries() ;
424 fPpsdRecPoints->Expand(size) ;
426 size = fTrackSegments->GetEntries() ;
427 fTrackSegments->Expand(size) ;
429 size = fRecParticles->GetEntries() ;
430 fRecParticles->Expand(size) ;
432 gAlice->TreeR()->Fill() ;
435 gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
437 // Deleting reconstructed objects
438 ResetReconstruction();
442 //____________________________________________________________________________
443 void AliPHOSv1::ResetHits()
445 // Reset hit tree for CPV in IHEP geometry
446 // Yuri Kharlov, 28 September 2000
448 AliDetector::ResetHits();
449 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
450 for (Int_t i=0; i<fGeom->GetNModules(); i++) ((CPVModule*)(*fCPVModules)[i]) -> Clear();
453 //____________________________________________________________________________
454 void AliPHOSv1::ResetDigits()
456 // May sound strange, but cumulative hits are store in digits Tree
457 AliDetector::ResetDigits();
463 //____________________________________________________________________________
464 void AliPHOSv1::ResetReconstruction()
466 // Deleting reconstructed objects
468 if ( fEmcRecPoints ) fEmcRecPoints->Delete();
469 if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
470 if ( fTrackSegments ) fTrackSegments->Delete();
471 if ( fRecParticles ) fRecParticles->Delete();
475 //____________________________________________________________________________
476 void AliPHOSv1::SetTreeAddress()
479 AliPHOS::SetTreeAddress();
481 // //Branch address for TreeR: RecPpsdRecPoint
482 // TTree *treeR = gAlice->TreeR();
483 // if ( treeR && fPpsdRecPoints ) {
484 // branch = treeR->GetBranch("PHOSPpsdRP");
485 // if (branch) branch->SetAddress(&fPpsdRecPoints) ;
488 // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
489 // Yuri Kharlov, 28 September 2000.
491 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
492 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).SetTreeAddress(i+1);
497 //____________________________________________________________________________
499 void AliPHOSv1::StepManager(void)
501 // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
503 Int_t relid[4] ; // (box, layer, row, column) indices
504 Int_t absid ; // absolute cell ID number
505 Float_t xyze[4] ; // position wrt MRS and energy deposited
506 TLorentzVector pos ; // Lorentz vector of the track current position
509 Int_t tracknumber = gAlice->CurrentTrack() ;
510 Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
511 TString name = fGeom->GetName() ;
513 if ( name == "GPS2" ) { // ======> CPV is a GPS' PPSD
515 if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
517 gMC->TrackPosition(pos) ;
521 xyze[3] = gMC->Edep() ;
523 if ( xyze[3] != 0 ) { // there is deposited energy
524 gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number
525 gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
526 // 1-> Geom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
527 // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
528 gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
529 gMC->CurrentVolID(relid[3]) ; // get the column number
531 // get the absolute Id number
533 fGeom->RelToAbsNumbering(relid, absid) ;
535 // add current hit to the hit list
536 AddHit(fIshunt, primary, tracknumber, absid, xyze);
538 } // there is deposited energy
539 } // We are inside the gas of the CPV
540 } // GPS2 configuration
542 else if ( name == "IHEP" ) { // ======> CPV is a IHEP's one
544 // Yuri Kharlov, 28 September 2000
546 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
547 gMC->IsTrackEntering() &&
548 gMC->TrackCharge() != 0) {
550 // Charged track has just entered to the CPV sensitive plane
552 AliPHOSv1 &PHOS = *(AliPHOSv1*)gAlice->GetModule("PHOS");
555 gMC->CurrentVolOffID(3,ModuleNumber);
558 // Current position of the hit in the CPV module ref. system
560 gMC -> TrackPosition(pos);
561 Float_t xyzm[3], xyzd[3], xyd[2];
562 for (Int_t i=0; i<3; i++) xyzm[i] = pos[i];
563 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
567 // Current momentum of the hit's track in the CPV module ref. system
570 gMC -> TrackMomentum(pmom);
571 Float_t pm[3], pd[3];
572 for (Int_t i=0; i<3; i++) pm[i] = pmom[i];
573 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
578 // Current particle type of the hit's track
580 Int_t ipart = gMC->TrackPid();
582 // Add the current particle in the list of the CPV hits.
584 PHOS.GetCPVModule(ModuleNumber).AddHit(pmom,xyd,ipart);
586 printf("CPV hit added to module #%2d: p = (% .4f, % .4f, % .4f, % .4f) GeV,\n",
587 ModuleNumber,pmom.Px(),pmom.Py(),pmom.Pz(),pmom.E());
588 printf( " xy = (%8.4f, %8.4f) cm, ipart = %d\n",
589 xyd[0],xyd[1],ipart);
591 // Digitize the current CPV hit:
593 // 1. find pad response and
595 TClonesArray *cpvDigits = new TClonesArray("CPVDigit",0); // array of digits for current hit
596 CPVDigitize(pmom,xyd,ModuleNumber,cpvDigits);
603 // 2. go through the current digit list and sum digits in pads
605 ndigits = cpvDigits->GetEntriesFast();
606 for (Int_t idigit=0; idigit<ndigits-1; idigit++) {
607 CPVDigit *cpvDigit1 = (CPVDigit*) cpvDigits->UncheckedAt(idigit);
608 Float_t x1 = cpvDigit1->GetXpad() ;
609 Float_t z1 = cpvDigit1->GetYpad() ;
610 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
611 CPVDigit *cpvDigit2 = (CPVDigit*) cpvDigits->UncheckedAt(jdigit);
612 Float_t x2 = cpvDigit2->GetXpad() ;
613 Float_t z2 = cpvDigit2->GetYpad() ;
614 if (x1==x2 && z1==z2) {
615 Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
616 cpvDigit2->SetQpad(qsum) ;
617 cpvDigits->RemoveAt(idigit) ;
621 cpvDigits->Compress() ;
623 // 3. add digits to temporary hit list fTmpHits
625 ndigits = cpvDigits->GetEntriesFast();
626 for (Int_t idigit=0; idigit<ndigits; idigit++) {
627 CPVDigit *cpvDigit = (CPVDigit*) cpvDigits->UncheckedAt(idigit);
628 relid[0] = ModuleNumber + 1 ; // CPV (or PHOS) module number
629 relid[1] =-1 ; // means CPV
630 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
631 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
633 // get the absolute Id number
634 fGeom->RelToAbsNumbering(relid, absid) ;
636 // add current digit to the temporary hit list
640 xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
641 primary = -1; // No need in primary for CPV
642 AddHit(fIshunt, primary, tracknumber, absid, xyze);
644 if (cpvDigit->GetQpad() > 0.02) {
645 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
646 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
647 qsum += cpvDigit->GetQpad();
652 } // end of IHEP configuration
654 if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
655 gMC->TrackPosition(pos) ;
659 xyze[3] = gMC->Edep() ;
661 if ( xyze[3] != 0 ) {
662 gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
663 relid[1] = 0 ; // means PBW04
664 gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
665 gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
667 // get the absolute Id number
669 fGeom->RelToAbsNumbering(relid, absid) ;
671 // add current hit to the hit list
673 AddHit(fIshunt, primary,tracknumber, absid, xyze);
675 } // there is deposited energy
676 } // we are inside a PHOS Xtal
679 //____________________________________________________________________________
680 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
682 // ------------------------------------------------------------------------
683 // Digitize one CPV hit:
684 // On input take exact 4-momentum p and position zxhit of the hit,
685 // find the pad response around this hit and
686 // put the amplitudes in the pads into array digits
688 // Author: Yuri Kharlov (after Serguei Sadovsky)
690 // ------------------------------------------------------------------------
692 const Float_t celWr = fGeom->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
693 const Float_t detR = 0.1; // Relative energy fluctuation in track for 100 e-
694 const Float_t dEdx = 4.0; // Average energy loss in CPV;
695 const Int_t ngamz = 5; // Ionization size in Z
696 const Int_t ngamx = 9; // Ionization size in Phi
697 const Float_t qNoise = 0.03; // charge noise in one pad
701 // Just a reminder on axes notation in the CPV module:
702 // axis Z goes along the beam
703 // axis X goes across the beam in the module plane
704 // axis Y is a normal to the module plane showing from the IP
706 Float_t hitX = zxhit[0];
707 Float_t hitZ =-zxhit[1];
710 Float_t pNorm = p.Py();
713 Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
714 Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
715 gRandom->Rannor(rnor1,rnor2);
716 E *= (1 + detR*rnor1) *
717 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
718 Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
719 Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
720 Float_t zhit2 = zhit1 + dZY;
721 Float_t xhit2 = xhit1 + dXY;
723 Int_t iwht1 = (Int_t) (xhit1 / celWr); // wire (x) coordinate "in"
724 Int_t iwht2 = (Int_t) (xhit2 / celWr); // wire (x) coordinate "out"
728 if (iwht1==iwht2) { // incline 1-wire hit
730 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
731 zxe[1][0] = (iwht1 + 0.5) * celWr;
733 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
734 zxe[1][1] = (iwht1 + 0.5) * celWr;
737 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
739 Int_t iwht3 = (iwht1 + iwht2) / 2;
740 Float_t xwht1 = (iwht1 + 0.5) * celWr; // wire 1
741 Float_t xwht2 = (iwht2 + 0.5) * celWr; // wire 2
742 Float_t xwht3 = (iwht3 + 0.5) * celWr; // wire 3
743 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
744 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
745 Float_t dxw1 = xhit1 - xwr13;
746 Float_t dxw2 = xhit2 - xwr23;
747 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
748 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
749 Float_t egm3 = celWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
750 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
752 zxe[2][0] = E * egm1;
753 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
755 zxe[2][1] = E * egm2;
756 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
758 zxe[2][2] = E * egm3;
760 else { // incline 2-wire hit
762 Float_t xwht1 = (iwht1 + 0.5) * celWr;
763 Float_t xwht2 = (iwht2 + 0.5) * celWr;
764 Float_t xwr12 = (xwht1 + xwht2) / 2;
765 Float_t dxw1 = xhit1 - xwr12;
766 Float_t dxw2 = xhit2 - xwr12;
767 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
768 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
769 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
771 zxe[2][0] = E * egm1;
772 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
774 zxe[2][1] = E * egm2;
777 // Finite size of ionization region
779 Int_t nCellZ = fGeom->GetNumberOfPadsZ();
780 Int_t nCellX = fGeom->GetNumberOfPadsPhi();
781 Int_t nz3 = (ngamz+1)/2;
782 Int_t nx3 = (ngamx+1)/2;
783 cpvDigits->Expand(nIter*ngamx*ngamz);
784 TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
786 for (Int_t iter=0; iter<nIter; iter++) {
788 Float_t zhit = zxe[0][iter];
789 Float_t xhit = zxe[1][iter];
790 Float_t qhit = zxe[2][iter];
791 Float_t zcell = zhit / fGeom->GetPadSizeZ();
792 Float_t xcell = xhit / fGeom->GetPadSizePhi();
793 if ( zcell<=0 || xcell<=0 ||
794 zcell>=nCellZ || xcell>=nCellX) return;
795 Int_t izcell = (Int_t) zcell;
796 Int_t ixcell = (Int_t) xcell;
797 Float_t zc = zcell - izcell - 0.5;
798 Float_t xc = xcell - ixcell - 0.5;
799 for (Int_t iz=1; iz<=ngamz; iz++) {
800 Int_t kzg = izcell + iz - nz3;
801 if (kzg<=0 || kzg>nCellZ) continue;
802 Float_t zg = (Float_t)(iz-nz3) - zc;
803 for (Int_t ix=1; ix<=ngamx; ix++) {
804 Int_t kxg = ixcell + ix - nx3;
805 if (kxg<=0 || kxg>nCellX) continue;
806 Float_t xg = (Float_t)(ix-nx3) - xc;
808 // Now calculate pad response
809 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
810 qpad += qNoise*rnor2;
811 if (qpad<0) continue;
813 // Fill the array with pad response ID and amplitude
814 new(ldigits[cpvDigits->GetEntriesFast()]) CPVDigit(kxg,kzg,qpad);
815 // printf("(%2d,%2d,%5.3f) ",kxg,kzg,qpad);
823 //____________________________________________________________________________
824 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
825 // ------------------------------------------------------------------------
826 // Calculate the amplitude in one CPV pad using the
827 // cumulative pad response function
828 // Author: Yuri Kharlov (after Serguei Sadovski)
830 // ------------------------------------------------------------------------
832 Double_t dz = fGeom->GetPadSizeZ() / 2;
833 Double_t dx = fGeom->GetPadSizePhi() / 2;
834 Double_t z = zhit * fGeom->GetPadSizeZ();
835 Double_t x = xhit * fGeom->GetPadSizePhi();
836 Double_t amplitude = qhit *
837 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
838 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
839 return (Float_t)amplitude;
842 //____________________________________________________________________________
843 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
844 // ------------------------------------------------------------------------
845 // Cumulative pad response function
846 // It includes several terms from the CF decomposition in electrostatics
847 // Note: this cumulative function is wrong since omits some terms
848 // but the cell amplitude obtained with it is correct because
849 // these omitting terms cancel
850 // Author: Yuri Kharlov (after Serguei Sadovski)
852 // ------------------------------------------------------------------------
854 const Double_t a=1.0;
855 const Double_t b=0.7;
857 Double_t r2 = x*x + y*y;
859 Double_t cumulPRF = 0;
860 for (Int_t i=0; i<=4; i++) {
861 Double_t b1 = (2*i + 1) * b;
862 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
864 cumulPRF *= a/(2*TMath::Pi());