]>
Commit | Line | Data |
---|---|---|
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 | |
50 | ClassImp(AliPHOSv1) | |
51 | ||
bea63bea | 52 | //____________________________________________________________________________ |
53 | AliPHOSv1::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 | //____________________________________________________________________________ |
70 | AliPHOSv1::AliPHOSv1(const char *name, const char *title): | |
e04976bd | 71 | AliPHOSv0(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 | //____________________________________________________________________________ | |
123 | AliPHOSv1::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 | 164 | AliPHOSv1::~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 | 195 | void 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 <mphits = *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 | //___________________________________________________________________________ |
240 | Int_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 | //___________________________________________________________________________ | |
251 | void 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 | //___________________________________________________________________________ | |
329 | void 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 | //_____________________________________________________________________________ |
360 | void 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 | //____________________________________________________________________________ |
443 | void 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 | //____________________________________________________________________________ |
454 | void 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 | //____________________________________________________________________________ | |
464 | void 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 | //____________________________________________________________________________ | |
476 | void 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 | 499 | void 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 | //____________________________________________________________________________ | |
680 | void 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 | //____________________________________________________________________________ | |
832 | Float_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 | //____________________________________________________________________________ |
851 | Double_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 | } |