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
24 // Layout EMC + CPV has name IHEP:
25 // Produces hits for CPV, cumulated hits
27 // Layout EMC + CPV + PPSD has name GPS:
28 // Produces hits for CPV, cumulated hits
30 //*-- Author: Yves Schutz (SUBATECH)
33 // --- ROOT system ---
41 // --- Standard library ---
46 #include <strstream.h>
48 // --- AliRoot header files ---
50 #include "AliPHOSv1.h"
51 #include "AliPHOSHit.h"
52 #include "AliPHOSDigit.h"
53 #include "AliPHOSReconstructioner.h"
60 //____________________________________________________________________________
61 AliPHOSv1::AliPHOSv1()
65 // Create an empty array of AliPHOSCPVModule to satisfy
66 // AliPHOSv1::Streamer when reading root file
68 fReconstructioner = 0;
69 fTrackSegmentMaker = 0;
71 fHits = new TClonesArray("AliPHOSHit",1000) ;
73 // if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
74 // Error("AliPHOSv1","Can not create array of EMC modules");
78 // if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
79 // Error("AliPHOSv1","Can not create array of CPV modules");
85 //____________________________________________________________________________
86 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
89 // ctor : title is used to identify the layout
90 // GPS2 = 5 modules (EMC + PPSD)
91 // IHEP = 5 modules (EMC + CPV )
92 // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
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).
100 fPinElectronicNoise = 0.010 ;
101 fDigitThreshold = 0.01 ; // 1 GeV
103 fDigitizeB = 10000000. ;
106 // We do not want to save in TreeH the raw hits
107 // But save the cumulated hits instead (need to create the branch myself)
108 // It is put in the Digit Tree because the TreeH is filled after each primary
109 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
111 fHits= new TClonesArray("AliPHOSHit",1000) ;
115 fReconstructioner = 0;
116 fTrackSegmentMaker = 0;
118 fIshunt = 1 ; // All hits are associated with primary particles
120 // Create array of EMC modules of the size of PHOS modules number
122 // if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
123 // Error("AliPHOSv1","Can not create array of EMC modules");
126 // TClonesArray &lemcmodule = *fEMCModules;
127 // for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
129 // Create array of CPV modules for the IHEP's version of CPV
131 // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
132 // // Create array of CPV modules of the size of PHOS modules number
134 // if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
135 // Error("AliPHOSv1","Can not create array of CPV modules");
138 // TClonesArray &lcpvmodule = *fCPVModules;
139 // for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
142 // // Create an empty array of AliPHOSCPVModule to satisfy
143 // // AliPHOSv1::Streamer when writing root file
145 // fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
150 //____________________________________________________________________________
151 AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
152 AliPHOSv0(name,title)
154 // ctor : title is used to identify the layout
155 // GPS2 = 5 modules (EMC + PPSD)
156 // We use 2 arrays of hits :
158 // - fHits (the "normal" one), which retains the hits associated with
159 // the current primary particle being tracked
160 // (this array is reset after each primary has been tracked).
162 // - fTmpHits, which retains all the hits of the current event. It
163 // is used for the digitization part.
165 fPinElectronicNoise = 0.010 ;
167 // We do not want to save in TreeH the raw hits
168 //fHits = new TClonesArray("AliPHOSHit",100) ;
171 fHits= new TClonesArray("AliPHOSHit",1000) ;
175 fIshunt = 1 ; // All hits are associated with primary particles
177 // gets an instance of the geometry parameters class
178 fGeom = AliPHOSGeometry::GetInstance(title, "") ;
180 if (fGeom->IsInitialized() )
181 cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
183 cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
185 // Defining the PHOS Reconstructioner
187 fReconstructioner = Reconstructioner ;
191 //____________________________________________________________________________
192 AliPHOSv1::~AliPHOSv1()
214 if ( fEmcRecPoints ) {
215 fEmcRecPoints->Delete() ;
216 delete fEmcRecPoints ;
220 if ( fPpsdRecPoints ) {
221 fPpsdRecPoints->Delete() ;
222 delete fPpsdRecPoints ;
226 if ( fTrackSegments ) {
227 fTrackSegments->Delete() ;
228 delete fTrackSegments ;
234 //____________________________________________________________________________
235 void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits, Int_t trackpid, TLorentzVector p, Float_t * lpos)
237 // Add a hit to the hit list.
238 // A PHOS hit is the sum of all hits in a single crystal
239 // or in a single PPSD gas cell
244 Bool_t deja = kFALSE ;
246 newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid, p, lpos) ;
248 for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
249 curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
250 if( *curHit == *newHit ) {
251 *curHit = *curHit + *newHit ;
257 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
264 //____________________________________________________________________________
265 void AliPHOSv1::Hits2SDigits(){
266 //Collects all hits in the same active volume into digit
271 AliPHOSDigit * newdigit ;
272 AliPHOSDigit * curdigit ;
273 Bool_t deja = kFALSE ;
277 for (itrack=0; itrack<gAlice->GetNtrack(); itrack++){
279 //=========== Get the Hits Tree for the Primary track itrack
281 gAlice->TreeH()->GetEvent(itrack);
284 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
285 hit = (AliPHOSHit*)fHits->At(i) ;
287 // Assign primary number only if contribution is significant
288 if( hit->GetEnergy() > fDigitThreshold)
289 newdigit = new AliPHOSDigit( hit->GetPrimary(), hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
291 newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
294 for ( j = 0 ; j < fnSdigits ; j++) {
295 curdigit = (AliPHOSDigit*) fSDigits->At(j) ;
296 if ( *curdigit == *newdigit) {
297 *curdigit = *curdigit + *newdigit ;
303 new((*fSDigits)[fnSdigits]) AliPHOSDigit(* newdigit) ;
310 } // loop over tracks
314 fnSdigits = fSDigits->GetEntries() ;
315 fSDigits->Expand(fnSdigits) ;
317 for (i = 0 ; i < fnSdigits ; i++) {
318 AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ;
319 digit->SetIndexInList(i) ;
322 gAlice->TreeS()->Fill() ;
324 TTree *ts = gAlice->TreeS();
325 TDirectory *dir = ts->GetDirectory() ;
326 cout << "dir name is " << dir->GetName() << endl ;
327 // ts->Write(0,TObject::kOverwrite) ;
329 gAlice->TreeS()->GetBranch("PHOS")->Fill();
330 gAlice->TreeS()->GetBranch("PHOS")->Write();
334 //____________________________________________________________________________
335 void AliPHOSv1::SDigits2Digits(){
336 //Adds noise to the summable digits and removes everething below thresholds
337 //Note, that sDigits should be SORTED in accordance with abs ID.
340 gAlice->TreeS()->GetEvent(0) ;
342 cout << "fSdigits " << fSDigits << " " << fSDigits->GetEntries() << endl ;
345 // Noise induced by the PIN diode of the PbWO crystals
346 Int_t iCurSDigit = 0 ;
347 //we assume, that there is al least one EMC digit...
348 Int_t idCurSDigit = ((AliPHOSDigit *)fSDigits->At(0))->GetId() ;
350 cout << "fDigits " << fDigits << " " <<idCurSDigit<< endl ;
353 for(absID = 1; absID < fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); absID++){
355 cout << "absID " << absID << " " << idCurSDigit << endl ;
357 Float_t noise = gRandom->Gaus(0., fPinElectronicNoise) ;
358 if(absID < idCurSDigit ){
359 cout << "In < idC " << noise << " " << fDigitThreshold << endl ;
360 if(noise >fDigitThreshold ){
361 cout << "noise " << absID << " " << noise << endl;
362 new((*fDigits)[fNdigits]) AliPHOSDigit( -1,absID,Digitize(noise) ) ;
363 cout << "FN " << fNdigits << endl ;
367 else{ //add noise and may be remove the true hit
368 cout << "correcting digit " << iCurSDigit << endl ;
369 Float_t signal = noise + Calibrate(((AliPHOSDigit *)fSDigits->At(iCurSDigit))->GetAmp()) ;
370 cout << "signal " << signal << endl ;
371 if( signal >fDigitThreshold ){
372 cout << "signal " << signal << endl ;
373 AliPHOSDigit * digit = (AliPHOSDigit*) fSDigits->At(iCurSDigit) ;
374 new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
375 ((AliPHOSDigit *)fDigits->At(fNdigits))->SetAmp(Digitize(signal));
376 cout << "fNdigits " << fNdigits << endl ;
380 if(iCurSDigit < fSDigits->GetEntries()-1){
382 idCurSDigit = ((AliPHOSDigit*)fSDigits->At(iCurSDigit))->GetId() ;
385 idCurSDigit = 10000000; //no real hits left
390 //remove PPSD/CPV digits below thresholds
392 for(idigit = iCurSDigit; idigit < fSDigits->GetEntries() ; idigit++){ //loop over CPV/PPSD digits
394 AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(idigit) ;
395 Float_t ene = Calibrate(digit->GetAmp()) ;
398 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
399 if ( relid[0] > fGeom->GetNCPVModules() ){ //ppsd
400 if ( ( (relid[1] > 0) && (ene > fPpsdEnergyThreshold)) || //PPSD digit
401 ( (relid[1] < 0) && (ene > fCpvEnergyThreshold ) ) ) //CPV digit
402 new((*fDigits)[fNdigits]) AliPHOSDigit( *digit ) ;
407 fDigits->Compress() ;
409 fNdigits = fDigits->GetEntries() ;
410 fDigits->Expand(fNdigits) ;
413 for (i = 0 ; i < fNdigits ; i++) {
414 AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ;
415 digit->SetIndexInList(i) ;
419 gAlice->TreeD()->Fill() ;
421 gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
425 //___________________________________________________________________________
426 void AliPHOSv1::MakeBranch(Option_t* opt, char *file)
431 // Create new branche in the current Root Tree in the digit Tree
432 AliDetector::MakeBranch(opt) ;
435 cH = strstr(opt,"S");
436 //Create a branch for SDigits
439 sprintf(branchname,"%s",GetName());
443 fSDigits = new TClonesArray("AliPHOSDigit",1000);
445 cout << " AliPHOSv1::MakeBranch : " << file << endl ;
446 gAlice->MakeBranchInTree(gAlice->TreeS(),branchname,&fSDigits,fBufferSize,file);
449 cH = strstr(opt,"D");
450 //Create a branch for Digits
453 sprintf(branchname,"%s",GetName());
457 fDigits = new TClonesArray("AliPHOSDigit",1000);
460 gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fSDigits,fBufferSize,file);
463 cH = strstr(opt,"R");
464 //Create a branch for Reconstruction
468 Int_t splitlevel = 0 ;
470 fEmcRecPoints->Delete() ;
472 if ( fEmcRecPoints && gAlice->TreeR() ) {
473 sprintf(branchname,"%sEmcRP",GetName()) ;
474 gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
477 fPpsdRecPoints->Delete() ;
479 if ( fPpsdRecPoints && gAlice->TreeR() ) {
480 sprintf(branchname,"%sPpsdRP",GetName()) ;
481 gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
484 fTrackSegments->Delete() ;
486 if ( fTrackSegments && gAlice->TreeR() ) {
487 sprintf(branchname,"%sTS",GetName()) ;
488 gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
491 fRecParticles->Delete() ;
493 if ( fRecParticles && gAlice->TreeR() ) {
494 sprintf(branchname,"%sRP",GetName()) ;
495 gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
504 //_____________________________________________________________________________
505 void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
507 // 1. Reinitializes the existing RecPoint, TrackSegment, and RecParticles Lists and
508 // 2. Creates TreeR with a branch for each list
509 // 3. Steers the reconstruction processes
510 // 4. Saves the 3 lists in TreeR
511 // 5. Write the Tree to File
513 fReconstructioner = Reconstructioner ;
517 // gAlice->MakeTree("R") ;
523 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
525 printf("Reconstruction: %d %d %d %d\n",
526 fEmcRecPoints->GetEntries(),fPpsdRecPoints->GetEntries(),
527 fTrackSegments->GetEntries(),fRecParticles->GetEntries());
529 // 4. Expand or Shrink the arrays to the proper size
533 size = fEmcRecPoints->GetEntries() ;
534 fEmcRecPoints->Expand(size) ;
536 size = fPpsdRecPoints->GetEntries() ;
537 fPpsdRecPoints->Expand(size) ;
539 size = fTrackSegments->GetEntries() ;
540 fTrackSegments->Expand(size) ;
542 size = fRecParticles->GetEntries() ;
543 fRecParticles->Expand(size) ;
545 gAlice->TreeR()->Fill() ;
548 gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
550 // Deleting reconstructed objects
551 ResetReconstruction();
555 //____________________________________________________________________________
556 void AliPHOSv1::ResetHits()
558 // Reset hit tree for EMC and CPV
559 // Yuri Kharlov, 28 September 2000
561 AliDetector::ResetHits();
562 // for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fEMCModules)[i]) -> Clear();
563 // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
564 // for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
568 //____________________________________________________________________________
569 void AliPHOSv1::ResetReconstruction()
571 // Deleting reconstructed objects
573 if ( fEmcRecPoints ) fEmcRecPoints ->Delete();
574 if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
575 if ( fTrackSegments ) fTrackSegments->Delete();
576 if ( fRecParticles ) fRecParticles ->Delete();
580 //____________________________________________________________________________
581 //void AliPHOSv1::SDigits2Digits() {
582 // // Adds the noise to the summable digits and keeps digits above a threshold
583 // // To make a digit
586 //____________________________________________________________________________
587 void AliPHOSv1::SetTreeAddress()
590 AliPHOS::SetTreeAddress();
592 // //Branch address for TreeR: RecPpsdRecPoint
593 // TTree *treeR = gAlice->TreeR();
594 // if ( treeR && fPpsdRecPoints ) {
595 // branch = treeR->GetBranch("PHOSPpsdRP");
596 // if (branch) branch->SetAddress(&fPpsdRecPoints) ;
599 // Set branch address for the Hits Tree for hits in EMC modules
600 // Yuri Kharlov, 23 November 2000.
602 // for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
604 // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
605 // Yuri Kharlov, 28 September 2000.
607 // if ( strcmp(fGeom->GetName(),"IHEP") == 0 || strcmp(fGeom->GetName(),"MIXT") == 0 ) {
608 // for( Int_t i=0; i<fGeom->GetNCPVModules(); i++ ) GetCPVModule(i).SetTreeAddress("CPV",i+1);
613 //____________________________________________________________________________
615 void AliPHOSv1::StepManager(void)
617 // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
619 // if (gMC->IsTrackEntering())
620 // cout << "Track enters the volume " << gMC->CurrentVolName() << endl;
621 // if (gMC->IsTrackExiting())
622 // cout << "Track leaves the volume " << gMC->CurrentVolName() << endl;
624 Int_t relid[4] ; // (box, layer, row, column) indices
625 Int_t absid ; // absolute cell ID number
626 Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited
627 TLorentzVector pos ; // Lorentz vector of the track current position
628 TLorentzVector pmom ; //momentum of the particle initiated hit
629 Float_t xyd[2] ; //local posiiton of the entering
630 Bool_t entered = kFALSE ;
633 Int_t tracknumber = gAlice->CurrentTrack() ;
634 Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
635 TString name = fGeom->GetName() ;
638 if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle,
639 // but may be without energy deposition
641 // Current position of the hit in the local ref. system
642 gMC -> TrackPosition(pos);
643 Float_t xyzm[3], xyzd[3] ;
645 for (i=0; i<3; i++) xyzm[i] = pos[i];
646 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
650 // Current momentum of the hit's track in the local ref. system
651 gMC -> TrackMomentum(pmom);
652 Float_t pm[3], pd[3];
653 for (i=0; i<3; i++) pm[i] = pmom[i];
654 gMC -> Gmtod (pm, pd, 2); // transform 3-momentum from master to daughter system
659 trackpid = gMC->TrackPid();
660 entered = kTRUE ; // Mark to create hit even withou energy deposition
665 if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD
667 if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
669 gMC->TrackPosition(pos) ;
673 xyze[3] = gMC->Edep() ;
675 if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering PPSD
676 gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number
677 if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
678 relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
680 gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
681 // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
682 // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
683 gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
684 gMC->CurrentVolID(relid[3]) ; // get the column number
686 // get the absolute Id number
688 fGeom->RelToAbsNumbering(relid, absid) ;
690 // add current hit to the hit list
691 AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
694 } // there is deposited energy
695 } // We are inside the gas of the CPV
696 } // GPS2 configuration
698 if ( name == "IHEP" || name == "MIXT" ) { // ======> CPV is a IHEP's one
700 // Yuri Kharlov, 28 September 2000
702 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
704 gMC->TrackCharge() != 0) {
706 // Digitize the current CPV hit:
708 // 1. find pad response and
711 gMC->CurrentVolOffID(3,moduleNumber);
715 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
716 CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
721 Int_t idigit,ndigits;
723 // 2. go through the current digit list and sum digits in pads
725 ndigits = cpvDigits->GetEntriesFast();
726 for (idigit=0; idigit<ndigits-1; idigit++) {
727 AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
728 Float_t x1 = cpvDigit1->GetXpad() ;
729 Float_t z1 = cpvDigit1->GetYpad() ;
730 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
731 AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
732 Float_t x2 = cpvDigit2->GetXpad() ;
733 Float_t z2 = cpvDigit2->GetYpad() ;
734 if (x1==x2 && z1==z2) {
735 Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
736 cpvDigit2->SetQpad(qsum) ;
737 cpvDigits->RemoveAt(idigit) ;
741 cpvDigits->Compress() ;
743 // 3. add digits to temporary hit list fTmpHits
745 ndigits = cpvDigits->GetEntriesFast();
746 for (idigit=0; idigit<ndigits; idigit++) {
747 AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
748 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
749 relid[1] =-1 ; // means CPV
750 relid[2] = cpvDigit->GetXpad() ; // column number of a pad
751 relid[3] = cpvDigit->GetYpad() ; // row number of a pad
753 // get the absolute Id number
754 fGeom->RelToAbsNumbering(relid, absid) ;
756 // add current digit to the temporary hit list
760 xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
761 primary = -1; // No need in primary for CPV
762 AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
764 if (cpvDigit->GetQpad() > 0.02) {
765 xmean += cpvDigit->GetQpad() * (cpvDigit->GetXpad() + 0.5);
766 zmean += cpvDigit->GetQpad() * (cpvDigit->GetYpad() + 0.5);
767 qsum += cpvDigit->GetQpad();
772 } // end of IHEP configuration
775 if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
776 gMC->TrackPosition(pos) ;
780 xyze[3] = gMC->Edep() ;
783 if ( (xyze[3] != 0) || entered ) { // Track is inside the crystal and deposits some energy or just entered
785 gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
787 if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
788 relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
790 relid[1] = 0 ; // means PBW04
791 gMC->CurrentVolOffID(4, relid[2]) ; // get the row number inside the module
792 gMC->CurrentVolOffID(3, relid[3]) ; // get the cell number inside the module
794 // get the absolute Id number
795 fGeom->RelToAbsNumbering(relid, absid) ;
797 // add current hit to the hit list
798 AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
801 } // there is deposited energy
802 } // we are inside a PHOS Xtal
807 //____________________________________________________________________________
808 void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
810 // ------------------------------------------------------------------------
811 // Digitize one CPV hit:
812 // On input take exact 4-momentum p and position zxhit of the hit,
813 // find the pad response around this hit and
814 // put the amplitudes in the pads into array digits
816 // Author: Yuri Kharlov (after Serguei Sadovsky)
818 // ------------------------------------------------------------------------
820 const Float_t kCelWr = fGeom->GetPadSizePhi()/2; // Distance between wires (2 wires above 1 pad)
821 const Float_t kDetR = 0.1; // Relative energy fluctuation in track for 100 e-
822 const Float_t kdEdx = 4.0; // Average energy loss in CPV;
823 const Int_t kNgamz = 5; // Ionization size in Z
824 const Int_t kNgamx = 9; // Ionization size in Phi
825 const Float_t kNoise = 0.03; // charge noise in one pad
829 // Just a reminder on axes notation in the CPV module:
830 // axis Z goes along the beam
831 // axis X goes across the beam in the module plane
832 // axis Y is a normal to the module plane showing from the IP
834 Float_t hitX = zxhit[0];
835 Float_t hitZ =-zxhit[1];
838 Float_t pNorm = p.Py();
839 Float_t eloss = kdEdx;
841 // cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
843 Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
844 Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
845 gRandom->Rannor(rnor1,rnor2);
846 eloss *= (1 + kDetR*rnor1) *
847 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
848 Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
849 Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
850 Float_t zhit2 = zhit1 + dZY;
851 Float_t xhit2 = xhit1 + dXY;
853 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
854 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
858 if (iwht1==iwht2) { // incline 1-wire hit
860 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
861 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
863 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
864 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
867 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
869 Int_t iwht3 = (iwht1 + iwht2) / 2;
870 Float_t xwht1 = (iwht1 + 0.5) * kCelWr; // wire 1
871 Float_t xwht2 = (iwht2 + 0.5) * kCelWr; // wire 2
872 Float_t xwht3 = (iwht3 + 0.5) * kCelWr; // wire 3
873 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
874 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
875 Float_t dxw1 = xhit1 - xwr13;
876 Float_t dxw2 = xhit2 - xwr23;
877 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
878 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
879 Float_t egm3 = kCelWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + kCelWr );
880 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
882 zxe[2][0] = eloss * egm1;
883 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
885 zxe[2][1] = eloss * egm2;
886 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
888 zxe[2][2] = eloss * egm3;
890 else { // incline 2-wire hit
892 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
893 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
894 Float_t xwr12 = (xwht1 + xwht2) / 2;
895 Float_t dxw1 = xhit1 - xwr12;
896 Float_t dxw2 = xhit2 - xwr12;
897 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
898 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
899 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
901 zxe[2][0] = eloss * egm1;
902 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
904 zxe[2][1] = eloss * egm2;
907 // Finite size of ionization region
909 Int_t nCellZ = fGeom->GetNumberOfCPVPadsZ();
910 Int_t nCellX = fGeom->GetNumberOfCPVPadsPhi();
911 Int_t nz3 = (kNgamz+1)/2;
912 Int_t nx3 = (kNgamx+1)/2;
913 cpvDigits->Expand(nIter*kNgamx*kNgamz);
914 TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
916 for (Int_t iter=0; iter<nIter; iter++) {
918 Float_t zhit = zxe[0][iter];
919 Float_t xhit = zxe[1][iter];
920 Float_t qhit = zxe[2][iter];
921 Float_t zcell = zhit / fGeom->GetPadSizeZ();
922 Float_t xcell = xhit / fGeom->GetPadSizePhi();
923 if ( zcell<=0 || xcell<=0 ||
924 zcell>=nCellZ || xcell>=nCellX) return;
925 Int_t izcell = (Int_t) zcell;
926 Int_t ixcell = (Int_t) xcell;
927 Float_t zc = zcell - izcell - 0.5;
928 Float_t xc = xcell - ixcell - 0.5;
929 for (Int_t iz=1; iz<=kNgamz; iz++) {
930 Int_t kzg = izcell + iz - nz3;
931 if (kzg<=0 || kzg>nCellZ) continue;
932 Float_t zg = (Float_t)(iz-nz3) - zc;
933 for (Int_t ix=1; ix<=kNgamx; ix++) {
934 Int_t kxg = ixcell + ix - nx3;
935 if (kxg<=0 || kxg>nCellX) continue;
936 Float_t xg = (Float_t)(ix-nx3) - xc;
938 // Now calculate pad response
939 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
940 qpad += kNoise*rnor2;
941 if (qpad<0) continue;
943 // Fill the array with pad response ID and amplitude
944 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
950 //____________________________________________________________________________
951 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
952 // ------------------------------------------------------------------------
953 // Calculate the amplitude in one CPV pad using the
954 // cumulative pad response function
955 // Author: Yuri Kharlov (after Serguei Sadovski)
957 // ------------------------------------------------------------------------
959 Double_t dz = fGeom->GetPadSizeZ() / 2;
960 Double_t dx = fGeom->GetPadSizePhi() / 2;
961 Double_t z = zhit * fGeom->GetPadSizeZ();
962 Double_t x = xhit * fGeom->GetPadSizePhi();
963 Double_t amplitude = qhit *
964 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
965 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
966 return (Float_t)amplitude;
969 //____________________________________________________________________________
970 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
971 // ------------------------------------------------------------------------
972 // Cumulative pad response function
973 // It includes several terms from the CF decomposition in electrostatics
974 // Note: this cumulative function is wrong since omits some terms
975 // but the cell amplitude obtained with it is correct because
976 // these omitting terms cancel
977 // Author: Yuri Kharlov (after Serguei Sadovski)
979 // ------------------------------------------------------------------------
981 const Double_t kA=1.0;
982 const Double_t kB=0.7;
984 Double_t r2 = x*x + y*y;
986 Double_t cumulPRF = 0;
987 for (Int_t i=0; i<=4; i++) {
988 Double_t b1 = (2*i + 1) * kB;
989 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
991 cumulPRF *= kA/(2*TMath::Pi());