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