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