Updated StepManager(), new functions for CPV
[u/mrichter/AliRoot.git] / PHOS / AliPHOSv1.cxx
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
16 /* $Id$ */
17
18 //_________________________________________________________________________
19 // Implementation version v1 of PHOS Manager class 
20 // Layout EMC + PPSD has name GPS2  
21 // Produces cumulated hits (no hits) and digits                  
22 //*-- Author: Yves Schutz (SUBATECH)
23
24
25 // --- ROOT system ---
26
27 #include "TBRIK.h"
28 #include "TNode.h"
29 #include "TRandom.h"
30 #include "TTree.h"
31
32
33 // --- Standard library ---
34
35 #include <stdio.h>
36 #include <string.h>
37 #include <stdlib.h>
38 #include <strstream.h>
39
40 // --- AliRoot header files ---
41
42 #include "AliPHOSv1.h"
43 #include "AliPHOSHit.h"
44 #include "AliPHOSDigit.h"
45 #include "AliPHOSReconstructioner.h"
46 #include "AliRun.h"
47 #include "AliConst.h"
48 #include "AliMC.h"
49
50 ClassImp(AliPHOSv1)
51
52 //____________________________________________________________________________
53 AliPHOSv1::AliPHOSv1()
54 {
55   // ctor
56   fNTmpHits = 0 ; 
57   fTmpHits  = 0 ;
58
59   // Create an empty array of CPVModule to satisfy
60   // AliPHOSv1::Streamer when reading root file
61
62   if ( NULL==(fCPVModules=new TClonesArray("CPVModule",0)) ) {
63     Error("AliPHOSv1","Can not create array of CPV modules");
64     exit(1);
65   }
66
67 }
68
69 //____________________________________________________________________________
70 AliPHOSv1::AliPHOSv1(const char *name, const char *title):
71 AliPHOSv0(name,title) 
72 {
73   // ctor : title is used to identify the layout
74   //        GPS2 = 5 modules (EMC + PPSD)
75   //        IHEP = 5 modules (EMC + CPV )
76   // We use 2 arrays of hits :
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.
84
85   fPinElectronicNoise = 0.010 ;
86   fDigitThreshold      = 0.1 ;   // 1 GeV 
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
91   // and the TreeD at the end of the event (branch is set in FinishEvent() ).
92   
93   fTmpHits= new TClonesArray("AliPHOSHit",1000) ;
94
95   fNTmpHits = fNhits = 0 ;
96
97   fDigits = new TClonesArray("AliPHOSDigit",1000) ;
98
99   fIshunt     =  1 ; // All hits are associated with primary particles
100   
101   // Create array of CPV modules for the IHEP's version of CPV
102
103   if ( strcmp(fGeom->GetName(),"IHEP") == 0 ) {
104     // Create array of CPVmodule of the size of PHOS modules number
105
106     if ( NULL==(fCPVModules=new TClonesArray("CPVModule",fGeom->GetNModules())) ) {
107       Error("AliPHOSv1","Can not create array of CPV modules");
108       exit(1);
109     }
110     TClonesArray &lcpvmodule = *fCPVModules;
111     for (Int_t i=0; i<fGeom->GetNModules(); i++) new(lcpvmodule[i]) CPVModule();
112   }
113   else {
114     // Create an empty array of CPVModule to satisfy
115     // AliPHOSv1::Streamer when writing root file
116
117     fCPVModules=new TClonesArray("CPVModule",0);
118
119   }
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() ) 
153     cout << "AliPHOS" << Version() << " : PHOS geometry intialized for " << fGeom->GetName() << endl ;
154   else
155     cout << "AliPHOS" << Version() << " : PHOS geometry initialization failed !" << endl ;   
156
157   // Defining the PHOS Reconstructioner
158  
159  fReconstructioner = Reconstructioner ;
160
161 }
162
163 //____________________________________________________________________________
164 AliPHOSv1::~AliPHOSv1()
165 {
166   // dtor
167
168   if ( fTmpHits) {
169     fTmpHits->Delete() ; 
170     delete fTmpHits ;
171     fTmpHits = 0 ; 
172   }
173
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 //   }
185   
186 //   if ( fTrackSegments ) {
187 //     fTrackSegments->Delete() ; 
188 //     delete fTrackSegments ;
189 //     fTrackSegments = 0 ; 
190 //   }
191
192 }
193
194 //____________________________________________________________________________
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.
198   // A PHOS hit is the sum of all hits in a single crystal
199   //   or in a single PPSD gas cell
200
201   Int_t hitCounter ;
202   TClonesArray &ltmphits = *fTmpHits ;
203   AliPHOSHit *newHit ;
204   AliPHOSHit *curHit ;
205   Bool_t deja = kFALSE ;
206
207   // In any case, fills the fTmpHit TClonesArray (with "accumulated hits")
208
209   newHit = new AliPHOSHit(shunt, primary, tracknumber, Id, hits) ;
210
211   // We do not want to save in TreeH the raw hits 
212   //  TClonesArray &lhits = *fHits;
213
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
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
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 }
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   
261   
262   Int_t i ;
263   Int_t relid[4];
264   Int_t j ; 
265   TClonesArray &lDigits = *fDigits ;
266   AliPHOSHit   * hit ;
267   AliPHOSDigit * newdigit ;
268   AliPHOSDigit * curdigit ;
269   Bool_t deja = kFALSE ; 
270   
271   for ( i = 0 ; i < fNTmpHits ; i++ ) {
272     hit = (AliPHOSHit*)fTmpHits->At(i) ;
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
278       newdigit = new AliPHOSDigit( -1               , hit->GetId(), Digitize( hit->GetEnergy() ) ) ;
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) ;
300
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) ; 
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());
324   }
325
326 }
327
328 //___________________________________________________________________________
329 void AliPHOSv1::MakeBranch(Option_t* opt)
330 {  
331   // Create new branche in the current Root Tree in the digit Tree
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   }
349
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
357 }
358
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   
377   fEmcRecPoints->Delete() ; 
378
379   if ( fEmcRecPoints && gAlice->TreeR() ) {
380     sprintf(branchname,"%sEmcRP",GetName()) ;
381     
382     gAlice->TreeR()->Branch(branchname, "TObjArray", &fEmcRecPoints, fBufferSize, splitlevel) ; 
383   }
384
385   fPpsdRecPoints->Delete() ; 
386
387
388   if ( fPpsdRecPoints && gAlice->TreeR() ) {
389     sprintf(branchname,"%sPpsdRP",GetName()) ;
390      
391     gAlice->TreeR()->Branch(branchname, "TObjArray", &fPpsdRecPoints, fBufferSize, splitlevel) ;
392   }
393
394   fTrackSegments->Delete() ; 
395
396   if ( fTrackSegments && gAlice->TreeR() ) { 
397     sprintf(branchname,"%sTS",GetName()) ;
398     gAlice->TreeR()->Branch(branchname, &fTrackSegments, fBufferSize) ;
399   }
400
401   fRecParticles->Delete() ; 
402
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     }
408   }
409   
410   // 3.
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);
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() ;
433   // 5.
434
435   gAlice->TreeR()->Write(0,TObject::kOverwrite) ;
436  
437   // Deleting reconstructed objects
438   ResetReconstruction();
439   
440 }
441
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 ) {
450     for (Int_t i=0; i<fGeom->GetNModules(); i++) ((CPVModule*)(*fCPVModules)[i]) -> Clear();
451   }
452 }  
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 }
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) ;
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
495 }
496
497 //____________________________________________________________________________
498
499 void AliPHOSv1::StepManager(void)
500 {
501   // Accumulates hits as long as the track stays in a single crystal or PPSD gas Cell
502
503   Int_t          relid[4] ;      // (box, layer, row, column) indices
504   Int_t          absid    ;      // absolute cell ID number
505   Float_t        xyze[4]  ;      // position wrt MRS and energy deposited
506   TLorentzVector pos      ;      // Lorentz vector of the track current position
507   Int_t          copy     ;
508
509   Int_t tracknumber =  gAlice->CurrentTrack() ; 
510   Int_t primary     =  gAlice->GetPrimary( gAlice->CurrentTrack() ); 
511   TString name      =  fGeom->GetName() ; 
512
513   if ( name == "GPS2" ) {                                       // ======> CPV is a GPS' PPSD
514
515     if( gMC->CurrentVolID(copy) == gMC->VolId("GCEL") ) // We are inside a gas cell 
516     {
517       gMC->TrackPosition(pos) ;
518       xyze[0] = pos[0] ;
519       xyze[1] = pos[1] ;
520       xyze[2] = pos[2] ;
521       xyze[3] = gMC->Edep() ; 
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
533         fGeom->RelToAbsNumbering(relid, absid) ; 
534
535         // add current hit to the hit list      
536         AddHit(fIshunt, primary, tracknumber, absid, xyze);
537
538       } // there is deposited energy 
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       
595       TClonesArray *cpvDigits = new TClonesArray("CPVDigit",0);   // array of digits for current hit
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++) {
607         CPVDigit  *cpvDigit1 = (CPVDigit*) cpvDigits->UncheckedAt(idigit);
608         Float_t x1 = cpvDigit1->GetXpad() ;
609         Float_t z1 = cpvDigit1->GetYpad() ;
610         for (Int_t jdigit=idigit+1; jdigit<ndigits; jdigit++) {
611           CPVDigit  *cpvDigit2 = (CPVDigit*) cpvDigits->UncheckedAt(jdigit);
612           Float_t x2 = cpvDigit2->GetXpad() ;
613           Float_t z2 = cpvDigit2->GetYpad() ;
614           if (x1==x2 && z1==z2) {
615             Float_t qsum = cpvDigit1->GetQpad() + cpvDigit2->GetQpad() ;
616             cpvDigit2->SetQpad(qsum) ;
617             cpvDigits->RemoveAt(idigit) ;
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++) {
627         CPVDigit  *cpvDigit = (CPVDigit*) cpvDigits->UncheckedAt(idigit);
628         relid[0] = ModuleNumber + 1 ;                             // CPV (or PHOS) module number
629         relid[1] =-1 ;                                            // means CPV
630         relid[2] = cpvDigit->GetXpad() ;                          // column number of a pad
631         relid[3] = cpvDigit->GetYpad() ;                          // row    number of a pad
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
653   
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
706   Float_t hitX  = zxhit[0];
707   Float_t hitZ  =-zxhit[1];
708   Float_t pX    = p.Px();
709   Float_t pZ    =-p.Pz();
710   Float_t pNorm = p.Py();
711   Float_t E     = dEdx;
712
713   Float_t dZY   = pZ/pNorm * fGeom->GetCPVGasThickness();
714   Float_t dXY   = pX/pNorm * fGeom->GetCPVGasThickness();
715   gRandom->Rannor(rnor1,rnor2);
716   E    *= (1 + detR*rnor1) *
717           TMath::Sqrt((1 + ( pow(dZY,2) + pow(dXY,2) ) / pow(fGeom->GetCPVGasThickness(),2)));
718   Float_t zhit1 = hitZ + fGeom->GetCPVActiveSize(1)/2 - dZY/2;
719   Float_t xhit1 = hitX + fGeom->GetCPVActiveSize(0)/2 - dXY/2;
720   Float_t zhit2 = zhit1 + dZY;
721   Float_t xhit2 = xhit1 + dXY;
722
723   Int_t   iwht1 = (Int_t) (xhit1 / celWr);           // wire (x) coordinate "in"
724   Int_t   iwht2 = (Int_t) (xhit2 / celWr);           // wire (x) coordinate "out"
725
726   Int_t   nIter;
727   Float_t zxe[3][5];
728   if (iwht1==iwht2) {                      // incline 1-wire hit
729     nIter = 2;
730     zxe[0][0] = (zhit1 + zhit2 - dZY*0.57735) / 2;
731     zxe[1][0] = (iwht1 + 0.5) * celWr;
732     zxe[2][0] =  E/2;
733     zxe[0][1] = (zhit1 + zhit2 + dZY*0.57735) / 2;
734     zxe[1][1] = (iwht1 + 0.5) * celWr;
735     zxe[2][1] =  E/2;
736   }
737   else if (TMath::Abs(iwht1-iwht2) != 1) { // incline 3-wire hit
738     nIter = 3;
739     Int_t iwht3 = (iwht1 + iwht2) / 2;
740     Float_t xwht1 = (iwht1 + 0.5) * celWr; // wire 1
741     Float_t xwht2 = (iwht2 + 0.5) * celWr; // wire 2
742     Float_t xwht3 = (iwht3 + 0.5) * celWr; // wire 3
743     Float_t xwr13 = (xwht1 + xwht3) / 2;   // center 13
744     Float_t xwr23 = (xwht2 + xwht3) / 2;   // center 23
745     Float_t dxw1  = xhit1 - xwr13;
746     Float_t dxw2  = xhit2 - xwr23;
747     Float_t egm1  = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
748     Float_t egm2  = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
749     Float_t egm3  =            celWr / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) + celWr );
750     zxe[0][0] = (dXY*(xwr13-xwht1)/dXY + zhit1 + zhit1) / 2;
751     zxe[1][0] =  xwht1;
752     zxe[2][0] =  E * egm1;
753     zxe[0][1] = (dXY*(xwr23-xwht1)/dXY + zhit1 + zhit2) / 2;
754     zxe[1][1] =  xwht2;
755     zxe[2][1] =  E * egm2;
756     zxe[0][2] =  dXY*(xwht3-xwht1)/dXY + zhit1;
757     zxe[1][2] =  xwht3;
758     zxe[2][2] =  E * egm3;
759   }
760   else {                                   // incline 2-wire hit
761     nIter = 2;
762     Float_t xwht1 = (iwht1 + 0.5) * celWr;
763     Float_t xwht2 = (iwht2 + 0.5) * celWr;
764     Float_t xwr12 = (xwht1 + xwht2) / 2;
765     Float_t dxw1  = xhit1 - xwr12;
766     Float_t dxw2  = xhit2 - xwr12;
767     Float_t egm1  = TMath::Abs(dxw1) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
768     Float_t egm2  = TMath::Abs(dxw2) / ( TMath::Abs(dxw1) + TMath::Abs(dxw2) );
769     zxe[0][0] = (zhit1 + zhit2 - dZY*egm1) / 2;
770     zxe[1][0] =  xwht1;
771     zxe[2][0] =  E * egm1;
772     zxe[0][1] = (zhit1 + zhit2 + dZY*egm2) / 2;
773     zxe[1][1] =  xwht2;
774     zxe[2][1] =  E * egm2;
775   }
776
777   // Finite size of ionization region
778
779   Int_t nCellZ  = fGeom->GetNumberOfPadsZ();
780   Int_t nCellX  = fGeom->GetNumberOfPadsPhi();
781   Int_t nz3     = (ngamz+1)/2;
782   Int_t nx3     = (ngamx+1)/2;
783   cpvDigits->Expand(nIter*ngamx*ngamz);
784   TClonesArray &ldigits = *(TClonesArray *)cpvDigits;
785
786   for (Int_t iter=0; iter<nIter; iter++) {
787
788     Float_t zhit = zxe[0][iter];
789     Float_t xhit = zxe[1][iter];
790     Float_t qhit = zxe[2][iter];
791     Float_t zcell = zhit / fGeom->GetPadSizeZ();
792     Float_t xcell = xhit / fGeom->GetPadSizePhi();
793     if ( zcell<=0      || xcell<=0 ||
794          zcell>=nCellZ || xcell>=nCellX) return;
795     Int_t izcell = (Int_t) zcell;
796     Int_t ixcell = (Int_t) xcell;
797     Float_t zc = zcell - izcell - 0.5;
798     Float_t xc = xcell - ixcell - 0.5;
799     for (Int_t iz=1; iz<=ngamz; iz++) {
800       Int_t kzg = izcell + iz - nz3;
801       if (kzg<=0 || kzg>nCellZ) continue;
802       Float_t zg = (Float_t)(iz-nz3) - zc;
803       for (Int_t ix=1; ix<=ngamx; ix++) {
804         Int_t kxg = ixcell + ix - nx3;
805         if (kxg<=0 || kxg>nCellX) continue;
806         Float_t xg = (Float_t)(ix-nx3) - xc;
807         
808         // Now calculate pad response
809         Float_t qpad = CPVPadResponseFunction(qhit,zg,xg);
810         qpad += qNoise*rnor2;
811         if (qpad<0) continue;
812         
813         // Fill the array with pad response ID and amplitude
814         new(ldigits[cpvDigits->GetEntriesFast()]) CPVDigit(kxg,kzg,qpad);
815 //      printf("(%2d,%2d,%5.3f) ",kxg,kzg,qpad);
816       }
817 //        cout << endl;
818     }
819 //      cout << endl;
820   }
821 }
822
823 //____________________________________________________________________________
824 Float_t AliPHOSv1::CPVPadResponseFunction(Float_t qhit, Float_t zhit, Float_t xhit) {
825   // ------------------------------------------------------------------------
826   // Calculate the amplitude in one CPV pad using the
827   // cumulative pad response function
828   // Author: Yuri Kharlov (after Serguei Sadovski)
829   // 3 October 2000
830   // ------------------------------------------------------------------------
831
832   Double_t dz = fGeom->GetPadSizeZ()   / 2;
833   Double_t dx = fGeom->GetPadSizePhi() / 2;
834   Double_t z  = zhit * fGeom->GetPadSizeZ();
835   Double_t x  = xhit * fGeom->GetPadSizePhi();
836   Double_t amplitude = qhit *
837     (CPVCumulPadResponse(z+dz,x+dx) - CPVCumulPadResponse(z+dz,x-dx) -
838      CPVCumulPadResponse(z-dz,x+dx) + CPVCumulPadResponse(z-dz,x-dx));
839   return (Float_t)amplitude;
840 }
841
842 //____________________________________________________________________________
843 Double_t AliPHOSv1::CPVCumulPadResponse(Double_t x, Double_t y) {
844   // ------------------------------------------------------------------------
845   // Cumulative pad response function
846   // It includes several terms from the CF decomposition in electrostatics
847   // Note: this cumulative function is wrong since omits some terms
848   //       but the cell amplitude obtained with it is correct because
849   //       these omitting terms cancel
850   // Author: Yuri Kharlov (after Serguei Sadovski)
851   // 3 October 2000
852   // ------------------------------------------------------------------------
853
854   const Double_t a=1.0;
855   const Double_t b=0.7;
856
857   Double_t r2       = x*x + y*y;
858   Double_t xy       = x*y;
859   Double_t cumulPRF = 0;
860   for (Int_t i=0; i<=4; i++) {
861     Double_t b1 = (2*i + 1) * b;
862     cumulPRF += TMath::Power(-1,i) * TMath::ATan( xy / (b1*TMath::Sqrt(b1*b1 + r2)) );
863   }
864   cumulPRF *= a/(2*TMath::Pi());
865   return cumulPRF;
866 }