Initialisation of some pointers
[u/mrichter/AliRoot.git] / MUON / AliMUON.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 $Log$
17 Revision 1.42  2001/01/17 20:53:40  hristov
18 Destructors corrected to avoid memory leaks
19
20 Revision 1.41  2000/12/21 22:12:40  morsch
21 Clean-up of coding rule violations,
22
23 Revision 1.40  2000/11/29 20:32:26  gosset
24 Digitize:
25 1. correction for array index out of bounds
26 2. one printout commented
27
28 Revision 1.39  2000/11/12 17:17:03  pcrochet
29 BuildGeometry of AliMUON for trigger chambers delegated to AliMUONSegmentationTriggerX (same strategy as for tracking chambers)
30
31 Revision 1.38  2000/11/06 09:20:43  morsch
32 AliMUON delegates part of BuildGeometry() to AliMUONSegmentation using the
33 Draw() method. This avoids code and parameter replication.
34
35 Revision 1.37  2000/10/26 09:53:37  pcrochet
36 put back trigger chambers in the display (there was a problem in buildgeometry)
37
38 Revision 1.36  2000/10/25 19:51:18  morsch
39 Correct x-position of chambers.
40
41 Revision 1.35  2000/10/24 19:46:21  morsch
42 BuildGeometry updated for slats in station 3-4.
43
44 Revision 1.34  2000/10/18 11:42:06  morsch
45 - AliMUONRawCluster contains z-position.
46 - Some clean-up of useless print statements during initialisations.
47
48 Revision 1.33  2000/10/09 14:01:57  morsch
49 Unused variables removed.
50
51 Revision 1.32  2000/10/06 09:08:10  morsch
52 Built geometry includes slat geometry for event display.
53
54 Revision 1.31  2000/10/02 21:28:08  fca
55 Removal of useless dependecies via forward declarations
56
57 Revision 1.30  2000/10/02 16:58:29  egangler
58 Cleaning of the code :
59 -> coding conventions
60 -> void Streamers
61 -> some useless includes removed or replaced by "class" statement
62
63 Revision 1.29  2000/07/28 13:49:38  morsch
64 SetAcceptance defines inner and outer chamber radii according to angular acceptance.
65 Can be used for simple acceptance studies.
66
67 Revision 1.28  2000/07/22 16:43:15  morsch
68 Same comment as before, but now done correctly I hope (sorry it's Saturday evening)
69
70 Revision 1.27  2000/07/22 16:36:50  morsch
71 Change order of indices in creation (new) of xhit and yhit
72
73 Revision 1.26  2000/07/03 11:54:57  morsch
74 AliMUONSegmentation and AliMUONHitMap have been replaced by AliSegmentation and AliHitMap in STEER
75 The methods GetPadIxy and GetPadXxy of AliMUONSegmentation have changed name to GetPadI and GetPadC.
76
77 Revision 1.25  2000/06/29 12:34:09  morsch
78 AliMUONSegmentation class has been made independent of AliMUONChamber. This makes
79 it usable with any other geometry class. The link to the object to which it belongs is
80 established via an index. This assumes that there exists a global geometry manager
81 from which the pointer to the parent object can be obtained (in our case gAlice).
82
83 Revision 1.24  2000/06/28 15:16:35  morsch
84 (1) Client code adapted to new method signatures in AliMUONSegmentation (see comments there)
85 to allow development of slat-muon chamber simulation and reconstruction code in the MUON
86 framework. The changes should have no side effects (mostly dummy arguments).
87 (2) Hit disintegration uses 3-dim hit coordinates to allow simulation
88 of chambers with overlapping modules (MakePadHits, Disintegration).
89
90 Revision 1.23  2000/06/28 12:19:17  morsch
91 More consequent seperation of global input data services (AliMUONClusterInput singleton) and the
92 cluster and hit reconstruction algorithms in AliMUONClusterFinderVS.
93 AliMUONClusterFinderVS becomes the base class for clustering and hit reconstruction.
94 It requires two cathode planes. Small modifications in the code will make it usable for
95 one cathode plane and, hence, more general (for test beam data).
96 AliMUONClusterFinder is now obsolete.
97
98 Revision 1.22  2000/06/28 08:06:10  morsch
99 Avoid global variables in AliMUONClusterFinderVS by seperating the input data for the fit from the
100 algorithmic part of the class. Input data resides inside the AliMUONClusterInput singleton.
101 It also naturally takes care of the TMinuit instance.
102
103 Revision 1.21  2000/06/27 08:54:41  morsch
104 Problems with on constant array sizes (in hitMap, nmuon, xhit, yhit) corrected.
105
106 Revision 1.20  2000/06/26 14:02:38  morsch
107 Add class AliMUONConstants with MUON specific constants using static memeber data and access methods.
108
109 Revision 1.19  2000/06/22 13:40:51  morsch
110 scope problem on HP, "i" declared once
111 pow changed to TMath::Power (PH, AM)
112
113 Revision 1.18  2000/06/15 07:58:48  morsch
114 Code from MUON-dev joined
115
116 Revision 1.14.4.17  2000/06/14 14:36:46  morsch
117 - add TriggerCircuit (PC)
118 - add GlobalTrigger and LocalTrigger and specific methods (PC)
119
120 Revision 1.14.4.16  2000/06/09 21:20:28  morsch
121 Most coding rule violations corrected
122
123 Revision 1.14.4.15  2000/05/02 09:54:32  morsch
124 RULE RN17 violations corrected
125
126 Revision 1.14.4.12  2000/04/26 12:25:02  morsch
127 Code revised by P. Crochet:
128 - Z position of TriggerChamber changed according to A.Tournaire Priv.Comm.
129 - ToF included in the method MakePadHits
130 - inner radius of flange between beam shielding and trigger corrected
131 - Trigger global volume updated (according to the new geometry)
132
133 Revision 1.14.4.11  2000/04/19 19:42:08  morsch
134 Some changes of variable names curing viols and methods concerning
135 correlated clusters removed.
136
137 Revision 1.14.4.10  2000/03/22 16:44:07  gosset
138 Memory leak suppressed in function Digitise:
139 p_adr->Delete() instead of Clear (I.Chevrot and A.Baldisseri)
140
141 Revision 1.14.4.9  2000/03/20 18:15:25  morsch
142 Positions of trigger chambers corrected (P.C.)
143
144 Revision 1.14.4.8  2000/02/21 15:38:01  morsch
145 Call to AddHitList introduced to make this version compatible with head.
146
147 Revision 1.14.4.7  2000/02/20 07:45:53  morsch
148 Bugs in Trigger part of BuildGeomemetry corrected (P.C)
149
150 Revision 1.14.4.6  2000/02/17 14:28:54  morsch
151 Trigger included into initialization and digitization
152
153 Revision 1.14.4.5  2000/02/15 10:02:58  morsch
154 Log messages of previous revisions added
155
156 Revision 1.14.4.2  2000/02/04 10:57:34  gosset
157 Z position of the chambers:
158 it was the Z position of the stations;
159 it is now really the Z position of the chambers.
160    !!!! WARNING: THE CALLS TO "AliMUONChamber::SetZPOS"
161    !!!!                   AND "AliMUONChamber::ZPosition"
162    !!!! HAVE TO BE CHANGED TO "AliMUONChamber::"SetZ"
163    !!!!                   AND "AliMUONChamber::Z"
164
165 Revision 1.14.4.3  2000/02/04 16:19:04  gosset
166 Correction for mis-spelling of NCH 
167
168 Revision 1.14.4.4  2000/02/15 09:43:38  morsch
169 Log message added
170
171 */
172
173
174 ///////////////////////////////////////////////
175 //  Manager and hits classes for set:MUON     //
176 ////////////////////////////////////////////////
177
178 #include <TTUBE.h>
179 #include <TBRIK.h>
180 #include <TRotMatrix.h>
181 #include <TGeometry.h>
182 #include <TNode.h> 
183 #include <TTree.h> 
184 #include <TRandom.h> 
185 #include <TObject.h>
186 #include <TVector.h>
187 #include <TObjArray.h>
188 #include <TMinuit.h>
189 #include <TParticle.h>
190 #include <TROOT.h>
191 #include <TFile.h>
192 #include <TNtuple.h>
193 #include <TCanvas.h>
194 #include <TPad.h>
195 #include <TDirectory.h>
196 #include <TObjectTable.h>
197 #include <AliPDG.h>
198 #include <TTUBE.h>
199
200 #include "AliMUON.h"
201 #include "AliMUONHit.h"
202 #include "AliMUONPadHit.h"
203 #include "AliMUONDigit.h"
204 #include "AliMUONTransientDigit.h"
205 #include "AliMUONRawCluster.h"
206 #include "AliMUONLocalTrigger.h"
207 #include "AliMUONGlobalTrigger.h"
208 #include "AliMUONTriggerCircuit.h"
209 #include "AliHitMap.h"
210 #include "AliMUONHitMapA1.h"
211 #include "AliMUONChamberTrigger.h"
212 #include "AliMUONConstants.h"
213 #include "AliMUONClusterFinderVS.h"
214 #include "AliMUONTriggerDecision.h"
215 #include "AliRun.h"
216 #include "AliMC.h"
217 #include "AliMUONClusterInput.h"
218 #include "iostream.h"
219 #include "AliCallf77.h" 
220 #include "AliConst.h" 
221
222 // Defaults parameters for Z positions of chambers
223 // taken from values for "stations" in AliMUON::AliMUON
224 //     const Float_t zch[7]={528, 690., 975., 1249., 1449., 1610, 1710.};
225 // and from array "dstation" in AliMUONv1::CreateGeometry
226 //          Float_t dstation[5]={20., 20., 20, 20., 20.};
227 //     for tracking chambers,
228 //          according to (Z1 = zch - dstation) and  (Z2 = zch + dstation)
229 //          for the first and second chambers in the station, respectively,
230 // and from "DTPLANES" in AliMUONv1::CreateGeometry
231 //           const Float_t DTPLANES = 15.;
232 //     for trigger chambers,
233 //          according to (Z1 = zch) and  (Z2 = zch + DTPLANES)
234 //          for the first and second chambers in the station, respectively
235
236 ClassImp(AliMUON)
237 //___________________________________________
238 AliMUON::AliMUON()
239 {
240 // Default Constructor
241 //
242    fIshunt       = 0;
243    fHits         = 0;
244    fPadHits      = 0;
245    fNPadHits     = 0;
246    fDchambers    = 0;
247    fTriggerCircuits = 0;     // cp new design of AliMUONTriggerDecision
248    fNdch         = 0;
249    fRawClusters  = 0;
250    fNrawch       = 0;
251    fGlobalTrigger   = 0;
252    fNLocalTrigger   = 0;
253    fLocalTrigger    = 0;
254    fNLocalTrigger   = 0;
255    fAccMin          = 0.;
256    fAccMax          = 0.;   
257    fAccCut          = kFALSE;
258    fChambers        = 0;
259    fHits2           = 0;
260    fPadHits2        = 0;
261    fFileName        = 0;
262    fTrH1            = 0;
263 }
264  
265 //___________________________________________
266 AliMUON::AliMUON(const char *name, const char *title)
267        : AliDetector(name,title)
268 {
269 //Begin_Html
270 /*
271 <img src="gif/alimuon.gif">
272 */
273 //End_Html
274    fHits2           = 0;
275    fPadHits2        = 0;
276    fFileName        = 0;
277    fTrH1            = 0;
278
279    fHits     = new TClonesArray("AliMUONHit",1000);
280    gAlice->AddHitList(fHits);
281    fPadHits = new TClonesArray("AliMUONPadHit",10000);
282    fNPadHits  =  0;
283    fIshunt     =  0;
284
285    fNdch      = new Int_t[AliMUONConstants::NCh()];
286
287    fDchambers = new TObjArray(AliMUONConstants::NCh());
288
289    Int_t i;
290    
291    for (i=0; i<AliMUONConstants::NCh() ;i++) {
292        (*fDchambers)[i] = new TClonesArray("AliMUONDigit",10000); 
293        fNdch[i]=0;
294    }
295
296    fNrawch      = new Int_t[AliMUONConstants::NTrackingCh()];
297
298    fRawClusters = new TObjArray(AliMUONConstants::NTrackingCh());
299
300    for (i=0; i<AliMUONConstants::NTrackingCh();i++) {
301        (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); 
302        fNrawch[i]=0;
303    }
304
305    fGlobalTrigger = new TClonesArray("AliMUONGlobalTrigger",1);    
306    fNGlobalTrigger = 0;
307    fLocalTrigger  = new TClonesArray("AliMUONLocalTrigger",234);    
308    fNLocalTrigger = 0;
309
310    SetMarkerColor(kRed);
311 //
312 //
313 //
314 //
315
316     Int_t ch;
317
318     fChambers = new TObjArray(AliMUONConstants::NCh());
319
320     // Loop over stations
321     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
322       // Loop over 2 chambers in the station
323         for (Int_t stCH = 0; stCH < 2; stCH++) {
324 //
325 //    
326 //    Default Parameters for Muon Tracking Stations
327
328
329             ch = 2 * st + stCH;
330 //
331             if (ch < AliMUONConstants::NTrackingCh()) {
332                 (*fChambers)[ch] = new AliMUONChamber(ch);
333             } else {
334                 (*fChambers)[ch] = new AliMUONChamberTrigger(ch);
335             }
336             
337             AliMUONChamber* chamber = (AliMUONChamber*) (*fChambers)[ch];
338             
339             chamber->SetGid(0);
340             // Default values for Z of chambers
341             chamber->SetZ(AliMUONConstants::DefaultChamberZ(ch));
342 //
343             chamber->InitGeo(AliMUONConstants::DefaultChamberZ(ch));
344 //          Set chamber inner and outer radius to default
345             chamber->SetRInner(AliMUONConstants::Dmin(st)/2);
346             chamber->SetROuter(AliMUONConstants::Dmax(st)/2);
347 //
348         } // Chamber stCH (0, 1) in 
349     }     // Station st (0...)
350     fMaxStepGas=0.01; 
351     fMaxStepAlu=0.1; 
352     fMaxDestepGas=-1;
353     fMaxDestepAlu=-1;
354 //
355    fMaxIterPad   = 0;
356    fCurIterPad   = 0;
357 //
358    fAccMin          = 0.;
359    fAccMax          = 0.;   
360    fAccCut          = kFALSE;
361
362    // cp new design of AliMUONTriggerDecision
363    fTriggerCircuits = new TObjArray(AliMUONConstants::NTriggerCircuit());
364    for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
365      (*fTriggerCircuits)[circ] = new AliMUONTriggerCircuit();     
366    }
367    // cp new design of AliMUONTriggerDecision
368
369 }
370  
371 //___________________________________________
372 AliMUON::AliMUON(const AliMUON& rMUON)
373 {
374 // Dummy copy constructor
375     ;
376     
377 }
378
379 AliMUON::~AliMUON()
380 {
381 // Destructor
382     printf("Calling AliMUON destructor !!!\n");
383     
384     Int_t i;
385     fIshunt  = 0;
386  
387     // Delete TObjArrays
388  
389     if (fChambers){
390       fChambers->Delete();
391       delete fChambers;
392     }
393  
394     if (fTriggerCircuits){
395       fTriggerCircuits->Delete();
396       delete fTriggerCircuits;
397     }
398  
399     if (fDchambers){
400       fDchambers->Delete();
401       delete fDchambers;
402     }
403  
404     if (fRawClusters){
405       fRawClusters->Delete();
406       delete fRawClusters;
407     }
408     for (i=0;i<AliMUONConstants::NTrackingCh();i++) {
409       fNrawch[i]=0;
410     }
411  
412     // Delete TClonesArrays
413  
414     if (fPadHits){
415       fPadHits->Delete();
416       delete fPadHits;
417     }
418  
419     if (fGlobalTrigger){
420       fGlobalTrigger->Delete();
421       delete fGlobalTrigger;
422     }
423     fNGlobalTrigger = 0;
424     
425     if (fLocalTrigger){
426       fLocalTrigger->Delete();
427       delete fLocalTrigger;
428     }
429     fNLocalTrigger = 0;
430
431     if (fHits2){
432       fHits2->Delete();
433       delete fHits2;
434     }
435
436     if (fPadHits2){
437       fPadHits2->Delete();
438       delete fPadHits2;
439     }
440
441     if (fHits) {
442       fHits->Delete();
443       delete fHits;
444     }
445
446     // Delete hits tree for background event
447
448     if (fTrH1) {
449       fTrH1->Delete();
450       delete fTrH1;
451     }
452 }
453  
454 //___________________________________________
455 void AliMUON::AddHit(Int_t track, Int_t *vol, Float_t *hits)
456 {
457   TClonesArray &lhits = *fHits;
458   new(lhits[fNhits++]) AliMUONHit(fIshunt,track,vol,hits);
459 }
460 //___________________________________________
461 void AliMUON::AddPadHit(Int_t *clhits)
462 {
463    TClonesArray &lclusters = *fPadHits;
464    new(lclusters[fNPadHits++]) AliMUONPadHit(clhits);
465 }
466 //_____________________________________________________________________________
467 void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
468 {
469     //
470     // Add a MUON digit to the list
471     //
472
473     TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
474     new(ldigits[fNdch[id]++]) AliMUONDigit(tracks,charges,digits);
475 }
476
477 //_____________________________________________________________________________
478 void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
479 {
480     //
481     // Add a MUON digit to the list
482     //
483
484     TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
485     new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
486 }
487
488 //___________________________________________
489 void AliMUON::AddGlobalTrigger(Int_t *singlePlus, Int_t *singleMinus,
490                                Int_t *singleUndef,
491                                Int_t *pairUnlike, Int_t *pairLike)
492 {
493 // add a MUON Global Trigger to the list (only one GlobalTrigger per event !)
494   TClonesArray &globalTrigger = *fGlobalTrigger;
495   new(globalTrigger[fNGlobalTrigger++]) 
496     AliMUONGlobalTrigger(singlePlus, singleMinus,  singleUndef, pairUnlike, 
497                          pairLike);
498 }
499 //___________________________________________
500 void AliMUON::AddLocalTrigger(Int_t *localtr)
501 {
502 // add a MUON Local Trigger to the list
503   TClonesArray &localTrigger = *fLocalTrigger;
504   new(localTrigger[fNLocalTrigger++]) AliMUONLocalTrigger(localtr);
505 }
506
507 //___________________________________________
508 void AliMUON::BuildGeometry()
509 {
510 // Geometry for event display
511   for (Int_t i=0; i<7; i++) {
512     for (Int_t j=0; j<2; j++) {
513       Int_t id=2*i+j+1;
514       this->Chamber(id-1).SegmentationModel(1)->Draw("eventdisplay");
515     }
516   }
517 }
518
519 //___________________________________________
520 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
521 {
522    return 9999;
523 }
524
525 //___________________________________________
526 void AliMUON::MakeBranch(Option_t* option)
527 {
528     // Create Tree branches for the MUON.
529     const Int_t kBufferSize = 4000;
530     char branchname[30];
531     sprintf(branchname,"%sCluster",GetName());
532     
533     AliDetector::MakeBranch(option);
534     
535     if (fPadHits   && gAlice->TreeH()) {
536         gAlice->TreeH()->Branch(branchname,&fPadHits, kBufferSize);
537         printf("Making Branch %s for clusters\n",branchname);
538     }
539     
540 // one branch for digits per chamber
541     Int_t i;
542     
543     for (i=0; i<AliMUONConstants::NCh() ;i++) {
544         sprintf(branchname,"%sDigits%d",GetName(),i+1);
545         
546         if (fDchambers   && gAlice->TreeD()) {
547             gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), kBufferSize);
548             printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
549         }       
550     }
551     
552     printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
553     
554 // one branch for raw clusters per chamber
555     for (i=0; i<AliMUONConstants::NTrackingCh() ;i++) {
556         sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
557         
558         if (fRawClusters   && gAlice->TreeR()) {
559             gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), kBufferSize);
560             printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
561         }       
562     }
563     
564 // one branch for global trigger
565     sprintf(branchname,"%sGlobalTrigger",GetName());
566     if (fGlobalTrigger && gAlice->TreeR()) {  
567         gAlice->TreeR()->Branch(branchname,&fGlobalTrigger,kBufferSize);
568         printf("Making Branch %s for Global Trigger\n",branchname);
569     }
570 // one branch for local trigger
571     sprintf(branchname,"%sLocalTrigger",GetName());
572     if (fLocalTrigger && gAlice->TreeR()) {  
573         gAlice->TreeR()->Branch(branchname,&fLocalTrigger,kBufferSize);
574         printf("Making Branch %s for Local Trigger\n",branchname);
575     }
576     
577 }
578
579 //___________________________________________
580 void AliMUON::SetTreeAddress()
581 {
582   // Set branch address for the Hits and Digits Tree.
583   char branchname[30];
584   AliDetector::SetTreeAddress();
585
586   TBranch *branch;
587   TTree *treeH = gAlice->TreeH();
588   TTree *treeD = gAlice->TreeD();
589   TTree *treeR = gAlice->TreeR();
590
591   if (treeH) {
592     if (fPadHits) {
593       branch = treeH->GetBranch("MUONCluster");
594       if (branch) branch->SetAddress(&fPadHits);
595     }
596   }
597
598   if (treeD) {
599       for (int i=0; i<AliMUONConstants::NCh(); i++) {
600           sprintf(branchname,"%sDigits%d",GetName(),i+1);
601           if (fDchambers) {
602               branch = treeD->GetBranch(branchname);
603               if (branch) branch->SetAddress(&((*fDchambers)[i]));
604           }
605       }
606   }
607
608   // printf("SetTreeAddress --- treeR address  %p \n",treeR);
609
610   if (treeR) {
611       for (int i=0; i<AliMUONConstants::NTrackingCh(); i++) {
612           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
613           if (fRawClusters) {
614               branch = treeR->GetBranch(branchname);
615               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
616           }
617       }
618
619       if (fLocalTrigger) {
620         branch = treeR->GetBranch("MUONLocalTrigger");
621         if (branch) branch->SetAddress(&fLocalTrigger);
622       }
623       if (fGlobalTrigger) {
624         branch = treeR->GetBranch("MUONGlobalTrigger");
625         if (branch) branch->SetAddress(&fGlobalTrigger);
626       }
627   }
628 }
629 //___________________________________________
630 void AliMUON::ResetHits()
631 {
632   // Reset number of clusters and the cluster array for this detector
633   AliDetector::ResetHits();
634   fNPadHits = 0;
635   if (fPadHits) fPadHits->Clear();
636 }
637
638 //____________________________________________
639 void AliMUON::ResetDigits()
640 {
641     //
642     // Reset number of digits and the digits array for this detector
643     //
644     for ( int i=0;i<AliMUONConstants::NCh();i++ ) {
645         if ((*fDchambers)[i])    ((TClonesArray*)(*fDchambers)[i])->Clear();
646         if (fNdch)  fNdch[i]=0;
647     }
648 }
649 //____________________________________________
650 void AliMUON::ResetRawClusters()
651 {
652     //
653     // Reset number of raw clusters and the raw clust array for this detector
654     //
655     for ( int i=0;i<AliMUONConstants::NTrackingCh();i++ ) {
656         if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
657         if (fNrawch)  fNrawch[i]=0;
658     }
659 }
660
661 //____________________________________________
662 void AliMUON::ResetTrigger()
663 {
664   //  Reset Local and Global Trigger 
665   fNGlobalTrigger = 0;
666   if (fGlobalTrigger) fGlobalTrigger->Clear();
667   fNLocalTrigger = 0;
668   if (fLocalTrigger) fLocalTrigger->Clear();
669 }
670
671 //____________________________________________
672 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
673 {
674 // Set the pad size for chamber id and cathode isec
675     Int_t i=2*(id-1);
676     ((AliMUONChamber*) (*fChambers)[i])  ->SetPadSize(isec,p1,p2);
677     ((AliMUONChamber*) (*fChambers)[i+1])->SetPadSize(isec,p1,p2);
678 }
679
680 //___________________________________________
681 void AliMUON::SetChambersZ(const Float_t *Z)
682 {
683   // Set Z values for all chambers (tracking and trigger)
684   // from the array pointed to by "Z"
685     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
686         ((AliMUONChamber*) ((*fChambers)[ch]))->SetZ(Z[ch]);
687     return;
688 }
689
690 //___________________________________________
691 void AliMUON::SetChambersZToDefault()
692 {
693   // Set Z values for all chambers (tracking and trigger)
694   // to default values
695   SetChambersZ(AliMUONConstants::DefaultChamberZ());
696   return;
697 }
698
699 //___________________________________________
700 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
701 {
702 // Set the inverse charge slope for chamber id
703     Int_t i=2*(id-1);
704     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSlope(p1);
705     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
706 }
707
708 //___________________________________________
709 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
710 {
711 // Set sigma of charge spread for chamber id
712     Int_t i=2*(id-1);
713     ((AliMUONChamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
714     ((AliMUONChamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
715 }
716
717 //___________________________________________
718 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
719 {
720 // Set integration limits for charge spread
721     Int_t i=2*(id-1);
722     ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
723     ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
724 }
725
726 //___________________________________________
727 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
728 {
729 // Set maximum number for ADCcounts (saturation)
730     Int_t i=2*(id-1);
731     ((AliMUONChamber*) (*fChambers)[i])->SetMaxAdc(p1);
732     ((AliMUONChamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
733 }
734
735 //___________________________________________
736 void AliMUON::SetMaxStepGas(Float_t p1)
737 {
738 // Set stepsize in gas
739      fMaxStepGas=p1;
740 }
741
742 //___________________________________________
743 void AliMUON::SetMaxStepAlu(Float_t p1)
744 {
745 // Set step size in Alu
746     fMaxStepAlu=p1;
747 }
748
749 //___________________________________________
750 void AliMUON::SetMaxDestepGas(Float_t p1)
751 {
752 // Set maximum step size in Gas
753     fMaxDestepGas=p1;
754 }
755
756 //___________________________________________
757 void AliMUON::SetMaxDestepAlu(Float_t p1)
758 {
759 // Set maximum step size in Alu
760     fMaxDestepAlu=p1;
761 }
762 //___________________________________________
763 void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
764 {
765 // Set acceptance cuts 
766    fAccCut=acc;
767    fAccMin=angmin*TMath::Pi()/180;
768    fAccMax=angmax*TMath::Pi()/180;
769    Int_t ch;
770    if (acc) {
771        for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
772            // Loop over 2 chambers in the station
773            for (Int_t stCH = 0; stCH < 2; stCH++) {
774                ch = 2 * st + stCH;
775 //         Set chamber inner and outer radius according to acceptance cuts
776                Chamber(ch).SetRInner(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMin));
777                Chamber(ch).SetROuter(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMax));
778            } // chamber loop
779        } // station loop
780    }
781 }
782 //___________________________________________
783 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
784 {
785 // Set the segmentation for chamber id cathode isec
786     ((AliMUONChamber*) (*fChambers)[id])->SetSegmentationModel(isec, segmentation);
787
788 }
789 //___________________________________________
790 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
791 {
792 // Set the response for chamber id
793     ((AliMUONChamber*) (*fChambers)[id])->SetResponseModel(response);
794 }
795
796 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
797 {
798 // Set ClusterFinder for chamber id
799     ((AliMUONChamber*) (*fChambers)[id])->SetReconstructionModel(reconst);
800 }
801
802 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
803 {
804 // Set number of segmented cathods for chamber id
805     ((AliMUONChamber*) (*fChambers)[id])->SetNsec(nsec);
806 }
807
808
809 //___________________________________________
810
811
812
813 void AliMUON::MakePadHits(Float_t xhit,Float_t yhit, Float_t zhit,
814                           Float_t eloss, Float_t tof,  Int_t idvol)
815 {
816 //
817 //  Calls the charge disintegration method of the current chamber and adds
818 //  the simulated cluster to the root treee 
819 //
820     Int_t clhits[7];
821     Float_t newclust[6][500];
822     Int_t nnew;
823     
824     
825 //
826 //  Integrated pulse height on chamber
827
828     
829     clhits[0]=fNhits+1;
830 //
831 //
832 //    if (idvol == 6) printf("\n ->Disintegration %f %f %f", xhit, yhit, eloss );
833     
834
835     ((AliMUONChamber*) (*fChambers)[idvol])
836         ->DisIntegration(eloss, tof, xhit, yhit, zhit, nnew, newclust);
837     Int_t ic=0;
838 //    if (idvol == 6) printf("\n nnew  %d \n", nnew);
839 //
840 //  Add new clusters
841     for (Int_t i=0; i<nnew; i++) {
842         if (Int_t(newclust[3][i]) > 0) {
843             ic++;
844 // Cathode plane
845             clhits[1] = Int_t(newclust[5][i]);
846 //  Cluster Charge
847             clhits[2] = Int_t(newclust[0][i]);
848 //  Pad: ix
849             clhits[3] = Int_t(newclust[1][i]);
850 //  Pad: iy 
851             clhits[4] = Int_t(newclust[2][i]);
852 //  Pad: charge
853             clhits[5] = Int_t(newclust[3][i]);
854 //  Pad: chamber sector
855             clhits[6] = Int_t(newclust[4][i]);
856             
857             AddPadHit(clhits);
858         }
859     }
860 }
861
862 //----------------------------------------------------------------------
863
864 void AliMUON::Digitise(Int_t nev,Int_t bgrEvent,Option_t *option,Option_t *opt,Text_t *filename)
865 {
866     // keep galice.root for signal and name differently the file for 
867     // background when add! otherwise the track info for signal will be lost !
868   
869     static Bool_t first=kTRUE;
870     static TFile *file;
871     char *addBackground = strstr(option,"Add");
872
873
874     AliMUONChamber*   iChamber;
875     AliSegmentation*  segmentation;
876
877     
878     Int_t trk[50];
879     Int_t chtrk[50];  
880     TObjArray *list=new TObjArray;
881     static TClonesArray *pAddress=0;
882     if(!pAddress) pAddress=new TClonesArray("TVector",1000);
883     Int_t digits[5]; 
884
885     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
886     AliHitMap** hitMap= new AliHitMap* [AliMUONConstants::NCh()];
887     for (Int_t i=0; i<AliMUONConstants::NCh(); i++) {hitMap[i]=0;}
888     if (addBackground ) {
889         if(first) {
890             fFileName=filename;
891             cout<<"filename"<<fFileName<<endl;
892             file=new TFile(fFileName);
893             cout<<"I have opened "<<fFileName<<" file "<<endl;
894             fHits2     = new TClonesArray("AliMUONHit",1000  );
895             fPadHits2 = new TClonesArray("AliMUONPadHit",10000);
896         }           
897         first=kFALSE;
898         file->cd();
899         //file->ls();
900         // Get Hits Tree header from file
901         if(fHits2) fHits2->Clear();
902         if(fPadHits2) fPadHits2->Clear();
903         if(fTrH1) delete fTrH1;
904         fTrH1=0;
905         
906         char treeName[20];
907         sprintf(treeName,"TreeH%d",bgrEvent);
908         fTrH1 = (TTree*)gDirectory->Get(treeName);
909         //printf("fTrH1 %p of treename %s for event %d \n",fTrH1,treeName,bgrEvent);
910         
911         if (!fTrH1) {
912             printf("ERROR: cannot find Hits Tree for event:%d\n",bgrEvent);
913         }
914         // Set branch addresses
915         TBranch *branch;
916         char branchname[20];
917         sprintf(branchname,"%s",GetName());
918         if (fTrH1 && fHits2) {
919             branch = fTrH1->GetBranch(branchname);
920             if (branch) branch->SetAddress(&fHits2);
921         }
922         if (fTrH1 && fPadHits2) {
923             branch = fTrH1->GetBranch("MUONCluster");
924             if (branch) branch->SetAddress(&fPadHits2);
925         }
926 // test
927         //Int_t ntracks1 =(Int_t)fTrH1->GetEntries();
928         //printf("background - ntracks1 - %d\n",ntracks1);
929     }
930     //
931     // loop over cathodes
932     //
933     AliHitMap* hm;
934     Int_t countadr=0;
935     for (int icat=0; icat<2; icat++) { 
936         Int_t counter=0;
937         Int_t * nmuon = new Int_t [AliMUONConstants::NCh()];
938         for (Int_t i =0; i<AliMUONConstants::NCh(); i++) {
939             iChamber=(AliMUONChamber*) (*fChambers)[i];
940             if (iChamber->Nsec()==1 && icat==1) {
941                 continue;
942             } else {
943                 segmentation=iChamber->SegmentationModel(icat+1);
944             }
945             hitMap[i] = new AliMUONHitMapA1(segmentation, list);
946             nmuon[i]=0;
947         }
948         //printf("Start loop over tracks \n");     
949 //
950 //   Loop over tracks
951 //
952
953         TTree *treeH = gAlice->TreeH();
954         Int_t ntracks =(Int_t) treeH->GetEntries();
955         Int_t jj;
956         Int_t nmaxmuon = 5;
957
958         Float_t ** xhit = new Float_t * [AliMUONConstants::NCh()];
959         for (jj=0; jj<AliMUONConstants::NCh(); jj++)
960           xhit[jj] = new Float_t[nmaxmuon];
961         Float_t ** yhit = new Float_t * [AliMUONConstants::NCh()];
962         for (jj=0; jj<AliMUONConstants::NCh(); jj++)
963           yhit[jj] = new Float_t[nmaxmuon];
964
965         for (Int_t track=0; track<ntracks; track++) {
966             gAlice->ResetHits();
967             treeH->GetEvent(track);
968 //
969 //   Loop over hits
970             for(AliMUONHit* mHit=(AliMUONHit*)pMUON->FirstHit(-1); 
971                 mHit;
972                 mHit=(AliMUONHit*)pMUON->NextHit()) 
973             {
974                 Int_t   nch   = mHit->fChamber-1;  // chamber number
975                 if (nch > AliMUONConstants::NCh()-1) continue;
976 //              if (nch > 9) continue;
977                 iChamber = &(pMUON->Chamber(nch));
978                 // new 17.07.99
979                 if (addBackground) {
980
981                     if (mHit->fParticle == kMuonPlus 
982                         || mHit->fParticle == kMuonMinus) {
983                       // mark muon hits
984                       if (nmuon[nch] < nmaxmuon) {
985                         xhit[nch][nmuon[nch]]=mHit->X();
986                         yhit[nch][nmuon[nch]]=mHit->Y();
987                         nmuon[nch]++;
988                       }
989                       // ignore muon if too many compared to nmaxmuon
990                       else printf("AliMUON::Digitize: nmuon %d ==> ignored\n",nmuon[nch]);
991                     }
992                 }
993                 
994
995
996                 
997 //
998 // Loop over pad hits
999                 for (AliMUONPadHit* mPad=
1000                          (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits);
1001                      mPad;
1002                      mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits))
1003                 {
1004                     Int_t cathode  = mPad->fCathode;    // cathode number
1005                     Int_t ipx      = mPad->fPadX;       // pad number on X
1006                     Int_t ipy      = mPad->fPadY;       // pad number on Y
1007                     Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
1008 //                  printf("\n Pad: %d %d %d %d", ipx, ipy, cathode,nch);
1009 //
1010 //
1011                     
1012                     if (cathode != (icat+1)) continue;
1013                     // fill the info array
1014 //                  Float_t thex, they, thez;
1015                     segmentation=iChamber->SegmentationModel(cathode);
1016 //                  segmentation->GetPadC(ipx,ipy,thex,they,thez);
1017 //                  Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1018 //                  if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1019
1020                     new((*pAddress)[countadr++]) TVector(2);
1021                     TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
1022                     trinfo(0)=(Float_t)track;
1023                     trinfo(1)=(Float_t)iqpad;
1024
1025                     digits[0]=ipx;
1026                     digits[1]=ipy;
1027                     digits[2]=iqpad;
1028                     digits[3]=iqpad;
1029                     if (mHit->fParticle == kMuonPlus ||
1030                         mHit->fParticle == kMuonMinus) {
1031                         digits[4]=mPad->fHitNumber;
1032                     } else digits[4]=-1;
1033
1034                     AliMUONTransientDigit* pdigit;
1035                     // build the list of fired pads and update the info
1036                     if (!hitMap[nch]->TestHit(ipx, ipy)) {
1037
1038                         list->AddAtAndExpand(
1039                             new AliMUONTransientDigit(nch,digits),counter);
1040                         
1041                         hitMap[nch]->SetHit(ipx, ipy, counter);
1042                         counter++;
1043                         pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
1044                         // list of tracks
1045                         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
1046                         trlist->Add(&trinfo);
1047                     } else {
1048                         pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
1049                         // update charge
1050                         (*pdigit).fSignal+=iqpad;
1051                         (*pdigit).fPhysics+=iqpad;                      
1052                         // update list of tracks
1053                         TObjArray* trlist=(TObjArray*)pdigit->TrackList();
1054                         Int_t lastEntry=trlist->GetLast();
1055                         TVector *pTrack=(TVector*)trlist->At(lastEntry);
1056                         TVector &ptrk=*pTrack;
1057                         Int_t lastTrack  = Int_t(ptrk(0));
1058                         Int_t lastCharge = Int_t(ptrk(1));
1059                         if (lastTrack==track) {
1060                             lastCharge+=iqpad;
1061                             trlist->RemoveAt(lastEntry);
1062                             trinfo(0)=lastTrack;
1063                             trinfo(1)=lastCharge;
1064                             trlist->AddAt(&trinfo,lastEntry);
1065                         } else {
1066                             trlist->Add(&trinfo);
1067                         }
1068                         // check the track list
1069                         Int_t nptracks=trlist->GetEntriesFast();
1070                         if (nptracks > 2) {
1071                             for (Int_t tr=0;tr<nptracks;tr++) {
1072                                 TVector *ppTrack=(TVector*)trlist->At(tr);
1073                                 TVector &pptrk=*ppTrack;
1074                                 trk[tr]   = Int_t(pptrk(0));
1075                                 chtrk[tr] = Int_t(pptrk(1));
1076                             }
1077                         } // end if nptracks
1078                     } //  end if pdigit
1079                 } //end loop over clusters
1080             } // hit loop
1081         } // track loop
1082
1083         // open the file with background
1084        
1085         if (addBackground) {
1086             ntracks =(Int_t)fTrH1->GetEntries();
1087 //
1088 //   Loop over tracks
1089 //
1090             for (Int_t track=0; track<ntracks; track++) {
1091
1092                 if (fHits2)       fHits2->Clear();
1093                 if (fPadHits2)   fPadHits2->Clear();
1094
1095                 fTrH1->GetEvent(track);
1096 //
1097 //   Loop over hits
1098                 AliMUONHit* mHit;
1099                 for(int i=0;i<fHits2->GetEntriesFast();++i) 
1100                 {       
1101                     mHit=(AliMUONHit*) (*fHits2)[i];
1102                     Int_t   nch   = mHit->fChamber-1;  // chamber number
1103                     if (nch >9) continue;
1104                     iChamber = &(pMUON->Chamber(nch));
1105 //                  Int_t rmin = (Int_t)iChamber->RInner();
1106 //                  Int_t rmax = (Int_t)iChamber->ROuter();
1107                     Float_t xbgr=mHit->X();
1108                     Float_t ybgr=mHit->Y();
1109                     Bool_t cond=kFALSE;
1110                     
1111                     for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
1112                         Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
1113                             +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
1114                         if (dist<100) cond=kTRUE;
1115                     }
1116                     if (!cond) continue;
1117                     
1118 //
1119 // Loop over pad hits
1120                     for (AliMUONPadHit* mPad=
1121                              (AliMUONPadHit*)pMUON->FirstPad(mHit,fPadHits2);
1122                          mPad;
1123                          mPad=(AliMUONPadHit*)pMUON->NextPad(fPadHits2))
1124                     {
1125                         //                  mPad = (AliMUONPadHit*) (*fPadHits2)[j];
1126                         Int_t cathode  = mPad->fCathode;    // cathode number
1127                         Int_t ipx      = mPad->fPadX;       // pad number on X
1128                         Int_t ipy      = mPad->fPadY;       // pad number on Y
1129                         Int_t iqpad    = Int_t(mPad->fQpad);// charge per pad
1130
1131                         if (cathode != (icat+1)) continue;
1132 //                      printf("\n Pad: %d %d %d", ipx, ipy, cathode);
1133                         
1134 //                      Float_t thex, they, thez;
1135 //                      segmentation=iChamber->SegmentationModel(cathode);
1136 //                      segmentation->GetPadC(ipx,ipy,thex,they,thez);
1137 //                      Float_t rpad=TMath::Sqrt(thex*thex+they*they);
1138 //                      if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
1139                         new((*pAddress)[countadr++]) TVector(2);
1140                         TVector &trinfo=*((TVector*) (*pAddress)[countadr-1]);
1141                         trinfo(0)=-1;  // tag background
1142                         trinfo(1)=-1;
1143                         
1144                         digits[0]=ipx;
1145                         digits[1]=ipy;
1146                         digits[2]=iqpad;
1147                         digits[3]=0;
1148                         digits[4]=-1;
1149                         
1150                         AliMUONTransientDigit* pdigit;
1151                         // build the list of fired pads and update the info
1152                         if (!hitMap[nch]->TestHit(ipx, ipy)) {
1153                             list->AddAtAndExpand(new AliMUONTransientDigit(nch,digits),counter);
1154                             
1155                             hitMap[nch]->SetHit(ipx, ipy, counter);
1156                             counter++;
1157                             
1158                             pdigit=(AliMUONTransientDigit*)list->At(list->GetLast());
1159                             // list of tracks
1160                             TObjArray *trlist=(TObjArray*)pdigit->
1161                                 TrackList();
1162                             trlist->Add(&trinfo);
1163                         } else {
1164                             pdigit=(AliMUONTransientDigit*) hitMap[nch]->GetHit(ipx, ipy);
1165                             // update charge
1166                             (*pdigit).fSignal+=iqpad;
1167                             
1168                             // update list of tracks
1169                             TObjArray* trlist=(TObjArray*)pdigit->
1170                                 TrackList();
1171                             Int_t lastEntry=trlist->GetLast();
1172                             TVector *pTrack=(TVector*)trlist->
1173                                 At(lastEntry);
1174                             TVector &ptrk=*pTrack;
1175                             Int_t lastTrack=Int_t(ptrk(0));
1176                             if (lastTrack==-1) {
1177                                 continue;
1178                             } else {
1179                                 trlist->Add(&trinfo);
1180                             }
1181                                 // check the track list
1182                             Int_t nptracks=trlist->GetEntriesFast();
1183                             if (nptracks > 0) {
1184                                 for (Int_t tr=0;tr<nptracks;tr++) {
1185                                     TVector *ppTrack=(TVector*)trlist->At(tr);
1186                                     TVector &pptrk=*ppTrack;
1187                                     trk[tr]=Int_t(pptrk(0));
1188                                     chtrk[tr]=Int_t(pptrk(1));
1189                                 }
1190                             } // end if nptracks
1191                         } //  end if pdigit
1192                     } //end loop over clusters
1193                 } // hit loop
1194             } // track loop
1195             //Int_t nentr2=list->GetEntriesFast();
1196             //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1197             TTree *fAli=gAlice->TreeK();
1198             TFile *file=NULL;
1199             
1200             if (fAli) file =fAli->GetCurrentFile();
1201             file->cd();
1202         } // if addBackground
1203         delete [] xhit;
1204         delete [] yhit;
1205
1206         Int_t tracks[10];
1207         Int_t charges[10];
1208         Int_t nentries=list->GetEntriesFast();
1209         
1210         for (Int_t nent=0;nent<nentries;nent++) {
1211             AliMUONTransientDigit *address=(AliMUONTransientDigit*)list->At(nent);
1212             if (address==0) continue; 
1213             Int_t ich=address->fChamber;
1214             Int_t q=address->fSignal; 
1215             iChamber=(AliMUONChamber*) (*fChambers)[ich];
1216 //
1217 //  Digit Response (noise, threshold, saturation, ...)
1218 //              if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise; 
1219             AliMUONResponse * response=iChamber->ResponseModel();
1220             q=response->DigitResponse(q);
1221             
1222             if (!q) continue;
1223             
1224             digits[0]=address->fPadX;
1225             digits[1]=address->fPadY;
1226             digits[2]=q;
1227             digits[3]=address->fPhysics;
1228             digits[4]=address->fHit;
1229             
1230             TObjArray* trlist=(TObjArray*)address->TrackList();
1231             Int_t nptracks=trlist->GetEntriesFast();
1232             //printf("nptracks, trlist   %d  %p\n",nptracks,trlist);
1233
1234             // this was changed to accomodate the real number of tracks
1235             if (nptracks > 10) {
1236                 cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1237                 nptracks=10;
1238             }
1239             if (nptracks > 2) {
1240                 printf("Attention - nptracks > 2  %d \n",nptracks);
1241                 printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1242             }
1243             for (Int_t tr=0;tr<nptracks;tr++) {
1244                 TVector *ppP=(TVector*)trlist->At(tr);
1245                 if(!ppP ) printf("ppP - %p\n",ppP);
1246                 TVector &pp  =*ppP;
1247                 tracks[tr]=Int_t(pp(0));
1248                 charges[tr]=Int_t(pp(1));
1249                 //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
1250             }      //end loop over list of tracks for one pad
1251             // Sort list of tracks according to charge
1252             if (nptracks > 1) {
1253                 SortTracks(tracks,charges,nptracks);
1254             }
1255             if (nptracks < 10 ) {
1256                 for (Int_t i=nptracks; i<10; i++) {
1257                     tracks[i]=0;
1258                     charges[i]=0;
1259                 }
1260             }
1261             
1262             // fill digits
1263             pMUON->AddDigits(ich,tracks,charges,digits);
1264             // delete trlist;
1265         }
1266         //cout<<"I'm out of the loops for digitisation"<<endl;
1267         //      gAlice->GetEvent(nev);
1268         gAlice->TreeD()->Fill();
1269         pMUON->ResetDigits();
1270         list->Delete();
1271
1272         
1273         for(Int_t ii=0;ii<AliMUONConstants::NCh();++ii) {
1274             if (hitMap[ii]) {
1275                 hm=hitMap[ii];
1276                 delete hm;
1277                 hitMap[ii]=0;
1278             }
1279         }
1280         delete [] nmuon;    
1281     } //end loop over cathodes
1282     delete [] hitMap;
1283     char hname[30];
1284     sprintf(hname,"TreeD%d",nev);
1285     gAlice->TreeD()->Write(hname);
1286     // reset tree
1287     gAlice->TreeD()->Reset();
1288     delete list;
1289     
1290     pAddress->Delete();
1291     // gObjectTable->Print();
1292 }
1293
1294 void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
1295 {
1296   //
1297   // Sort the list of tracks contributing to a given digit
1298   // Only the 3 most significant tracks are acctually sorted
1299   //
1300   
1301   //
1302   //  Loop over signals, only 3 times
1303   //
1304   
1305   Int_t qmax;
1306   Int_t jmax;
1307   Int_t idx[3] = {-2,-2,-2};
1308   Int_t jch[3] = {-2,-2,-2};
1309   Int_t jtr[3] = {-2,-2,-2};
1310   Int_t i,j,imax;
1311   
1312   if (ntr<3) imax=ntr;
1313   else imax=3;
1314   for(i=0;i<imax;i++){
1315     qmax=0;
1316     jmax=0;
1317     
1318     for(j=0;j<ntr;j++){
1319       
1320       if((i == 1 && j == idx[i-1]) 
1321          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
1322       
1323       if(charges[j] > qmax) {
1324         qmax = charges[j];
1325         jmax=j;
1326       }       
1327     } 
1328     
1329     if(qmax > 0) {
1330       idx[i]=jmax;
1331       jch[i]=charges[jmax]; 
1332       jtr[i]=tracks[jmax]; 
1333     }
1334     
1335   } 
1336   
1337   for(i=0;i<3;i++){
1338     if (jtr[i] == -2) {
1339          charges[i]=0;
1340          tracks[i]=0;
1341     } else {
1342          charges[i]=jch[i];
1343          tracks[i]=jtr[i];
1344     }
1345   }
1346
1347 }
1348
1349 //___________________________________________
1350 void AliMUON::Trigger(Int_t nev){
1351 // call the Trigger Algorithm and fill TreeR
1352
1353   Int_t singlePlus[3]  = {0,0,0}; 
1354   Int_t singleMinus[3] = {0,0,0}; 
1355   Int_t singleUndef[3] = {0,0,0};
1356   Int_t pairUnlike[3]  = {0,0,0}; 
1357   Int_t pairLike[3]    = {0,0,0};
1358
1359   ResetTrigger();
1360   AliMUONTriggerDecision* decision= new AliMUONTriggerDecision(1);
1361   decision->Trigger();   
1362   decision->GetGlobalTrigger(singlePlus, singleMinus, singleUndef,
1363                              pairUnlike, pairLike);
1364 // add a local trigger in the list 
1365   AddGlobalTrigger(singlePlus, singleMinus, singleUndef, pairUnlike, pairLike);
1366   Int_t i;
1367   
1368   for (Int_t icirc=0; icirc<AliMUONConstants::NTriggerCircuit(); icirc++) { 
1369       if(decision->GetITrigger(icirc)==1) {
1370           Int_t localtr[7]={0,0,0,0,0,0,0};      
1371           Int_t loLpt[2]={0,0}; Int_t loHpt[2]={0,0}; Int_t loApt[2]={0,0};
1372           decision->GetLutOutput(icirc, loLpt, loHpt, loApt);
1373           localtr[0] = icirc;
1374           localtr[1] = decision->GetStripX11(icirc);
1375           localtr[2] = decision->GetDev(icirc);
1376           localtr[3] = decision->GetStripY11(icirc);
1377           for (i=0; i<2; i++) {    // convert the Lut output in 1 digit 
1378               localtr[4] = localtr[4]+Int_t(loLpt[i]*TMath::Power(2,i));
1379               localtr[5] = localtr[5]+Int_t(loHpt[i]*TMath::Power(2,i));
1380               localtr[6] = localtr[6]+Int_t(loApt[i]*TMath::Power(2,i));
1381           }
1382           AddLocalTrigger(localtr);  // add a local trigger in the list
1383       }
1384   }
1385   delete decision;
1386
1387   gAlice->TreeR()->Fill();
1388   ResetTrigger();
1389   char hname[30];
1390   sprintf(hname,"TreeR%d",nev);
1391   gAlice->TreeR()->Write(hname);
1392   gAlice->TreeR()->Reset();
1393   printf("\n End of trigger for event %d", nev);
1394 }
1395
1396
1397 //____________________________________________
1398 void AliMUON::FindClusters(Int_t nev,Int_t lastEntry)
1399 {
1400 //
1401 //  Perform cluster finding
1402 //
1403     TClonesArray *dig1, *dig2;
1404     Int_t ndig, k;
1405     dig1 = new TClonesArray("AliMUONDigit",1000);
1406     dig2 = new TClonesArray("AliMUONDigit",1000);
1407     AliMUONDigit *digit;
1408 //
1409 // Loop on chambers and on cathode planes
1410 //
1411     
1412     for (Int_t ich=0;ich<10;ich++) {
1413         AliMUONChamber* iChamber=(AliMUONChamber*) (*fChambers)[ich];
1414         AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();    
1415         gAlice->ResetDigits();
1416         gAlice->TreeD()->GetEvent(lastEntry);
1417         TClonesArray *muonDigits  = this->DigitsAddress(ich);
1418         ndig=muonDigits->GetEntriesFast();
1419         printf("\n 1 Found %d digits in %p %d", ndig, muonDigits,ich);
1420         TClonesArray &lhits1 = *dig1;
1421         Int_t n=0;
1422         for (k=0; k<ndig; k++) {
1423             digit=(AliMUONDigit*) muonDigits->UncheckedAt(k);
1424             if (rec->TestTrack(digit->fTracks[0]))
1425                 new(lhits1[n++]) AliMUONDigit(*digit);
1426         }
1427         gAlice->ResetDigits();
1428         gAlice->TreeD()->GetEvent(lastEntry+1);
1429         muonDigits  = this->DigitsAddress(ich);
1430         ndig=muonDigits->GetEntriesFast();
1431         printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
1432         TClonesArray &lhits2 = *dig2;
1433         n=0;
1434         
1435         for (k=0; k<ndig; k++) {
1436             digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
1437             if (rec->TestTrack(digit->fTracks[0]))
1438             new(lhits2[n++]) AliMUONDigit(*digit);
1439         }
1440
1441         if (rec) {       
1442             AliMUONClusterInput::Instance()->SetDigits(ich, dig1, dig2);
1443             rec->FindRawClusters();
1444         }
1445         dig1->Delete();
1446         dig2->Delete();
1447     } // for ich
1448     gAlice->TreeR()->Fill();
1449     ResetRawClusters();
1450     char hname[30];
1451     sprintf(hname,"TreeR%d",nev);
1452     gAlice->TreeR()->Write(hname);
1453     gAlice->TreeR()->Reset();
1454     printf("\n End of cluster finding for event %d", nev);
1455     
1456     delete dig1;
1457     delete dig2;
1458     //gObjectTable->Print();
1459 }
1460  
1461
1462 void AliMUON::Streamer(TBuffer &R__b)
1463 {
1464    // Stream an object of class AliMUON.
1465       AliMUONChamber        *iChamber;
1466       AliMUONTriggerCircuit *iTriggerCircuit;
1467       AliSegmentation       *segmentation;
1468       AliMUONResponse       *response;
1469       TClonesArray          *digitsaddress;
1470       TClonesArray          *rawcladdress;
1471       Int_t i;
1472       if (R__b.IsReading()) {
1473           Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1474           AliDetector::Streamer(R__b);
1475           R__b >> fNPadHits;
1476           R__b >> fPadHits; // diff
1477           R__b >> fNLocalTrigger;       
1478           R__b >> fLocalTrigger;       
1479           R__b >> fNGlobalTrigger;       
1480           R__b >> fGlobalTrigger;   
1481           R__b >> fDchambers;
1482           R__b >> fRawClusters;
1483           R__b.ReadArray(fNdch);
1484           R__b.ReadArray(fNrawch);
1485           R__b >> fAccCut;
1486           R__b >> fAccMin;
1487           R__b >> fAccMax; 
1488           R__b >> fChambers;
1489           R__b >> fTriggerCircuits;
1490           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
1491               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
1492               iTriggerCircuit->Streamer(R__b);
1493           }
1494 // Stream chamber related information
1495           for (i =0; i<AliMUONConstants::NCh(); i++) {
1496               iChamber=(AliMUONChamber*) (*fChambers)[i];
1497               iChamber->Streamer(R__b);
1498               if (iChamber->Nsec()==1) {
1499                   segmentation=iChamber->SegmentationModel(1);
1500                   if (segmentation)
1501                       segmentation->Streamer(R__b);
1502               } else {
1503                   segmentation=iChamber->SegmentationModel(1);
1504                   if (segmentation)
1505                       segmentation->Streamer(R__b);
1506                   if (segmentation)
1507                       segmentation=iChamber->SegmentationModel(2);
1508                   segmentation->Streamer(R__b);
1509               }
1510               response=iChamber->ResponseModel();
1511               if (response)
1512                   response->Streamer(R__b);       
1513               digitsaddress=(TClonesArray*) (*fDchambers)[i];
1514               digitsaddress->Streamer(R__b);
1515               if (i < AliMUONConstants::NTrackingCh()) {
1516                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1517                   rawcladdress->Streamer(R__b);
1518               }
1519           }
1520           
1521       } else {
1522           R__b.WriteVersion(AliMUON::IsA());
1523           AliDetector::Streamer(R__b);
1524           R__b << fNPadHits;
1525           R__b << fPadHits; // diff
1526           R__b << fNLocalTrigger;       
1527           R__b << fLocalTrigger;       
1528           R__b << fNGlobalTrigger;       
1529           R__b << fGlobalTrigger; 
1530           R__b << fDchambers;
1531           R__b << fRawClusters;
1532           R__b.WriteArray(fNdch, AliMUONConstants::NCh());
1533           R__b.WriteArray(fNrawch, AliMUONConstants::NTrackingCh());
1534           
1535           R__b << fAccCut;
1536           R__b << fAccMin;
1537           R__b << fAccMax; 
1538           
1539           R__b << fChambers;
1540           R__b << fTriggerCircuits;
1541           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
1542               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
1543               iTriggerCircuit->Streamer(R__b);
1544           }
1545           for (i =0; i<AliMUONConstants::NCh(); i++) {
1546               iChamber=(AliMUONChamber*) (*fChambers)[i];
1547               iChamber->Streamer(R__b);
1548               if (iChamber->Nsec()==1) {
1549                   segmentation=iChamber->SegmentationModel(1);
1550                   if (segmentation)
1551                       segmentation->Streamer(R__b);
1552               } else {
1553                   segmentation=iChamber->SegmentationModel(1);
1554                   if (segmentation)
1555                       segmentation->Streamer(R__b);
1556                   segmentation=iChamber->SegmentationModel(2);
1557                   if (segmentation)
1558                       segmentation->Streamer(R__b);
1559               }
1560               response=iChamber->ResponseModel();
1561               if (response)
1562                   response->Streamer(R__b);
1563               digitsaddress=(TClonesArray*) (*fDchambers)[i];
1564               digitsaddress->Streamer(R__b);
1565               if (i < AliMUONConstants::NTrackingCh()) {
1566                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1567                   rawcladdress->Streamer(R__b);
1568               }
1569           }
1570       }
1571 }
1572 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters) 
1573 {
1574 //
1575     // Initialise the pad iterator
1576     // Return the address of the first padhit for hit
1577     TClonesArray *theClusters = clusters;
1578     Int_t nclust = theClusters->GetEntriesFast();
1579     if (nclust && hit->fPHlast > 0) {
1580         AliMUON::fMaxIterPad=hit->fPHlast;
1581         AliMUON::fCurIterPad=hit->fPHfirst;
1582         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
1583     } else {
1584         return 0;
1585     }
1586 }
1587
1588 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters) 
1589 {
1590 // Get next pad (in iterator) 
1591 //
1592     AliMUON::fCurIterPad++;
1593     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
1594         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
1595     } else {
1596         return 0;
1597     }
1598 }
1599
1600
1601 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
1602 {
1603 //
1604 //  Return rawcluster (icluster) for chamber ichamber and cathode icathod
1605 //  Obsolete ??
1606     TClonesArray *muonRawCluster  = RawClustAddress(ichamber);
1607     ResetRawClusters();
1608     TTree *treeR = gAlice->TreeR();
1609     Int_t nent=(Int_t)treeR->GetEntries();
1610     treeR->GetEvent(nent-2+icathod-1);
1611     //treeR->GetEvent(icathod);
1612     //Int_t nrawcl = (Int_t)muonRawCluster->GetEntriesFast();
1613
1614     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)muonRawCluster->UncheckedAt(icluster);
1615     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
1616     
1617     return  mRaw;
1618 }
1619
1620 AliMUON& AliMUON::operator = (const AliMUON& rhs)
1621 {
1622 // copy operator
1623 // dummy version
1624     return *this;
1625 }
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644