Chamber half-planes of stations 3-5 at different z-positions.
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
CommitLineData
7587f5a5 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
b2a60966 16/* $Id$ */
5f20d3fb 17
7587f5a5 18//_________________________________________________________________________
5f20d3fb 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)
b2a60966 23
7587f5a5 24
25// --- ROOT system ---
bea63bea 26
27#include "TBRIK.h"
28#include "TNode.h"
7587f5a5 29#include "TRandom.h"
94de3818 30#include "TTree.h"
7587f5a5 31
5f20d3fb 32
7587f5a5 33// --- Standard library ---
34
de9ec31b 35#include <stdio.h>
36#include <string.h>
37#include <stdlib.h>
38#include <strstream.h>
7587f5a5 39
40// --- AliRoot header files ---
41
42#include "AliPHOSv1.h"
43#include "AliPHOSHit.h"
44#include "AliPHOSDigit.h"
bea63bea 45#include "AliPHOSReconstructioner.h"
7587f5a5 46#include "AliRun.h"
47#include "AliConst.h"
94de3818 48#include "AliMC.h"
7587f5a5 49
50ClassImp(AliPHOSv1)
51
bea63bea 52//____________________________________________________________________________
53AliPHOSv1::AliPHOSv1()
54{
5f20d3fb 55 // ctor
bea63bea 56 fNTmpHits = 0 ;
fa412d9b 57 fTmpHits = 0 ;
fc879520 58
3d402178 59 // Create an empty array of AliPHOSCPVModule to satisfy
fa412d9b 60 // AliPHOSv1::Streamer when reading root file
61
3d402178 62 if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
fa412d9b 63 Error("AliPHOSv1","Can not create array of CPV modules");
64 exit(1);
65 }
fc879520 66
bea63bea 67}
68
7587f5a5 69//____________________________________________________________________________
70AliPHOSv1::AliPHOSv1(const char *name, const char *title):
e04976bd 71AliPHOSv0(name,title)
7587f5a5 72{
5f20d3fb 73 // ctor : title is used to identify the layout
fa412d9b 74 // GPS2 = 5 modules (EMC + PPSD)
75 // IHEP = 5 modules (EMC + CPV )
5f20d3fb 76 // We use 2 arrays of hits :
77 //
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).
81 //
82 // - fTmpHits, which retains all the hits of the current event. It
83 // is used for the digitization part.
fa412d9b 84
5f20d3fb 85 fPinElectronicNoise = 0.010 ;
b95514c0 86 fDigitThreshold = 0.1 ; // 1 GeV
5f20d3fb 87
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
fa412d9b 91 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
5f20d3fb 92
93 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
94
95 fNTmpHits = fNhits = 0 ;
96
97 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
98
5f20d3fb 99 fIshunt = 1 ; // All hits are associated with primary particles
fa412d9b 100
101 // Create array of CPV modules for the IHEP's version of CPV
102
103 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
3d402178 104 // Create array of AliPHOSCPVmodule of the size of PHOS modules number
fa412d9b 105
3d402178 106 if ( NULL==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
fa412d9b 107 Error("AliPHOSv1","Can not create array of CPV modules");
108 exit(1);
109 }
110 TClonesArray &lcpvmodule = *fCPVModules;
3d402178 111 for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
fa412d9b 112 }
113 else {
3d402178 114 // Create an empty array of AliPHOSCPVModule to satisfy
fa412d9b 115 // AliPHOSv1::Streamer when writing root file
116
3d402178 117 fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
fa412d9b 118
119 }
5f20d3fb 120}
121
122//____________________________________________________________________________
123AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
124 AliPHOSv0(name,title)
125{
126 // ctor : title is used to identify the layout
127 // GPS2 = 5 modules (EMC + PPSD)
128 // We use 2 arrays of hits :
129 //
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).
133 //
134 // - fTmpHits, which retains all the hits of the current event. It
135 // is used for the digitization part.
136
137 fPinElectronicNoise = 0.010 ;
138
139 // We do not want to save in TreeH the raw hits
140 //fHits = new TClonesArray("AliPHOSHit",100) ;
141
142 fDigits = new TClonesArray("AliPHOSDigit",1000) ;
143 fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
144
145 fNTmpHits = fNhits = 0 ;
146
147 fIshunt = 1 ; // All hits are associated with primary particles
148
149 // gets an instance of the geometry parameters class
150 fGeom = AliPHOSGeometry::GetInstance(title, "") ;
151
152 if (fGeom->IsInitialized() )
88bdfa12 153 cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
5f20d3fb 154 else
fc879520 155 cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
5f20d3fb 156
157 // Defining the PHOS Reconstructioner
158
159 fReconstructioner = Reconstructioner ;
160
7587f5a5 161}
b2a60966 162
7587f5a5 163//____________________________________________________________________________
bea63bea 164AliPHOSv1::~AliPHOSv1()
b2a60966 165{
bea63bea 166 // dtor
5f20d3fb 167
8dfa469d 168 if ( fTmpHits) {
169 fTmpHits->Delete() ;
170 delete fTmpHits ;
171 fTmpHits = 0 ;
172 }
5f20d3fb 173
fc879520 174// if ( fEmcRecPoints ) {
175// fEmcRecPoints->Delete() ;
176// delete fEmcRecPoints ;
177// fEmcRecPoints = 0 ;
178// }
179
180// if ( fPpsdRecPoints ) {
181// fPpsdRecPoints->Delete() ;
182// delete fPpsdRecPoints ;
183// fPpsdRecPoints = 0 ;
184// }
5f20d3fb 185
fc879520 186// if ( fTrackSegments ) {
187// fTrackSegments->Delete() ;
188// delete fTrackSegments ;
189// fTrackSegments = 0 ;
190// }
5f20d3fb 191
7587f5a5 192}
193
7587f5a5 194//____________________________________________________________________________
bea63bea 195void AliPHOSv1::AddHit(Int_t shunt, Int_t primary, Int_t tracknumber, Int_t Id, Float_t * hits)
196{
197 // Add a hit to the hit list.
5f20d3fb 198 // A PHOS hit is the sum of all hits in a single crystal
199 // or in a single PPSD gas cell
bea63bea 200
5f20d3fb 201 Int_t hitCounter ;
202 TClonesArray &ltmphits = *fTmpHits ;
bea63bea 203 AliPHOSHit *newHit ;
5f20d3fb 204 AliPHOSHit *curHit ;
205 Bool_t deja = kFALSE ;
bea63bea 206
5f20d3fb 207 // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
bea63bea 208
209 newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
210
5f20d3fb 211 // We do not want to save in TreeH the raw hits
bea63bea 212 // TClonesArray &lhits = *fHits;
bea63bea 213
5f20d3fb 214 for ( hitCounter = 0 ; hitCounter < fNTmpHits && !deja ; hitCounter++ ) {
215 curHit = (AliPHOSHit*) ltmphits[hitCounter] ;
216 if( *curHit == *newHit ) {
217 *curHit = *curHit + *newHit ;
218 deja = kTRUE ;
219 }
220 }
221
222 if ( !deja ) {
223 new(ltmphits[fNTmpHits]) AliPHOSHit(*newHit) ;
224 fNTmpHits++ ;
225 }
226
bea63bea 227 // We do not want to save in TreeH the raw hits
228 // new(lhits[fNhits]) AliPHOSHit(*newHit) ;
229 // fNhits++ ;
230
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).
234
235 delete newHit;
236
237}
238
5f20d3fb 239//___________________________________________________________________________
240Int_t AliPHOSv1::Digitize(Float_t Energy)
241{
242 // Applies the energy calibration
243
244 Float_t fB = 100000000. ;
245 Float_t fA = 0. ;
246 Int_t chan = Int_t(fA + Energy*fB ) ;
247 return chan ;
248}
bea63bea 249
250//___________________________________________________________________________
251void AliPHOSv1::FinishEvent()
252{
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
256
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.
260
5f20d3fb 261
bea63bea 262 Int_t i ;
263 Int_t relid[4];
264 Int_t j ;
265 TClonesArray &lDigits = *fDigits ;
fa412d9b 266 AliPHOSHit * hit ;
bea63bea 267 AliPHOSDigit * newdigit ;
268 AliPHOSDigit * curdigit ;
269 Bool_t deja = kFALSE ;
270
5f20d3fb 271 for ( i = 0 ; i < fNTmpHits ; i++ ) {
272 hit = (AliPHOSHit*)fTmpHits->At(i) ;
4a2ca5e9 273
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() ) ) ;
277 else
fa412d9b 278 newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
bea63bea 279 deja =kFALSE ;
280 for ( j = 0 ; j < fNdigits ; j++) {
281 curdigit = (AliPHOSDigit*) lDigits[j] ;
282 if ( *curdigit == *newdigit) {
283 *curdigit = *curdigit + *newdigit ;
284 deja = kTRUE ;
285 }
286 }
287 if ( !deja ) {
288 new(lDigits[fNdigits]) AliPHOSDigit(* newdigit) ;
289 fNdigits++ ;
290 }
291
292 delete newdigit ;
293 }
294
295 // Noise induced by the PIN diode of the PbWO crystals
296
297 Float_t energyandnoise ;
298 for ( i = 0 ; i < fNdigits ; i++ ) {
299 newdigit = (AliPHOSDigit * ) fDigits->At(i) ;
fa412d9b 300
bea63bea 301 fGeom->AbsToRelNumbering(newdigit->GetId(), relid) ;
302
303 if (relid[1]==0){ // Digits belong to EMC (PbW0_4 crystals)
304 energyandnoise = newdigit->GetAmp() + Digitize(gRandom->Gaus(0., fPinElectronicNoise)) ;
305
306 if (energyandnoise < 0 )
307 energyandnoise = 0 ;
308
309 if ( newdigit->GetAmp() < fDigitThreshold ) // if threshold not surpassed, remove digit from list
310 fDigits->RemoveAt(i) ;
311 }
312 }
313
314 fDigits->Compress() ;
315
316 fNdigits = fDigits->GetEntries() ;
317 for (i = 0 ; i < fNdigits ; i++) {
318 newdigit = (AliPHOSDigit *) fDigits->At(i) ;
319 newdigit->SetIndexInList(i) ;
fa412d9b 320
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());
bea63bea 324 }
fa412d9b 325
5f20d3fb 326}
327
328//___________________________________________________________________________
329void AliPHOSv1::MakeBranch(Option_t* opt)
330{
331 // Create new branche in the current Root Tree in the digit Tree
5f20d3fb 332 AliDetector::MakeBranch(opt) ;
333
334 char branchname[10];
335 sprintf(branchname,"%s",GetName());
336 char *cdD = strstr(opt,"D");
337 if (fDigits && gAlice->TreeD() && cdD) {
338 gAlice->TreeD()->Branch(branchname, &fDigits, fBufferSize);
339 }
340
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) ;
347 }
348 }
bea63bea 349
fa412d9b 350 // Create new branches CPV<i> for hits in CPV modules for IHEP geometry
351 // Yuri Kharlov, 28 September 2000.
352
353 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
354 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).MakeBranch(i+1);
355 }
356
bea63bea 357}
358
5f20d3fb 359//_____________________________________________________________________________
360void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
361{
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
367
368 fReconstructioner = Reconstructioner ;
369
370 char branchname[10] ;
371
372 // 1.
373
374 // gAlice->MakeTree("R") ;
375 Int_t splitlevel = 0 ;
376
fc879520 377 fEmcRecPoints->Delete() ;
5f20d3fb 378
379 if ( fEmcRecPoints && gAlice->TreeR() ) {
380 sprintf(branchname,"%sEmcRP",GetName()) ;
381
5f20d3fb 382 gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
383 }
384
fa412d9b 385 fPpsdRecPoints->Delete() ;
5f20d3fb 386
5f20d3fb 387
388 if ( fPpsdRecPoints && gAlice->TreeR() ) {
389 sprintf(branchname,"%sPpsdRP",GetName()) ;
390
5f20d3fb 391 gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
392 }
393
fc879520 394 fTrackSegments->Delete() ;
5f20d3fb 395
5f20d3fb 396 if ( fTrackSegments && gAlice->TreeR() ) {
397 sprintf(branchname,"%sTS",GetName()) ;
398 gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
399 }
400
fa412d9b 401 fRecParticles->Delete() ;
fc879520 402
fa412d9b 403 if (strcmp(fGeom->GetName(),"GPS2") == 0) {
404 if ( fRecParticles && gAlice->TreeR() ) {
405 sprintf(branchname,"%sRP",GetName()) ;
406 gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
407 }
5f20d3fb 408 }
409
410 // 3.
fa412d9b 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);
5f20d3fb 415
416 // 4. Expand or Shrink the arrays to the proper size
417
418 Int_t size ;
419
420 size = fEmcRecPoints->GetEntries() ;
421 fEmcRecPoints->Expand(size) ;
422
423 size = fPpsdRecPoints->GetEntries() ;
424 fPpsdRecPoints->Expand(size) ;
425
426 size = fTrackSegments->GetEntries() ;
427 fTrackSegments->Expand(size) ;
428
429 size = fRecParticles->GetEntries() ;
430 fRecParticles->Expand(size) ;
431
432 gAlice->TreeR()->Fill() ;
5f20d3fb 433 // 5.
434
4a2ca5e9 435 gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
5f20d3fb 436
437 // Deleting reconstructed objects
438 ResetReconstruction();
439
440}
441
fa412d9b 442//____________________________________________________________________________
443void AliPHOSv1::ResetHits()
444{
445 // Reset hit tree for CPV in IHEP geometry
446 // Yuri Kharlov, 28 September 2000
447
448 AliDetector::ResetHits();
449 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
3d402178 450 for (Int_t i=0; i<fGeom->GetNModules(); i++) ((AliPHOSCPVModule*)(*fCPVModules)[i]) -> Clear();
fa412d9b 451 }
452}
5f20d3fb 453//____________________________________________________________________________
454void AliPHOSv1::ResetDigits()
455{
456 // May sound strange, but cumulative hits are store in digits Tree
457 AliDetector::ResetDigits();
458 if( fTmpHits ) {
459 fTmpHits->Delete();
460 fNTmpHits = 0 ;
461 }
462}
463//____________________________________________________________________________
464void AliPHOSv1::ResetReconstruction()
465{
466 // Deleting reconstructed objects
467
468 if ( fEmcRecPoints ) fEmcRecPoints->Delete();
469 if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
470 if ( fTrackSegments ) fTrackSegments->Delete();
471 if ( fRecParticles ) fRecParticles->Delete();
472
473}
5f20d3fb 474
475//____________________________________________________________________________
476void AliPHOSv1::SetTreeAddress()
477{
478 // TBranch *branch;
479 AliPHOS::SetTreeAddress();
480
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) ;
fa412d9b 486// }
487
488 // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
489 // Yuri Kharlov, 28 September 2000.
490
491 if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
492 for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetCPVModule(i).SetTreeAddress(i+1);
493 }
494
5f20d3fb 495}
496
497//____________________________________________________________________________
498
7587f5a5 499void AliPHOSv1::StepManager(void)
500{
b2a60966 501 // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
b2a60966 502
7587f5a5 503 Int_t relid[4] ; // (box, layer, row, column) indices
fa412d9b 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
507 Int_t copy ;
7587f5a5 508
bea63bea 509 Int_t tracknumber = gAlice->CurrentTrack() ;
fa412d9b 510 Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
511 TString name = fGeom->GetName() ;
512
513 if ( name == "GPS2" ) { // ======> CPV is a GPS' PPSD
514
b2a60966 515 if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
7587f5a5 516 {
517 gMC->TrackPosition(pos) ;
518 xyze[0] = pos[0] ;
519 xyze[1] = pos[1] ;
520 xyze[2] = pos[2] ;
bea63bea 521 xyze[3] = gMC->Edep() ;
7587f5a5 522
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
530
531 // get the absolute Id number
532
bea63bea 533 fGeom->RelToAbsNumbering(relid, absid) ;
7587f5a5 534
bea63bea 535 // add current hit to the hit list
536 AddHit(fIshunt, primary, tracknumber, absid, xyze);
7587f5a5 537
538 } // there is deposited energy
fa412d9b 539 } // We are inside the gas of the CPV
540 } // GPS2 configuration
541
542 else if ( name == "IHEP" ) { // ======> CPV is a IHEP's one
543
544 // Yuri Kharlov, 28 September 2000
545
546 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
547 gMC->IsTrackEntering() &&
548 gMC->TrackCharge() != 0) {
549
550 // Charged track has just entered to the CPV sensitive plane
551
552 AliPHOSv1 &PHOS = *(AliPHOSv1*)gAlice->GetModule("PHOS");
553
554 Int_t ModuleNumber;
555 gMC->CurrentVolOffID(3,ModuleNumber);
556 ModuleNumber--;
557
558 // Current position of the hit in the CPV module ref. system
559
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
564 xyd[0] = xyzd[0];
565 xyd[1] =-xyzd[2];
566
567 // Current momentum of the hit's track in the CPV module ref. system
568
569 TLorentzVector pmom;
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
574 pmom[0] = pd[0];
575 pmom[1] =-pd[1];
576 pmom[2] =-pd[2];
577
578 // Current particle type of the hit's track
579
580 Int_t ipart = gMC->TrackPid();
581
582 // Add the current particle in the list of the CPV hits.
583
584 PHOS.GetCPVModule(ModuleNumber).AddHit(pmom,xyd,ipart);
585
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);
590
591 // Digitize the current CPV hit:
592
593 // 1. find pad response and
594
3d402178 595 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
fa412d9b 596 CPVDigitize(pmom,xyd,ModuleNumber,cpvDigits);
597
598 Float_t xmean = 0;
599 Float_t zmean = 0;
600 Float_t qsum = 0;
601 Int_t ndigits;
602
603 // 2. go through the current digit list and sum digits in pads
604
605 ndigits = cpvDigits->GetEntriesFast();
606 for (Int_t idigit=0; idigit<ndigits-1; idigit++) {
3d402178 607 AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
fa412d9b 608 Float_t x1 = cpvDigit1->GetXpad() ;
609 Float_t z1 = cpvDigit1->GetYpad() ;
610 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
3d402178 611 AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
fa412d9b 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) ;
618 }
619 }
620 }
621 cpvDigits->Compress() ;
622
623 // 3. add digits to temporary hit list fTmpHits
624
625 ndigits = cpvDigits->GetEntriesFast();
626 for (Int_t idigit=0; idigit<ndigits; idigit++) {
3d402178 627 AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
fa412d9b 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
632
633 // get the absolute Id number
634 fGeom->RelToAbsNumbering(relid, absid) ;
635
636 // add current digit to the temporary hit list
637 xyze[0] = 0. ;
638 xyze[1] = 0. ;
639 xyze[2] = 0. ;
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);
643
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();
648 }
649 }
650 delete cpvDigits;
651 }
652 } // end of IHEP configuration
7587f5a5 653
fa412d9b 654 if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
655 gMC->TrackPosition(pos) ;
656 xyze[0] = pos[0] ;
657 xyze[1] = pos[1] ;
658 xyze[2] = pos[2] ;
659 xyze[3] = gMC->Edep() ;
660
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
666
667 // get the absolute Id number
668
669 fGeom->RelToAbsNumbering(relid, absid) ;
670
671 // add current hit to the hit list
672
673 AddHit(fIshunt, primary,tracknumber, absid, xyze);
674
675 } // there is deposited energy
676 } // we are inside a PHOS Xtal
677}
678
679//____________________________________________________________________________
680void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
681{
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
687 //
688 // Author: Yuri Kharlov (after Serguei Sadovsky)
689 // 2 October 2000
690 // ------------------------------------------------------------------------
691
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
698
699 Float_t rnor1,rnor2;
700
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
705
3d402178 706// cout << __PRETTY_FUNCTION__ << ": YVK : Start digitization\n";
707
fa412d9b 708 Float_t hitX = zxhit[0];
709 Float_t hitZ =-zxhit[1];
710 Float_t pX = p.Px();
711 Float_t pZ =-p.Pz();
712 Float_t pNorm = p.Py();
713 Float_t E = dEdx;
714
3d402178 715// cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
716
fa412d9b 717 Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
718 Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
719 gRandom->Rannor(rnor1,rnor2);
720 E *= (1 + detR*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;
726
3d402178 727// cout << "CPVDigitize: YVK : "<<xhit1<<" "<<xhit2<<" | "<<zhit1<<" "<<zhit2<<" | "
728// << (Int_t)(xhit1/fGeom->GetPadSizePhi())<<" "<<(Int_t)(zhit1/fGeom->GetPadSizeZ())<<" "<<endl;
729
fa412d9b 730 Int_t iwht1 = (Int_t) (xhit1 / celWr); // wire (x) coordinate "in"
731 Int_t iwht2 = (Int_t) (xhit2 / celWr); // wire (x) coordinate "out"
732
733 Int_t nIter;
734 Float_t zxe[3][5];
735 if (iwht1==iwht2) { // incline 1-wire hit
736 nIter = 2;
737 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
738 zxe[1][0] = (iwht1 + 0.5) * celWr;
739 zxe[2][0] = E/2;
740 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
741 zxe[1][1] = (iwht1 + 0.5) * celWr;
742 zxe[2][1] = E/2;
743 }
744 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
745 nIter = 3;
746 Int_t iwht3 = (iwht1 + iwht2) / 2;
747 Float_t xwht1 = (iwht1 + 0.5) * celWr; // wire 1
748 Float_t xwht2 = (iwht2 + 0.5) * celWr; // wire 2
749 Float_t xwht3 = (iwht3 + 0.5) * celWr; // wire 3
750 Float_t xwr13 = (xwht1 + xwht3) / 2; // center 13
751 Float_t xwr23 = (xwht2 + xwht3) / 2; // center 23
752 Float_t dxw1 = xhit1 - xwr13;
753 Float_t dxw2 = xhit2 - xwr23;
754 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
755 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
756 Float_t egm3 = celWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
757 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
758 zxe[1][0] = xwht1;
759 zxe[2][0] = E * egm1;
760 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
761 zxe[1][1] = xwht2;
762 zxe[2][1] = E * egm2;
763 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
764 zxe[1][2] = xwht3;
765 zxe[2][2] = E * egm3;
766 }
767 else { // incline 2-wire hit
768 nIter = 2;
769 Float_t xwht1 = (iwht1 + 0.5) * celWr;
770 Float_t xwht2 = (iwht2 + 0.5) * celWr;
771 Float_t xwr12 = (xwht1 + xwht2) / 2;
772 Float_t dxw1 = xhit1 - xwr12;
773 Float_t dxw2 = xhit2 - xwr12;
774 Float_t egm1 = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
775 Float_t egm2 = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
776 zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
777 zxe[1][0] = xwht1;
778 zxe[2][0] = E * egm1;
779 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
780 zxe[1][1] = xwht2;
781 zxe[2][1] = E * egm2;
782 }
bea63bea 783
fa412d9b 784 // Finite size of ionization region
785
786 Int_t nCellZ = fGeom->GetNumberOfPadsZ();
787 Int_t nCellX = fGeom->GetNumberOfPadsPhi();
788 Int_t nz3 = (ngamz+1)/2;
789 Int_t nx3 = (ngamx+1)/2;
790 cpvDigits->Expand(nIter*ngamx*ngamz);
791 TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
792
793 for (Int_t iter=0; iter<nIter; iter++) {
3d402178 794// cout << "CPVDigitize: YVK : iter="<<iter<<endl;
fa412d9b 795
796 Float_t zhit = zxe[0][iter];
797 Float_t xhit = zxe[1][iter];
798 Float_t qhit = zxe[2][iter];
799 Float_t zcell = zhit / fGeom->GetPadSizeZ();
800 Float_t xcell = xhit / fGeom->GetPadSizePhi();
801 if ( zcell<=0 || xcell<=0 ||
802 zcell>=nCellZ || xcell>=nCellX) return;
803 Int_t izcell = (Int_t) zcell;
804 Int_t ixcell = (Int_t) xcell;
805 Float_t zc = zcell - izcell - 0.5;
806 Float_t xc = xcell - ixcell - 0.5;
807 for (Int_t iz=1; iz<=ngamz; iz++) {
808 Int_t kzg = izcell + iz - nz3;
809 if (kzg<=0 || kzg>nCellZ) continue;
810 Float_t zg = (Float_t)(iz-nz3) - zc;
811 for (Int_t ix=1; ix<=ngamx; ix++) {
812 Int_t kxg = ixcell + ix - nx3;
813 if (kxg<=0 || kxg>nCellX) continue;
814 Float_t xg = (Float_t)(ix-nx3) - xc;
815
816 // Now calculate pad response
817 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
818 qpad += qNoise*rnor2;
819 if (qpad<0) continue;
820
821 // Fill the array with pad response ID and amplitude
3d402178 822 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
fa412d9b 823// printf("(%2d,%2d,%5.3f) ",kxg,kzg,qpad);
824 }
825// cout << endl;
826 }
827// cout << endl;
828 }
829}
830
831//____________________________________________________________________________
832Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
833 // ------------------------------------------------------------------------
834 // Calculate the amplitude in one CPV pad using the
835 // cumulative pad response function
836 // Author: Yuri Kharlov (after Serguei Sadovski)
837 // 3 October 2000
838 // ------------------------------------------------------------------------
839
840 Double_t dz = fGeom->GetPadSizeZ() / 2;
841 Double_t dx = fGeom->GetPadSizePhi() / 2;
842 Double_t z = zhit * fGeom->GetPadSizeZ();
843 Double_t x = xhit * fGeom->GetPadSizePhi();
844 Double_t amplitude = qhit *
845 (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
846 CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
847 return (Float_t)amplitude;
7587f5a5 848}
849
fa412d9b 850//____________________________________________________________________________
851Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
852 // ------------------------------------------------------------------------
853 // Cumulative pad response function
854 // It includes several terms from the CF decomposition in electrostatics
855 // Note: this cumulative function is wrong since omits some terms
856 // but the cell amplitude obtained with it is correct because
857 // these omitting terms cancel
858 // Author: Yuri Kharlov (after Serguei Sadovski)
859 // 3 October 2000
860 // ------------------------------------------------------------------------
861
862 const Double_t a=1.0;
863 const Double_t b=0.7;
864
865 Double_t r2 = x*x + y*y;
866 Double_t xy = x*y;
867 Double_t cumulPRF = 0;
868 for (Int_t i=0; i<=4; i++) {
869 Double_t b1 = (2*i + 1) * b;
870 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
871 }
872 cumulPRF *= a/(2*TMath::Pi());
873 return cumulPRF;
874}