parameters have been redistributed; Hits2SDigits etc ... introduce
[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
a3dfe79c 20//---
21// Layout EMC + PPSD has name GPS2:
ed4205d8 22// Produces cumulated hits
a3dfe79c 23//---
24// Layout EMC + CPV has name IHEP:
ed4205d8 25// Produces hits for CPV, cumulated hits
26//---
27// Layout EMC + CPV + PPSD has name GPS:
28// Produces hits for CPV, cumulated hits
29//---
5f20d3fb 30//*-- Author: Yves Schutz (SUBATECH)
b2a60966 31
7587f5a5 32
33// --- ROOT system ---
bea63bea 34
35#include "TBRIK.h"
36#include "TNode.h"
7587f5a5 37#include "TRandom.h"
94de3818 38#include "TTree.h"
7587f5a5 39
5f20d3fb 40
7587f5a5 41// --- Standard library ---
42
de9ec31b 43#include <stdio.h>
44#include <string.h>
45#include <stdlib.h>
46#include <strstream.h>
7587f5a5 47
48// --- AliRoot header files ---
49
50#include "AliPHOSv1.h"
51#include "AliPHOSHit.h"
52#include "AliPHOSDigit.h"
bea63bea 53#include "AliPHOSReconstructioner.h"
7587f5a5 54#include "AliRun.h"
55#include "AliConst.h"
94de3818 56#include "AliMC.h"
7587f5a5 57
58ClassImp(AliPHOSv1)
59
bea63bea 60//____________________________________________________________________________
61AliPHOSv1::AliPHOSv1()
62{
5f20d3fb 63 // ctor
fc879520 64
3d402178 65 // Create an empty array of AliPHOSCPVModule to satisfy
fa412d9b 66 // AliPHOSv1::Streamer when reading root file
ed4205d8 67
68 fReconstructioner = 0;
69 fTrackSegmentMaker = 0;
037cc66d 70
71 fHits = new TClonesArray("AliPHOSHit",1000) ;
ed4205d8 72
037cc66d 73 // if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
74 // Error("AliPHOSv1","Can not create array of EMC modules");
75 // exit(1);
76 // }
fa412d9b 77
037cc66d 78 // if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",0)) ) {
79 // Error("AliPHOSv1","Can not create array of CPV modules");
80 // exit(1);
81 // }
fc879520 82
bea63bea 83}
84
7587f5a5 85//____________________________________________________________________________
86AliPHOSv1::AliPHOSv1(const char *name, const char *title):
e04976bd 87AliPHOSv0(name,title)
7587f5a5 88{
5f20d3fb 89 // ctor : title is used to identify the layout
fa412d9b 90 // GPS2 = 5 modules (EMC + PPSD)
91 // IHEP = 5 modules (EMC + CPV )
ed4205d8 92 // MIXT = 4 modules (EMC + CPV ) and 1 module (EMC + PPSD)
5f20d3fb 93 //
ed4205d8 94 // We store hits :
5f20d3fb 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).
98 //
fa412d9b 99
5f20d3fb 100 fPinElectronicNoise = 0.010 ;
037cc66d 101 fDigitThreshold = 0.01 ; // 1 GeV
102 fDigitizeA= 0. ;
103 fDigitizeB = 10000000. ;
104
5f20d3fb 105
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
fa412d9b 109 // and the TreeD at the end of the event (branch is set in FinishEvent() ).
5f20d3fb 110
ed4205d8 111 fHits= new TClonesArray("AliPHOSHit",1000) ;
5f20d3fb 112
ed4205d8 113 fNhits = 0 ;
5f20d3fb 114
ed4205d8 115 fReconstructioner = 0;
116 fTrackSegmentMaker = 0;
5f20d3fb 117
5f20d3fb 118 fIshunt = 1 ; // All hits are associated with primary particles
fa412d9b 119
ed4205d8 120 // Create array of EMC modules of the size of PHOS modules number
121
037cc66d 122// if ( 0==(fEMCModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNModules())) ) {
123// Error("AliPHOSv1","Can not create array of EMC modules");
124// exit(1);
125// }
126// TClonesArray &lemcmodule = *fEMCModules;
127// for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lemcmodule[i]) AliPHOSCPVModule();
ed4205d8 128
fa412d9b 129 // Create array of CPV modules for the IHEP's version of CPV
130
037cc66d 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
fa412d9b 133
037cc66d 134// if ( 0==(fCPVModules=new TClonesArray("AliPHOSCPVModule",fGeom->GetNCPVModules())) ) {
135// Error("AliPHOSv1","Can not create array of CPV modules");
136// exit(1);
137// }
138// TClonesArray &lcpvmodule = *fCPVModules;
139// for (Int_t i=0; i<fGeom->GetNCPVModules(); i++) new(lcpvmodule[i]) AliPHOSCPVModule();
140// }
141// else {
142// // Create an empty array of AliPHOSCPVModule to satisfy
143// // AliPHOSv1::Streamer when writing root file
fa412d9b 144
037cc66d 145// fCPVModules=new TClonesArray("AliPHOSCPVModule",0);
fa412d9b 146
037cc66d 147// }
5f20d3fb 148}
149
150//____________________________________________________________________________
151AliPHOSv1::AliPHOSv1(AliPHOSReconstructioner * Reconstructioner, const char *name, const char *title):
152 AliPHOSv0(name,title)
153{
154 // ctor : title is used to identify the layout
155 // GPS2 = 5 modules (EMC + PPSD)
156 // We use 2 arrays of hits :
157 //
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).
161 //
162 // - fTmpHits, which retains all the hits of the current event. It
163 // is used for the digitization part.
164
165 fPinElectronicNoise = 0.010 ;
166
167 // We do not want to save in TreeH the raw hits
168 //fHits = new TClonesArray("AliPHOSHit",100) ;
169
ed4205d8 170 fDigits = 0 ;
171 fHits= new TClonesArray("AliPHOSHit",1000) ;
b73f246d 172
ed4205d8 173 fNhits = 0 ;
5f20d3fb 174
175 fIshunt = 1 ; // All hits are associated with primary particles
176
177 // gets an instance of the geometry parameters class
178 fGeom = AliPHOSGeometry::GetInstance(title, "") ;
179
180 if (fGeom->IsInitialized() )
88bdfa12 181 cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
5f20d3fb 182 else
fc879520 183 cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;
5f20d3fb 184
185 // Defining the PHOS Reconstructioner
186
187 fReconstructioner = Reconstructioner ;
188
7587f5a5 189}
b2a60966 190
7587f5a5 191//____________________________________________________________________________
bea63bea 192AliPHOSv1::~AliPHOSv1()
b2a60966 193{
bea63bea 194 // dtor
5f20d3fb 195
ed4205d8 196 if ( fHits) {
197 fHits->Delete() ;
198 delete fHits ;
199 fHits = 0 ;
8dfa469d 200 }
5f20d3fb 201
037cc66d 202 if ( fSDigits) {
203 fSDigits->Delete() ;
204 delete fSDigits ;
205 fSDigits = 0 ;
206 }
207
ed4205d8 208 if ( fDigits) {
209 fDigits->Delete() ;
210 delete fDigits ;
211 fDigits = 0 ;
212 }
5f20d3fb 213
ed4205d8 214 if ( fEmcRecPoints ) {
215 fEmcRecPoints->Delete() ;
216 delete fEmcRecPoints ;
217 fEmcRecPoints = 0 ;
218 }
219
220 if ( fPpsdRecPoints ) {
221 fPpsdRecPoints->Delete() ;
222 delete fPpsdRecPoints ;
223 fPpsdRecPoints = 0 ;
224 }
fad3e5b9 225
ed4205d8 226 if ( fTrackSegments ) {
227 fTrackSegments->Delete() ;
228 delete fTrackSegments ;
229 fTrackSegments = 0 ;
230 }
231
7587f5a5 232}
233
7587f5a5 234//____________________________________________________________________________
037cc66d 235void 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)
bea63bea 236{
237 // Add a hit to the hit list.
5f20d3fb 238 // A PHOS hit is the sum of all hits in a single crystal
239 // or in a single PPSD gas cell
bea63bea 240
5f20d3fb 241 Int_t hitCounter ;
bea63bea 242 AliPHOSHit *newHit ;
5f20d3fb 243 AliPHOSHit *curHit ;
244 Bool_t deja = kFALSE ;
bea63bea 245
037cc66d 246 newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits, trackpid, p, lpos) ;
bea63bea 247
ed4205d8 248 for ( hitCounter = 0 ; hitCounter < fNhits && !deja ; hitCounter++ ) {
249 curHit = (AliPHOSHit*) (*fHits)[hitCounter] ;
250 if( *curHit == *newHit ) {
251 *curHit = *curHit + *newHit ;
252 deja = kTRUE ;
5f20d3fb 253 }
254 }
255
256 if ( !deja ) {
ed4205d8 257 new((*fHits)[fNhits]) AliPHOSHit(*newHit) ;
258 fNhits++ ;
5f20d3fb 259 }
260
bea63bea 261 delete newHit;
bea63bea 262}
263
ed4205d8 264//____________________________________________________________________________
037cc66d 265void AliPHOSv1::Hits2SDigits(){
266 //Collects all hits in the same active volume into digit
ed4205d8 267
bea63bea 268 Int_t i ;
bea63bea 269 Int_t j ;
fa412d9b 270 AliPHOSHit * hit ;
bea63bea 271 AliPHOSDigit * newdigit ;
272 AliPHOSDigit * curdigit ;
273 Bool_t deja = kFALSE ;
274
037cc66d 275
ed4205d8 276 Int_t itrack ;
037cc66d 277 for (itrack=0; itrack<gAlice->GetNtrack(); itrack++){
278
ed4205d8 279 //=========== Get the Hits Tree for the Primary track itrack
280 gAlice->ResetHits();
281 gAlice->TreeH()->GetEvent(itrack);
282
037cc66d 283
ed4205d8 284 for ( i = 0 ; i < fHits->GetEntries() ; i++ ) {
285 hit = (AliPHOSHit*)fHits->At(i) ;
286
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() ) ) ;
290 else
291 newdigit = new AliPHOSDigit( -1 , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
292 deja =kFALSE ;
293
037cc66d 294 for ( j = 0 ; j < fnSdigits ; j++) {
295 curdigit = (AliPHOSDigit*) fSDigits->At(j) ;
ed4205d8 296 if ( *curdigit == *newdigit) {
297 *curdigit = *curdigit + *newdigit ;
298 deja = kTRUE ;
299 }
bea63bea 300 }
bea63bea 301
ed4205d8 302 if ( !deja ) {
037cc66d 303 new((*fSDigits)[fnSdigits]) AliPHOSDigit(* newdigit) ;
304 fnSdigits++ ;
ed4205d8 305 }
306
307 delete newdigit ;
308 }
309
310 } // loop over tracks
037cc66d 311
312 fSDigits->Sort() ;
313
314 fnSdigits = fSDigits->GetEntries() ;
315 fSDigits->Expand(fnSdigits) ;
316
317 for (i = 0 ; i < fnSdigits ; i++) {
318 AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(i) ;
319 digit->SetIndexInList(i) ;
320 }
321
322 gAlice->TreeS()->Fill() ;
ed4205d8 323
037cc66d 324 TTree *ts = gAlice->TreeS();
325 TDirectory *dir = ts->GetDirectory() ;
326 cout << "dir name is " << dir->GetName() << endl ;
327 // ts->Write(0,TObject::kOverwrite) ;
328
329gAlice->TreeS()->GetBranch("PHOS")->Fill();
330gAlice->TreeS()->GetBranch("PHOS")->Write();
331
332
333}
334//____________________________________________________________________________
335void 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.
338
339
340 gAlice->TreeS()->GetEvent(0) ;
341
342 cout << "fSdigits " << fSDigits << " " << fSDigits->GetEntries() << endl ;
343
344
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() ;
349
350 cout << "fDigits " << fDigits << " " <<idCurSDigit<< endl ;
351
352 Int_t absID ;
353 for(absID = 1; absID < fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ(); absID++){
ed4205d8 354
037cc66d 355 cout << "absID " << absID << " " << idCurSDigit << endl ;
356
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 ;
364 fNdigits++ ;
365 }
366 }
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 ;
377 fNdigits++ ;
378 }
379
380 if(iCurSDigit < fSDigits->GetEntries()-1){
381 iCurSDigit++ ;
382 idCurSDigit = ((AliPHOSDigit*)fSDigits->At(iCurSDigit))->GetId() ;
383 }
384 else
385 idCurSDigit = 10000000; //no real hits left
386 }
ed4205d8 387
037cc66d 388 }
389
390 //remove PPSD/CPV digits below thresholds
391 Int_t idigit ;
392 for(idigit = iCurSDigit; idigit < fSDigits->GetEntries() ; idigit++){ //loop over CPV/PPSD digits
393
394 AliPHOSDigit * digit = (AliPHOSDigit *) fSDigits->At(idigit) ;
395 Float_t ene = Calibrate(digit->GetAmp()) ;
396
397 Int_t relid[4] ;
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 ) ;
403 fNdigits++ ;
bea63bea 404 }
037cc66d 405 }
ed4205d8 406
bea63bea 407 fDigits->Compress() ;
ed4205d8 408
409 fNdigits = fDigits->GetEntries() ;
410 fDigits->Expand(fNdigits) ;
bea63bea 411
037cc66d 412 Int_t i ;
bea63bea 413 for (i = 0 ; i < fNdigits ; i++) {
037cc66d 414 AliPHOSDigit * digit = (AliPHOSDigit *) fDigits->At(i) ;
415 digit->SetIndexInList(i) ;
bea63bea 416 }
fa412d9b 417
037cc66d 418
ed4205d8 419 gAlice->TreeD()->Fill() ;
5f20d3fb 420
ed4205d8 421 gAlice->TreeD()->Write(0,TObject::kOverwrite) ;
422
423}
2ab0c725 424
5f20d3fb 425//___________________________________________________________________________
2ab0c725 426void AliPHOSv1::MakeBranch(Option_t* opt, char *file)
037cc66d 427{
428
429
430 char *cH ;
5f20d3fb 431 // Create new branche in the current Root Tree in the digit Tree
5f20d3fb 432 AliDetector::MakeBranch(opt) ;
037cc66d 433
434
435 cH = strstr(opt,"S");
436 //Create a branch for SDigits
437 if( cH ){
438 char branchname[20];
439 sprintf(branchname,"%s",GetName());
440 if(fSDigits)
441 fSDigits->Clear();
442 else
443 fSDigits = new TClonesArray("AliPHOSDigit",1000);
444 fnSdigits = 0 ;
445 cout << " AliPHOSv1::MakeBranch : " << file << endl ;
446 gAlice->MakeBranchInTree(gAlice->TreeS(),branchname,&fSDigits,fBufferSize,file);
fa412d9b 447 }
ed4205d8 448
037cc66d 449 cH = strstr(opt,"D");
450 //Create a branch for Digits
451 if( cH ){
452 char branchname[20];
453 sprintf(branchname,"%s",GetName());
454 if(fSDigits)
455 fDigits->Clear();
456 else
457 fDigits = new TClonesArray("AliPHOSDigit",1000);
458 fNdigits = 0 ;
459
460 gAlice->MakeBranchInTree(gAlice->TreeD(),branchname,&fSDigits,fBufferSize,file);
461 }
462
463 cH = strstr(opt,"R");
464 //Create a branch for Reconstruction
465 if( cH ){
466 char branchname[20];
467
468 Int_t splitlevel = 0 ;
469
470 fEmcRecPoints->Delete() ;
471
472 if ( fEmcRecPoints && gAlice->TreeR() ) {
473 sprintf(branchname,"%sEmcRP",GetName()) ;
474 gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ;
475 }
476
477 fPpsdRecPoints->Delete() ;
478
479 if ( fPpsdRecPoints && gAlice->TreeR() ) {
480 sprintf(branchname,"%sPpsdRP",GetName()) ;
481 gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
482 }
483
484 fTrackSegments->Delete() ;
485
486 if ( fTrackSegments && gAlice->TreeR() ) {
487 sprintf(branchname,"%sTS",GetName()) ;
488 gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
489 }
490
491 fRecParticles->Delete() ;
492
493 if ( fRecParticles && gAlice->TreeR() ) {
494 sprintf(branchname,"%sRP",GetName()) ;
495 gAlice->TreeR()->Branch(branchname, &fRecParticles, fBufferSize) ;
496 }
497
498 }
499
500
501
bea63bea 502}
503
5f20d3fb 504//_____________________________________________________________________________
505void AliPHOSv1::Reconstruction(AliPHOSReconstructioner * Reconstructioner)
506{
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
512
513 fReconstructioner = Reconstructioner ;
037cc66d 514
5f20d3fb 515 // 1.
516
517 // gAlice->MakeTree("R") ;
5f20d3fb 518
037cc66d 519 MakeBranch("R") ;
5f20d3fb 520
521 // 3.
b73f246d 522
fad3e5b9 523 fReconstructioner->Make(fDigits, fEmcRecPoints, fPpsdRecPoints, fTrackSegments, fRecParticles);
b73f246d 524
fad3e5b9 525 printf("Reconstruction: %d %d %d %d\n",
526 fEmcRecPoints->GetEntries(),fPpsdRecPoints->GetEntries(),
b73f246d 527 fTrackSegments->GetEntries(),fRecParticles->GetEntries());
5f20d3fb 528
529 // 4. Expand or Shrink the arrays to the proper size
530
531 Int_t size ;
532
533 size = fEmcRecPoints->GetEntries() ;
534 fEmcRecPoints->Expand(size) ;
ed4205d8 535
5f20d3fb 536 size = fPpsdRecPoints->GetEntries() ;
537 fPpsdRecPoints->Expand(size) ;
538
539 size = fTrackSegments->GetEntries() ;
540 fTrackSegments->Expand(size) ;
541
542 size = fRecParticles->GetEntries() ;
543 fRecParticles->Expand(size) ;
544
545 gAlice->TreeR()->Fill() ;
5f20d3fb 546 // 5.
547
4a2ca5e9 548 gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
5f20d3fb 549
550 // Deleting reconstructed objects
551 ResetReconstruction();
552
553}
554
fa412d9b 555//____________________________________________________________________________
556void AliPHOSv1::ResetHits()
c4f224e7 557{
ed4205d8 558 // Reset hit tree for EMC and CPV
fa412d9b 559 // Yuri Kharlov, 28 September 2000
560
561 AliDetector::ResetHits();
037cc66d 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();
565 // }
ed4205d8 566
5f20d3fb 567}
568//____________________________________________________________________________
569void AliPHOSv1::ResetReconstruction()
570{
571 // Deleting reconstructed objects
572
b73f246d 573 if ( fEmcRecPoints ) fEmcRecPoints ->Delete();
5f20d3fb 574 if ( fPpsdRecPoints ) fPpsdRecPoints->Delete();
575 if ( fTrackSegments ) fTrackSegments->Delete();
b73f246d 576 if ( fRecParticles ) fRecParticles ->Delete();
5f20d3fb 577
578}
5f20d3fb 579
037cc66d 580//____________________________________________________________________________
581//void AliPHOSv1::SDigits2Digits() {
582// // Adds the noise to the summable digits and keeps digits above a threshold
583// // To make a digit
584//}
585
5f20d3fb 586//____________________________________________________________________________
587void AliPHOSv1::SetTreeAddress()
588{
589 // TBranch *branch;
590 AliPHOS::SetTreeAddress();
591
ed4205d8 592// //Branch address for TreeR: RecPpsdRecPoint
5f20d3fb 593// TTree *treeR = gAlice->TreeR();
594// if ( treeR && fPpsdRecPoints ) {
595// branch = treeR->GetBranch("PHOSPpsdRP");
596// if (branch) branch->SetAddress(&fPpsdRecPoints) ;
fa412d9b 597// }
598
ed4205d8 599 // Set branch address for the Hits Tree for hits in EMC modules
600 // Yuri Kharlov, 23 November 2000.
601
037cc66d 602 // for( Int_t i=0; i<fGeom->GetNModules(); i++ ) GetEMCModule(i).SetTreeAddress("EMC",i+1);
ed4205d8 603
fa412d9b 604 // Set branch address for the Hits Tree for hits in CPV modules for IHEP geometry
605 // Yuri Kharlov, 28 September 2000.
606
037cc66d 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);
609 // }
fa412d9b 610
5f20d3fb 611}
612
613//____________________________________________________________________________
614
7587f5a5 615void AliPHOSv1::StepManager(void)
616{
b2a60966 617 // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
b2a60966 618
ed4205d8 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;
623
7587f5a5 624 Int_t relid[4] ; // (box, layer, row, column) indices
fa412d9b 625 Int_t absid ; // absolute cell ID number
2168f43b 626 Float_t xyze[4]={0,0,0,0} ; // position wrt MRS and energy deposited
fa412d9b 627 TLorentzVector pos ; // Lorentz vector of the track current position
037cc66d 628 TLorentzVector pmom ; //momentum of the particle initiated hit
629 Float_t xyd[2] ; //local posiiton of the entering
630 Bool_t entered = kFALSE ;
fa412d9b 631 Int_t copy ;
7587f5a5 632
bea63bea 633 Int_t tracknumber = gAlice->CurrentTrack() ;
fa412d9b 634 Int_t primary = gAlice->GetPrimary( gAlice->CurrentTrack() );
635 TString name = fGeom->GetName() ;
037cc66d 636 Int_t trackpid = 0 ;
637
638 if( gMC->IsTrackEntering() ){ // create hit with position and momentum of new particle,
639 // but may be without energy deposition
640
641 // Current position of the hit in the local ref. system
642 gMC -> TrackPosition(pos);
643 Float_t xyzm[3], xyzd[3] ;
644 Int_t i;
645 for (i=0; i<3; i++) xyzm[i] = pos[i];
646 gMC -> Gmtod (xyzm, xyzd, 1); // transform coordinate from master to daughter system
647 xyd[0] = xyzd[0];
648 xyd[1] =-xyzd[2];
649
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
655 pmom[0] = pd[0];
656 pmom[1] =-pd[1];
657 pmom[2] =-pd[2];
658
659 trackpid = gMC->TrackPid();
660 entered = kTRUE ; // Mark to create hit even withou energy deposition
661
662 }
663
664
ed4205d8 665 if ( name == "GPS2" || name == "MIXT" ) { // ======> CPV is a GPS' PPSD
fa412d9b 666
b2a60966 667 if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell
7587f5a5 668 {
669 gMC->TrackPosition(pos) ;
670 xyze[0] = pos[0] ;
671 xyze[1] = pos[1] ;
672 xyze[2] = pos[2] ;
bea63bea 673 xyze[3] = gMC->Edep() ;
7587f5a5 674
037cc66d 675 if ( (xyze[3] != 0) || entered ) { // there is deposited energy or new particle entering PPSD
7587f5a5 676 gMC->CurrentVolOffID(5, relid[0]) ; // get the PHOS Module number
fad3e5b9 677 if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(5),"PHO1") == 0 ){
ed4205d8 678 relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
fad3e5b9 679 }
7587f5a5 680 gMC->CurrentVolOffID(3, relid[1]) ; // get the Micromegas Module number
ed4205d8 681 // 1-> fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() upper
682 // > fGeom->GetNumberOfModulesPhi() * fGeom->GetNumberOfModulesZ() lower
7587f5a5 683 gMC->CurrentVolOffID(1, relid[2]) ; // get the row number of the cell
684 gMC->CurrentVolID(relid[3]) ; // get the column number
685
686 // get the absolute Id number
687
bea63bea 688 fGeom->RelToAbsNumbering(relid, absid) ;
7587f5a5 689
bea63bea 690 // add current hit to the hit list
037cc66d 691 AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
692
7587f5a5 693
694 } // there is deposited energy
fa412d9b 695 } // We are inside the gas of the CPV
696 } // GPS2 configuration
697
ed4205d8 698 if ( name == "IHEP" || name == "MIXT" ) { // ======> CPV is a IHEP's one
fa412d9b 699
700 // Yuri Kharlov, 28 September 2000
701
702 if( gMC->CurrentVolID(copy) == gMC->VolId("CPVQ") &&
037cc66d 703 entered &&
704 gMC->TrackCharge() != 0) {
fa412d9b 705
037cc66d 706 // Digitize the current CPV hit:
707
708 // 1. find pad response and
fa412d9b 709
a3dfe79c 710 Int_t moduleNumber;
711 gMC->CurrentVolOffID(3,moduleNumber);
712 moduleNumber--;
fa412d9b 713
fa412d9b 714
3d402178 715 TClonesArray *cpvDigits = new TClonesArray("AliPHOSCPVDigit",0); // array of digits for current hit
a3dfe79c 716 CPVDigitize(pmom,xyd,moduleNumber,cpvDigits);
fa412d9b 717
718 Float_t xmean = 0;
719 Float_t zmean = 0;
720 Float_t qsum = 0;
cd461ab8 721 Int_t idigit,ndigits;
fa412d9b 722
723 // 2. go through the current digit list and sum digits in pads
724
725 ndigits = cpvDigits->GetEntriesFast();
cd461ab8 726 for (idigit=0; idigit<ndigits-1; idigit++) {
3d402178 727 AliPHOSCPVDigit *cpvDigit1 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
fa412d9b 728 Float_t x1 = cpvDigit1->GetXpad() ;
729 Float_t z1 = cpvDigit1->GetYpad() ;
730 for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
3d402178 731 AliPHOSCPVDigit *cpvDigit2 = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(jdigit);
fa412d9b 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) ;
738 }
739 }
740 }
741 cpvDigits->Compress() ;
742
743 // 3. add digits to temporary hit list fTmpHits
744
745 ndigits = cpvDigits->GetEntriesFast();
cd461ab8 746 for (idigit=0; idigit<ndigits; idigit++) {
3d402178 747 AliPHOSCPVDigit *cpvDigit = (AliPHOSCPVDigit*) cpvDigits->UncheckedAt(idigit);
a3dfe79c 748 relid[0] = moduleNumber + 1 ; // CPV (or PHOS) module number
fa412d9b 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
752
753 // get the absolute Id number
754 fGeom->RelToAbsNumbering(relid, absid) ;
755
756 // add current digit to the temporary hit list
757 xyze[0] = 0. ;
758 xyze[1] = 0. ;
759 xyze[2] = 0. ;
760 xyze[3] = cpvDigit->GetQpad() ; // amplitude in a pad
761 primary = -1; // No need in primary for CPV
037cc66d 762 AddHit(fIshunt, primary, tracknumber, absid, xyze, trackpid, pmom, xyd);
fa412d9b 763
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();
768 }
769 }
770 delete cpvDigits;
771 }
772 } // end of IHEP configuration
7587f5a5 773
037cc66d 774
fa412d9b 775 if(gMC->CurrentVolID(copy) == gMC->VolId("PXTL") ) { // We are inside a PBWO crystal
776 gMC->TrackPosition(pos) ;
777 xyze[0] = pos[0] ;
778 xyze[1] = pos[1] ;
779 xyze[2] = pos[2] ;
780 xyze[3] = gMC->Edep() ;
ed4205d8 781
037cc66d 782
783 if ( (xyze[3] != 0) || entered ) { // Track is inside the crystal and deposits some energy or just entered
ed4205d8 784
fa412d9b 785 gMC->CurrentVolOffID(10, relid[0]) ; // get the PHOS module number ;
037cc66d 786
fad3e5b9 787 if ( name == "MIXT" && strcmp(gMC->CurrentVolOffName(10),"PHO1") == 0 )
788 relid[0] += fGeom->GetNModules() - fGeom->GetNPPSDModules();
789
fa412d9b 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
793
794 // get the absolute Id number
fa412d9b 795 fGeom->RelToAbsNumbering(relid, absid) ;
ed4205d8 796
fa412d9b 797 // add current hit to the hit list
037cc66d 798 AddHit(fIshunt, primary,tracknumber, absid, xyze, trackpid,pmom, xyd);
799
ed4205d8 800
fa412d9b 801 } // there is deposited energy
802 } // we are inside a PHOS Xtal
037cc66d 803
804
fa412d9b 805}
806
807//____________________________________________________________________________
808void AliPHOSv1::CPVDigitize (TLorentzVector p, Float_t *zxhit, Int_t moduleNumber, TClonesArray *cpvDigits)
809{
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
815 //
816 // Author: Yuri Kharlov (after Serguei Sadovsky)
817 // 2 October 2000
818 // ------------------------------------------------------------------------
819
a3dfe79c 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
fa412d9b 826
827 Float_t rnor1,rnor2;
828
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
833
834 Float_t hitX = zxhit[0];
835 Float_t hitZ =-zxhit[1];
836 Float_t pX = p.Px();
837 Float_t pZ =-p.Pz();
838 Float_t pNorm = p.Py();
a3dfe79c 839 Float_t eloss = kdEdx;
3d402178 840
7eb9d12d 841// cout << "CPVDigitize: YVK : "<<hitX<<" "<<hitZ<<" | "<<pX<<" "<<pZ<<" "<<pNorm<<endl;
842
fa412d9b 843 Float_t dZY = pZ/pNorm * fGeom->GetCPVGasThickness();
844 Float_t dXY = pX/pNorm * fGeom->GetCPVGasThickness();
845 gRandom->Rannor(rnor1,rnor2);
a3dfe79c 846 eloss *= (1 + kDetR*rnor1) *
847 TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
fa412d9b 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;
852
a3dfe79c 853 Int_t iwht1 = (Int_t) (xhit1 / kCelWr); // wire (x) coordinate "in"
854 Int_t iwht2 = (Int_t) (xhit2 / kCelWr); // wire (x) coordinate "out"
fa412d9b 855
856 Int_t nIter;
857 Float_t zxe[3][5];
858 if (iwht1==iwht2) { // incline 1-wire hit
859 nIter = 2;
860 zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
a3dfe79c 861 zxe[1][0] = (iwht1 + 0.5) * kCelWr;
862 zxe[2][0] = eloss/2;
fa412d9b 863 zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
a3dfe79c 864 zxe[1][1] = (iwht1 + 0.5) * kCelWr;
865 zxe[2][1] = eloss/2;
fa412d9b 866 }
867 else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
868 nIter = 3;
869 Int_t iwht3 = (iwht1 + iwht2) / 2;
a3dfe79c 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
fa412d9b 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;
a3dfe79c 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 );
fa412d9b 880 zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
881 zxe[1][0] = xwht1;
a3dfe79c 882 zxe[2][0] = eloss * egm1;
fa412d9b 883 zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
884 zxe[1][1] = xwht2;
a3dfe79c 885 zxe[2][1] = eloss * egm2;
fa412d9b 886 zxe[0][2] = dXY*(xwht3-xwht1)/dXY + zhit1;
887 zxe[1][2] = xwht3;
a3dfe79c 888 zxe[2][2] = eloss * egm3;
fa412d9b 889 }
890 else { // incline 2-wire hit
891 nIter = 2;
a3dfe79c 892 Float_t xwht1 = (iwht1 + 0.5) * kCelWr;
893 Float_t xwht2 = (iwht2 + 0.5) * kCelWr;
fa412d9b 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;
900 zxe[1][0] = xwht1;
a3dfe79c 901 zxe[2][0] = eloss * egm1;
fa412d9b 902 zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
903 zxe[1][1] = xwht2;
a3dfe79c 904 zxe[2][1] = eloss * egm2;
fa412d9b 905 }
bea63bea 906
fa412d9b 907 // Finite size of ionization region
908
ed4205d8 909 Int_t nCellZ = fGeom->GetNumberOfCPVPadsZ();
910 Int_t nCellX = fGeom->GetNumberOfCPVPadsPhi();
a3dfe79c 911 Int_t nz3 = (kNgamz+1)/2;
912 Int_t nx3 = (kNgamx+1)/2;
913 cpvDigits->Expand(nIter*kNgamx*kNgamz);
fa412d9b 914 TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
915
916 for (Int_t iter=0; iter<nIter; iter++) {
917
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;
a3dfe79c 929 for (Int_t iz=1; iz<=kNgamz; iz++) {
fa412d9b 930 Int_t kzg = izcell + iz - nz3;
931 if (kzg<=0 || kzg>nCellZ) continue;
932 Float_t zg = (Float_t)(iz-nz3) - zc;
a3dfe79c 933 for (Int_t ix=1; ix<=kNgamx; ix++) {
fa412d9b 934 Int_t kxg = ixcell + ix - nx3;
935 if (kxg<=0 || kxg>nCellX) continue;
936 Float_t xg = (Float_t)(ix-nx3) - xc;
937
938 // Now calculate pad response
939 Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
a3dfe79c 940 qpad += kNoise*rnor2;
fa412d9b 941 if (qpad<0) continue;
942
943 // Fill the array with pad response ID and amplitude
3d402178 944 new(ldigits[cpvDigits->GetEntriesFast()]) AliPHOSCPVDigit(kxg,kzg,qpad);
fa412d9b 945 }
fa412d9b 946 }
fa412d9b 947 }
948}
949
950//____________________________________________________________________________
951Float_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)
956 // 3 October 2000
957 // ------------------------------------------------------------------------
958
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;
7587f5a5 967}
968
fa412d9b 969//____________________________________________________________________________
970Double_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)
978 // 3 October 2000
979 // ------------------------------------------------------------------------
980
a3dfe79c 981 const Double_t kA=1.0;
982 const Double_t kB=0.7;
fa412d9b 983
984 Double_t r2 = x*x + y*y;
985 Double_t xy = x*y;
986 Double_t cumulPRF = 0;
987 for (Int_t i=0; i<=4; i++) {
a3dfe79c 988 Double_t b1 = (2*i + 1) * kB;
fa412d9b 989 cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
990 }
a3dfe79c 991 cumulPRF *= kA/(2*TMath::Pi());
fa412d9b 992 return cumulPRF;
993}
7eb9d12d 994