Introduce TLorentzVector in AliMC::GetSecondary. Thanks to I.Hrivnacova
[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 /*
17 $Log$
18 Revision 1.12  1999/10/07 21:08:10  fca
19 Corrections by G.Chabratova
20
21 Revision 1.11  1999/10/05 17:15:45  fca
22 Minor syntax for the Alpha OSF
23
24 Revision 1.10  1999/10/01 09:24:40  fca
25 Protect against no current file in FinishEvent
26
27 Revision 1.9  1999/09/29 09:24:20  fca
28 Introduction of the Copyright and cvs Log
29
30 */
31
32 ////////////////////////////////////////////////
33 //  Manager and hits classes for set:MUON     //
34 ////////////////////////////////////////////////
35
36 #include <TTUBE.h>
37 #include <TBRIK.h>
38 #include <TRotMatrix.h>
39 #include <TNode.h> 
40 #include <TTree.h> 
41 #include <TRandom.h> 
42 #include <TObject.h>
43 #include <TVector.h>
44 #include <TObjArray.h>
45 #include <TMinuit.h>
46 #include <TParticle.h>
47 #include <TROOT.h>
48 #include <TFile.h>
49 #include <TNtuple.h>
50 #include <TCanvas.h>
51 #include <TPad.h>
52 #include <TDirectory.h>
53 #include <TObjectTable.h>
54 #include <AliPDG.h>
55
56 #include "AliMUON.h"
57 #include "TTUBE.h"
58 #include "AliMUONClusterFinder.h"
59 #include "AliRun.h"
60 #include "AliMC.h"
61 #include "iostream.h"
62 #include "AliCallf77.h" 
63
64 #ifndef WIN32 
65 # define reco_init       reco_init_
66 # define cutpxz          cutpxz_
67 # define sigmacut        sigmacut_
68 # define xpreci          xpreci_
69 # define ypreci          ypreci_
70 # define reconstmuon     reconstmuon_
71 # define trackf_read_geant     trackf_read_geant_
72 # define trackf_read_spoint     trackf_read_spoint_
73 # define chfill          chfill_
74 # define chfill2         chfill2_
75 # define chf1            chf1_
76 # define chfnt           chfnt_
77 # define hist_create     hist_create_
78 # define hist_closed     hist_closed_
79 # define rndm            rndm_
80 # define fcn             fcn_
81 # define trackf_fit      trackf_fit_
82 # define prec_fit        prec_fit_
83 # define fcnfit          fcnfit_
84 # define reco_term       reco_term_
85 #else 
86 # define reco_init       RECO_INIT
87 # define cutpxz          CUTPXZ
88 # define sigmacut        SIGMACUT
89 # define xpreci          XPRECI
90 # define ypreci          YPRECI
91 # define reconstmuon     RECONSTMUON
92 # define trackf_read_geant     TRACKF_READ_GEANT
93 # define trackf_read_spoint     TRACKF_READ_SPOINT
94 # define chfill          CHFILL
95 # define chfill2         CHFILL2
96 # define chf1            CHF1
97 # define chfnt           CHFNT
98 # define hist_create     HIST_CREATE
99 # define hist_closed     HIST_CLOSED
100 # define rndm            RNDM
101 # define fcn             FCN
102 # define trackf_fit      TRACKF_FIT
103 # define prec_fit        PREC_FIT
104 # define fcnfit          FCNFIT
105 # define reco_term       RECO_TERM
106 #endif 
107
108 extern "C"
109 {
110 void type_of_call reco_init(Double_t &, Double_t &, Double_t &);
111 void type_of_call reco_term();
112 void type_of_call cutpxz(Double_t &);
113 void type_of_call sigmacut(Double_t &);
114 void type_of_call xpreci(Double_t &);
115 void type_of_call ypreci(Double_t &);
116 void type_of_call reconstmuon(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &);
117 void type_of_call trackf_read_geant(Int_t *, Double_t *, Double_t *, Double_t *, Int_t *, Int_t *, Double_t *, Double_t *, Double_t *, Double_t *,Int_t &, Double_t *, Double_t *, Double_t *, Int_t &, Int_t &, Double_t *, Double_t *, Double_t *, Double_t *); 
118 void type_of_call trackf_read_spoint(Int_t *, Double_t *, Double_t *, Double_t *, Int_t *, Int_t *, Double_t *, Double_t *, Double_t *, Double_t *,Int_t &, Double_t *, Double_t *, Double_t *, Int_t &, Int_t &, Double_t *, Double_t *, Double_t *, Double_t *); 
119 void type_of_call chfill(Int_t &, Float_t &, Float_t &, Float_t &);
120 void type_of_call chfill2(Int_t &, Float_t &, Float_t &, Float_t &);
121 void type_of_call chf1(Int_t &, Float_t &, Float_t &);
122 void type_of_call chfnt(Int_t &, Int_t &, Int_t *, Int_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *, Float_t *);
123 void type_of_call hist_create();
124 void type_of_call hist_closed();
125 void type_of_call fcnf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
126 void type_of_call fcn(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
127 void type_of_call trackf_fit(Int_t &, Double_t *, Double_t *, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
128 void type_of_call prec_fit(Double_t &, Double_t &, Double_t &, Double_t &, Double_t&, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
129 void type_of_call fcnfitf(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
130 void type_of_call fcnfit(Int_t &, Double_t *, Double_t &, Double_t *, Int_t &, Int_t &);
131 Float_t type_of_call rndm() {return gRandom->Rndm();}
132 }
133
134 void fcnfwrap(Int_t &i1, Double_t *d1, Double_t &d2,
135         Double_t *d3, Int_t i2)
136 {
137    fcnf(i1,d1,d2,d3,i2);
138 }
139
140 void fcnfitfwrap(Int_t &i1, Double_t *d1, Double_t &d2,
141         Double_t *d3, Int_t i2)
142 {
143    fcnfitf(i1,d1,d2,d3,i2);
144 }
145
146
147 // Static variables for the pad-hit iterator routines
148 static Int_t sMaxIterPad=0;
149 static Int_t sCurIterPad=0;
150 static TTree *TrH1;
151 static TTree *TK1;
152 static TClonesArray *fHits2;        //Listof hits for one track only
153 static TClonesArray *fClusters2;    //List of clusters for one track only
154 static TClonesArray *fParticles2;   //List of particles in the Kine tree
155 ClassImp(AliMUON)
156 //___________________________________________
157 AliMUON::AliMUON()
158 {
159    fIshunt     = 0;
160    fHits       = 0;
161    fClusters   = 0;
162    fNclusters  = 0;
163    fDchambers  = 0;
164    fNdch       = 0;
165    fRawClusters= 0;
166    fNrawch     = 0;
167    fCathCorrel= 0;
168    fNcorch     = 0;
169    fTreeC = 0;
170
171    // modifs perso
172    fSPxzCut    = 0;
173    fSSigmaCut  = 0;
174    fSXPrec     = 0; 
175    fSYPrec     = 0;
176 }
177  
178 //___________________________________________
179 AliMUON::AliMUON(const char *name, const char *title)
180        : AliDetector(name,title)
181 {
182 //Begin_Html
183 /*
184 <img src="gif/alimuon.gif">
185 */
186 //End_Html
187  
188    fHits     = new TClonesArray("AliMUONhit",1000);
189    fClusters = new TClonesArray("AliMUONcluster",10000);
190    fNclusters  =  0;
191    fIshunt     =  0;
192
193    fNdch      = new Int_t[10];
194
195    fDchambers = new TObjArray(10);
196
197    Int_t i;
198    
199    for (i=0; i<10 ;i++) {
200        (*fDchambers)[i] = new TClonesArray("AliMUONdigit",10000); 
201        fNdch[i]=0;
202    }
203
204    fNrawch      = new Int_t[10];
205
206    fRawClusters = new TObjArray(10);
207
208    for (i=0; i<10 ;i++) {
209        (*fRawClusters)[i] = new TClonesArray("AliMUONRawCluster",10000); 
210        fNrawch[i]=0;
211    }
212
213    fNcorch      = new Int_t[10];
214    fCathCorrel = new TObjArray(10);
215    for (i=0; i<10 ;i++) {
216        (*fCathCorrel)[i] = new TClonesArray("AliMUONcorrelation",1000); 
217        fNcorch[i]=0;
218    }
219
220    fTreeC = 0;
221
222 //   
223 // Transport angular cut
224    fAccCut=0;
225    fAccMin=2;
226    fAccMax=9;
227
228    // modifs perso
229    fSPxzCut   = 3.0;
230    fSSigmaCut = 1.0;
231    fSXPrec    = 0.01; 
232    fSYPrec    = 0.144;
233
234    SetMarkerColor(kRed);
235 }
236  
237 //___________________________________________
238 AliMUON::~AliMUON()
239 {
240
241     printf("Calling AliMUON destructor !!!\n");
242     
243   Int_t i;
244   fIshunt  = 0;
245   delete fHits;
246   delete fClusters;
247   delete fTreeC;
248
249   for (i=0;i<10;i++) {
250       delete (*fDchambers)[i];
251       fNdch[i]=0;
252   }
253   delete fDchambers;
254
255   for (i=0;i<10;i++) {
256       delete (*fRawClusters)[i];
257       fNrawch[i]=0;
258   }
259   delete fRawClusters;
260
261   for (i=0;i<10;i++) {
262       delete (*fCathCorrel)[i];
263       fNcorch[i]=0;
264   }
265   delete fCathCorrel;
266 }
267  
268 //___________________________________________
269 void AliMUON::AddHit(Int_t track, Int_t *vol, Float_t *hits)
270 {
271   TClonesArray &lhits = *fHits;
272   new(lhits[fNhits++]) AliMUONhit(fIshunt,track,vol,hits);
273 }
274 //___________________________________________
275 void AliMUON::AddCluster(Int_t *clhits)
276 {
277    TClonesArray &lclusters = *fClusters;
278    new(lclusters[fNclusters++]) AliMUONcluster(clhits);
279 }
280 //_____________________________________________________________________________
281 void AliMUON::AddDigits(Int_t id, Int_t *tracks, Int_t *charges, Int_t *digits)
282 {
283     //
284     // Add a MUON digit to the list
285     //
286
287     TClonesArray &ldigits = *((TClonesArray*)(*fDchambers)[id]);
288     new(ldigits[fNdch[id]++]) AliMUONdigit(tracks,charges,digits);
289 }
290
291 //_____________________________________________________________________________
292 void AliMUON::AddRawCluster(Int_t id, const AliMUONRawCluster& c)
293 {
294     //
295     // Add a MUON digit to the list
296     //
297
298     TClonesArray &lrawcl = *((TClonesArray*)(*fRawClusters)[id]);
299     new(lrawcl[fNrawch[id]++]) AliMUONRawCluster(c);
300 }
301 //_____________________________________________________________________________
302 void AliMUON::AddCathCorrel(Int_t id, Int_t *idx, Float_t *x, Float_t *y)
303 {
304     //
305     // Add a MUON digit to the list
306     //
307
308     TClonesArray &lcorrel = *((TClonesArray*)(*fCathCorrel)[id]);
309     new(lcorrel[fNcorch[id]++]) AliMUONcorrelation(idx,x,y);
310 }
311
312 //___________________________________________
313 void AliMUON::BuildGeometry()
314 {
315     TNode *Node, *NodeF, *Top;
316     const int kColorMUON = kBlue;
317     //
318     Top=gAlice->GetGeometry()->GetNode("alice");
319 // MUON
320 //
321 //  z-Positions of Chambers
322     const Float_t cz[5]={511., 686., 971., 1245., 1445.};
323 //
324 //  inner diameter
325     const Float_t dmi[5]={ 35.,  47.,  67.,   86.,  100.};
326 //
327 //  outer diameter
328     const Float_t dma[5]={183., 245., 346.,  520.,  520.};
329
330     TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90,  0, 90, 90, 0, 0);
331     TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
332     TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
333     TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90,  0, 0, 0);
334     
335
336     float rmin, rmax, dx, dy, dz, dr, zpos;
337     float dzc=4.;
338     char NameChamber[9], NameSense[9], NameFrame[9], NameNode[7];
339     for (Int_t i=0; i<5; i++) {
340         for (Int_t j=0; j<2; j++) {
341             Int_t id=2*i+j+1;
342             if (j==0) {
343                 zpos=cz[i]-dzc;
344             } else {
345                 zpos=cz[i]+dzc;
346             }
347             
348             
349             sprintf(NameChamber,"C_MUON%d",id);
350             sprintf(NameSense,"S_MUON%d",id);
351             sprintf(NameFrame,"F_MUON%d",id);   
352             rmin = dmi[i]/2.-3;
353             rmax = dma[i]/2.+3;
354             new TTUBE(NameChamber,"Mother","void",rmin,rmax,0.25,1.);
355             rmin = dmi[i]/2.;
356             rmax = dma[i]/2.;
357             new TTUBE(NameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
358             dx=(rmax-rmin)/2;
359             dy=3.;
360             dz=0.25;
361             TBRIK* FMUON = new TBRIK(NameFrame,"Frame","void",dx,dy,dz);
362             Top->cd();
363             sprintf(NameNode,"MUON%d",100+id);
364             Node = new TNode(NameNode,"ChamberNode",NameChamber,0,0,zpos,"");
365             Node->SetLineColor(kColorMUON);
366             fNodes->Add(Node);
367             Node->cd();
368             sprintf(NameNode,"MUON%d",200+id);
369             Node = new TNode(NameNode,"Sens. Region Node",NameSense,0,0,0,"");
370             Node->SetLineColor(kColorMUON);
371             fNodes->Add(Node);
372             Node->cd();
373             dr=dx+rmin;
374             sprintf(NameNode,"MUON%d",300+id);
375             NodeF = new TNode(NameNode,"Frame0",FMUON,dr, 0, 0,rot000,"");
376             NodeF->SetLineColor(kColorMUON);
377             fNodes->Add(NodeF);
378             Node->cd();
379             sprintf(NameNode,"MUON%d",400+id);
380             NodeF = new TNode(NameNode,"Frame1",FMUON,0 ,dr,0,rot090,"");
381             NodeF->SetLineColor(kColorMUON);
382             fNodes->Add(NodeF);
383             Node->cd();
384             sprintf(NameNode,"MUON%d",500+id);
385             NodeF = new TNode(NameNode,"Frame2",FMUON,-dr,0,0,rot180,"");
386             NodeF->SetLineColor(kColorMUON);
387             fNodes->Add(NodeF);
388             Node  ->cd();
389             sprintf(NameNode,"MUON%d",600+id);   
390             NodeF = new TNode(NameNode,"Frame3",FMUON,0,-dr,0,rot270,"");
391             NodeF->SetLineColor(kColorMUON);
392             fNodes->Add(NodeF);
393         }
394     }
395 }
396
397
398 //___________________________________________
399 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
400 {
401    return 9999;
402 }
403
404 //___________________________________________
405 void AliMUON::MakeBranch(Option_t* option)
406 {
407   // Create Tree branches for the MUON.
408   
409   const Int_t buffersize = 4000;
410   char branchname[30];
411   sprintf(branchname,"%sCluster",GetName());
412
413   AliDetector::MakeBranch(option);
414
415   if (fClusters   && gAlice->TreeH()) {
416     gAlice->TreeH()->Branch(branchname,&fClusters, buffersize);
417     printf("Making Branch %s for clusters\n",branchname);
418   }
419
420 // one branch for digits per chamber
421   Int_t i;
422   
423   for (i=0; i<10 ;i++) {
424       sprintf(branchname,"%sDigits%d",GetName(),i+1);
425       
426       if (fDchambers   && gAlice->TreeD()) {
427           gAlice->TreeD()->Branch(branchname,&((*fDchambers)[i]), buffersize);
428           printf("Making Branch %s for digits in chamber %d\n",branchname,i+1);
429       } 
430   }
431
432   //printf("Make Branch - TreeR address %p\n",gAlice->TreeR());
433
434 // one branch for raw clusters per chamber
435   for (i=0; i<10 ;i++) {
436       sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
437       
438       if (fRawClusters   && gAlice->TreeR()) {
439          gAlice->TreeR()->Branch(branchname,&((*fRawClusters)[i]), buffersize);
440          printf("Making Branch %s for raw clusters in chamber %d\n",branchname,i+1);
441       } 
442   }
443
444 }
445
446 //___________________________________________
447 void AliMUON::SetTreeAddress()
448 {
449   // Set branch address for the Hits and Digits Tree.
450   char branchname[30];
451   AliDetector::SetTreeAddress();
452
453   TBranch *branch;
454   TTree *treeH = gAlice->TreeH();
455   TTree *treeD = gAlice->TreeD();
456   TTree *treeR = gAlice->TreeR();
457
458   if (treeH) {
459     if (fClusters) {
460       branch = treeH->GetBranch("MUONCluster");
461       if (branch) branch->SetAddress(&fClusters);
462     }
463   }
464
465   if (treeD) {
466       for (int i=0; i<10; i++) {
467           sprintf(branchname,"%sDigits%d",GetName(),i+1);
468           if (fDchambers) {
469               branch = treeD->GetBranch(branchname);
470               if (branch) branch->SetAddress(&((*fDchambers)[i]));
471           }
472       }
473   }
474
475   // printf("SetTreeAddress --- treeR address  %p \n",treeR);
476
477   if (treeR) {
478       for (int i=0; i<10; i++) {
479           sprintf(branchname,"%sRawClusters%d",GetName(),i+1);
480           if (fRawClusters) {
481               branch = treeR->GetBranch(branchname);
482               if (branch) branch->SetAddress(&((*fRawClusters)[i]));
483           }
484       }
485   }
486
487 }
488 //___________________________________________
489 void AliMUON::ResetHits()
490 {
491   // Reset number of clusters and the cluster array for this detector
492   AliDetector::ResetHits();
493   fNclusters = 0;
494   if (fClusters) fClusters->Clear();
495 }
496
497 //____________________________________________
498 void AliMUON::ResetDigits()
499 {
500     //
501     // Reset number of digits and the digits array for this detector
502     //
503     for ( int i=0;i<10;i++ ) {
504         if ((*fDchambers)[i])    ((TClonesArray*)(*fDchambers)[i])->Clear();
505         if (fNdch)  fNdch[i]=0;
506     }
507 }
508 //____________________________________________
509 void AliMUON::ResetRawClusters()
510 {
511     //
512     // Reset number of raw clusters and the raw clust array for this detector
513     //
514     for ( int i=0;i<10;i++ ) {
515         if ((*fRawClusters)[i])    ((TClonesArray*)(*fRawClusters)[i])->Clear();
516         if (fNrawch)  fNrawch[i]=0;
517     }
518 }
519 //____________________________________________
520 void AliMUON::ResetCorrelation()
521 {
522     //
523     // Reset number of correl clusters and the correl clust array for 
524     // this detector
525     //
526     for ( int i=0;i<10;i++ ) {
527         if ((*fCathCorrel)[i])   ((TClonesArray*)(*fCathCorrel)[i])->Clear();
528         if (fNcorch)  fNcorch[i]=0;
529     }
530 }
531
532 //___________________________________________
533
534 void AliMUON::SetPADSIZ(Int_t id, Int_t isec, Float_t p1, Float_t p2)
535 {
536     Int_t i=2*(id-1);
537     ((AliMUONchamber*) (*fChambers)[i])  ->SetPADSIZ(isec,p1,p2);
538     ((AliMUONchamber*) (*fChambers)[i+1])->SetPADSIZ(isec,p1,p2);
539 }
540
541 //___________________________________________
542 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
543 {
544     Int_t i=2*(id-1);
545     ((AliMUONchamber*) (*fChambers)[i])->SetChargeSlope(p1);
546     ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSlope(p1);
547 }
548
549 //___________________________________________
550 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
551 {
552     Int_t i=2*(id-1);
553     ((AliMUONchamber*) (*fChambers)[i])->SetChargeSpread(p1,p2);
554     ((AliMUONchamber*) (*fChambers)[i+1])->SetChargeSpread(p1,p2);
555 }
556
557 //___________________________________________
558 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
559 {
560     Int_t i=2*(id-1);
561     ((AliMUONchamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
562     ((AliMUONchamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
563 }
564
565 //___________________________________________
566 void AliMUON::SetMaxAdc(Int_t id, Float_t p1)
567 {
568     Int_t i=2*(id-1);
569     ((AliMUONchamber*) (*fChambers)[i])->SetMaxAdc(p1);
570     ((AliMUONchamber*) (*fChambers)[i+1])->SetMaxAdc(p1);
571 }
572
573 //___________________________________________
574 void AliMUON::SetMaxStepGas(Float_t p1)
575 {
576      fMaxStepGas=p1;
577 }
578
579 //___________________________________________
580 void AliMUON::SetMaxStepAlu(Float_t p1)
581 {
582     fMaxStepAlu=p1;
583 }
584
585 //___________________________________________
586 void AliMUON::SetMaxDestepGas(Float_t p1)
587 {
588     fMaxDestepGas=p1;
589 }
590
591 //___________________________________________
592 void AliMUON::SetMaxDestepAlu(Float_t p1)
593 {
594     fMaxDestepAlu=p1;
595 }
596 //___________________________________________
597 void AliMUON::SetMuonAcc(Bool_t acc, Float_t angmin, Float_t angmax)
598 {
599    fAccCut=acc;
600    fAccMin=angmin;
601    fAccMax=angmax;
602 }
603 //___________________________________________
604 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliMUONsegmentation *segmentation)
605 {
606     ((AliMUONchamber*) (*fChambers)[id])->SegmentationModel(isec, segmentation);
607
608 }
609 //___________________________________________
610 void   AliMUON::SetResponseModel(Int_t id, AliMUONresponse *response)
611 {
612     ((AliMUONchamber*) (*fChambers)[id])->ResponseModel(response);
613 }
614
615 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinder *reconst)
616 {
617     ((AliMUONchamber*) (*fChambers)[id])->ReconstructionModel(reconst);
618 }
619
620 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
621 {
622     ((AliMUONchamber*) (*fChambers)[id])->SetNsec(nsec);
623 }
624
625
626 //___________________________________________
627
628 void AliMUON::StepManager()
629 {
630     printf("Dummy version of muon step -- it should never happen!!\n");
631     /*
632     const Float_t kRaddeg = 180/TMath::Pi();
633     Int_t nsec, ipart;
634     TLorentzVector x, p;
635     Float_t pt, th0, th2;
636     char proc[5];
637     if(fAccCut) {
638         if((nsec=gMC->NSecondaries())>0) {
639             gMC->ProdProcess(proc);
640             if((gMC->TrackPid()==443 || gMC->TrackPid()==553) && !strcmp(proc,"DCAY")) {
641                 //
642                 // Check angular acceptance
643                 // --- and have muons from resonance decays in the wanted window --- 
644                 if(nsec != 2) {
645                     printf(" AliMUON::StepManager: Strange resonance Decay into %d particles\n",nsec);
646                     gMC->StopEvent();
647                 } else {
648                     gMC->GetSecondary(0,ipart,x,p);
649                     pt  = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
650                     th0 = TMath::ATan2(pt,p[2])*kRaddeg;
651                     gMC->GetSecondary(1,ipart,x,p);
652                     pt  = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
653                     th2 = TMath::ATan2(pt,p[2])*kRaddeg;
654                     if(!(fAccMin < th0 && th0 < fAccMax) ||
655                        !(fAccMin < th2 && th2 < fAccMax)) 
656                         gMC->StopEvent();
657                 }
658             }
659         }
660     }
661     */
662 }
663
664 void AliMUON::MakePadHits(Float_t xhit,Float_t yhit,Float_t eloss, Int_t idvol)
665 {
666 //
667 //  Calls the charge disintegration method of the current chamber and adds
668 //  the simulated cluster to the root treee 
669 //
670     Int_t clhits[7];
671     Float_t newclust[6][500];
672     Int_t nnew;
673     
674     
675 //
676 //  Integrated pulse height on chamber
677
678     
679     clhits[0]=fNhits+1;
680 //
681 //
682     ((AliMUONchamber*) (*fChambers)[idvol])->DisIntegration(eloss, xhit, yhit, nnew, newclust);
683 //    printf("\n Add new clusters %d %f \n", nnew, eloss*1.e9);
684     Int_t ic=0;
685     
686 //
687 //  Add new clusters
688     for (Int_t i=0; i<nnew; i++) {
689         if (Int_t(newclust[3][i]) > 0) {
690             ic++;
691 // Cathode plane
692             clhits[1] = Int_t(newclust[5][i]);
693 //  Cluster Charge
694             clhits[2] = Int_t(newclust[0][i]);
695 //  Pad: ix
696             clhits[3] = Int_t(newclust[1][i]);
697 //  Pad: iy 
698             clhits[4] = Int_t(newclust[2][i]);
699 //  Pad: charge
700             clhits[5] = Int_t(newclust[3][i]);
701 //  Pad: chamber sector
702             clhits[6] = Int_t(newclust[4][i]);
703             
704             AddCluster(clhits);
705         }
706     }
707 //    printf("\n %d new clusters added", ic);
708 }
709
710 void AliMUON::Digitise(Int_t nev,Int_t bgr_ev,Option_t *option, Option_t *,Text_t *filename)
711 {
712     // keep galice.root for signal and name differently the file for 
713     // background when add! otherwise the track info for signal will be lost !
714   
715     static Bool_t first=kTRUE;
716 //    static TTree *TrH1;
717     static TFile *File;
718     char *Add = strstr(option,"Add");
719     //char *listoftracks = strstr(opt,"listoftracks");
720
721     AliMUONchamber*  iChamber;
722     AliMUONsegmentation*  segmentation;
723
724     
725     Int_t trk[50];
726     Int_t chtrk[50];  
727     TObjArray *list=new TObjArray;
728     static TClonesArray *p_adr=0;
729     if(!p_adr) p_adr=new TClonesArray("TVector",1000);
730     Int_t digits[5]; 
731
732     AliMUON *MUON  = (AliMUON *) gAlice->GetModule("MUON");
733     AliMUONHitMap * HitMap[10];
734     for (Int_t i=0; i<10; i++) {HitMap[i]=0;}
735     if (Add ) {
736         if(first) {
737             fFileName=filename;
738             cout<<"filename"<<fFileName<<endl;
739             File=new TFile(fFileName);
740             cout<<"I have opened "<<fFileName<<" file "<<endl;
741             fHits2     = new TClonesArray("AliMUONhit",1000  );
742             fClusters2 = new TClonesArray("AliMUONcluster",10000);
743         }           
744         first=kFALSE;
745         File->cd();
746         //File->ls();
747         // Get Hits Tree header from file
748         if(fHits2) fHits2->Clear();
749         if(fClusters2) fClusters2->Clear();
750         if(TrH1) delete TrH1;
751         TrH1=0;
752         
753         char treeName[20];
754         sprintf(treeName,"TreeH%d",bgr_ev);
755         TrH1 = (TTree*)gDirectory->Get(treeName);
756         //printf("TrH1 %p of treename %s for event %d \n",TrH1,treeName,bgr_ev);
757         
758         if (!TrH1) {
759             printf("ERROR: cannot find Hits Tree for event:%d\n",bgr_ev);
760         }
761         // Set branch addresses
762         TBranch *branch;
763         char branchname[20];
764         sprintf(branchname,"%s",GetName());
765         if (TrH1 && fHits2) {
766             branch = TrH1->GetBranch(branchname);
767             if (branch) branch->SetAddress(&fHits2);
768         }
769         if (TrH1 && fClusters2) {
770             branch = TrH1->GetBranch("MUONCluster");
771             if (branch) branch->SetAddress(&fClusters2);
772         }
773 // test
774         //Int_t ntracks1 =(Int_t)TrH1->GetEntries();
775         //printf("background - ntracks1 - %d\n",ntracks1);
776     }
777     //
778     // loop over cathodes
779     //
780     AliMUONHitMap* hm;
781     Int_t countadr=0;
782     for (int icat=0; icat<2; icat++) { 
783         Int_t counter=0;
784         for (Int_t i =0; i<10; i++) {
785             iChamber=(AliMUONchamber*) (*fChambers)[i];
786             if (iChamber->Nsec()==1 && icat==1) {
787                 continue;
788             } else {
789                 segmentation=iChamber->GetSegmentationModel(icat+1);
790             }
791             HitMap[i] = new AliMUONHitMapA1(segmentation, list);
792         }
793         //printf("Start loop over tracks \n");     
794 //
795 //   Loop over tracks
796 //
797
798         TTree *TH = gAlice->TreeH();
799         Int_t ntracks =(Int_t) TH->GetEntries();
800         //printf("signal - ntracks %d\n",ntracks);
801         Int_t nmuon[10]={0,0,0,0,0,0,0,0,0,0};
802         Float_t xhit[10][2];
803         Float_t yhit[10][2];
804         
805         for (Int_t track=0; track<ntracks; track++) {
806             gAlice->ResetHits();
807             TH->GetEvent(track);
808             
809 //
810 //   Loop over hits
811             for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); 
812                 mHit;
813                 mHit=(AliMUONhit*)MUON->NextHit()) 
814             {
815                 Int_t   nch   = mHit->fChamber-1;  // chamber number
816                 if (nch >9) continue;
817                 iChamber = &(MUON->Chamber(nch));
818                 Int_t rmin = (Int_t)iChamber->RInner();
819                 Int_t rmax = (Int_t)iChamber->ROuter();
820                 // new 17.07.99
821                 if (Add) {
822
823                   if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
824                     xhit[nch][nmuon[nch]]=mHit->fX;
825                     yhit[nch][nmuon[nch]]=mHit->fY;
826                     nmuon[nch]++;
827                     if (nmuon[nch] >2) printf("nmuon %d\n",nmuon[nch]);
828                     
829                   }
830                 }
831
832
833
834                 
835 //
836 // Loop over pad hits
837                 for (AliMUONcluster* mPad=
838                          (AliMUONcluster*)MUON->FirstPad(mHit,fClusters);
839                      mPad;
840                      mPad=(AliMUONcluster*)MUON->NextPad(fClusters))
841                 {
842                     Int_t cathode  = mPad->fCathode;    // cathode number
843                     Int_t ipx      = mPad->fPadX;       // pad number on X
844                     Int_t ipy      = mPad->fPadY;       // pad number on Y
845                     Int_t iqpad    = Int_t(mPad->fQpad*kScale);// charge per pad
846 //                  Int_t iqpad    = mPad->fQpad;       // charge per pad
847 //
848 //
849                     
850                     if (cathode != (icat+1)) continue;
851                     // fill the info array
852                     Float_t thex, they;
853                     segmentation=iChamber->GetSegmentationModel(cathode);
854                     segmentation->GetPadCxy(ipx,ipy,thex,they);
855                     Float_t rpad=TMath::Sqrt(thex*thex+they*they);
856                     if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
857
858                     new((*p_adr)[countadr++]) TVector(2);
859                     TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
860                     trinfo(0)=(Float_t)track;
861                     trinfo(1)=(Float_t)iqpad;
862
863                     digits[0]=ipx;
864                     digits[1]=ipy;
865                     digits[2]=iqpad;
866                     digits[3]=iqpad;
867                     if (mHit->fParticle == kMuonPlus || mHit->fParticle == kMuonMinus) {
868                     digits[4]=mPad->fHitNumber;
869                     } else digits[4]=-1;
870
871                     AliMUONlist* pdigit;
872                     // build the list of fired pads and update the info
873                     if (!HitMap[nch]->TestHit(ipx, ipy)) {
874
875                         list->AddAtAndExpand(
876                             new AliMUONlist(nch,digits),counter);
877                         
878                         HitMap[nch]->SetHit(ipx, ipy, counter);
879                         counter++;
880                         pdigit=(AliMUONlist*)list->At(list->GetLast());
881                         // list of tracks
882                         TObjArray *trlist=(TObjArray*)pdigit->TrackList();
883                         trlist->Add(&trinfo);
884                     } else {
885                         pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
886                         // update charge
887                         (*pdigit).fSignal+=iqpad;
888                         (*pdigit).fPhysics+=iqpad;                      
889                         // update list of tracks
890                         TObjArray* trlist=(TObjArray*)pdigit->TrackList();
891                         Int_t last_entry=trlist->GetLast();
892                         TVector *ptrk_p=(TVector*)trlist->At(last_entry);
893                         TVector &ptrk=*ptrk_p;
894                         Int_t last_track=Int_t(ptrk(0));
895                         Int_t last_charge=Int_t(ptrk(1));
896                         if (last_track==track) {
897                             last_charge+=iqpad;
898                             trlist->RemoveAt(last_entry);
899                             trinfo(0)=last_track;
900                             trinfo(1)=last_charge;
901                             trlist->AddAt(&trinfo,last_entry);
902                         } else {
903                             trlist->Add(&trinfo);
904                         }
905                         // check the track list
906                         Int_t nptracks=trlist->GetEntriesFast();
907                         if (nptracks > 2) {
908                             for (Int_t tr=0;tr<nptracks;tr++) {
909                                 TVector *pptrk_p=(TVector*)trlist->At(tr);
910                                 TVector &pptrk=*pptrk_p;
911                                 trk[tr]=Int_t(pptrk(0));
912                                 chtrk[tr]=Int_t(pptrk(1));
913                             }
914                         } // end if nptracks
915                     } //  end if pdigit
916                 } //end loop over clusters
917             } // hit loop
918         } // track loop
919         
920         //Int_t nentr1=list->GetEntriesFast();
921         //printf(" \n counter, nentr1 %d %d\n",counter,nentr1);
922
923         // open the file with background
924        
925         if (Add ) {
926             ntracks =(Int_t)TrH1->GetEntries();
927             //printf("background - icat,ntracks1  %d %d\n",icat,ntracks);
928             //printf("background - Start loop over tracks \n");     
929 //
930 //   Loop over tracks
931 //
932             for (Int_t track=0; track<ntracks; track++) {
933
934                 if (fHits2)       fHits2->Clear();
935                 if (fClusters2)   fClusters2->Clear();
936
937                 TrH1->GetEvent(track);
938 //
939 //   Loop over hits
940                 AliMUONhit* mHit;
941                 for(int i=0;i<fHits2->GetEntriesFast();++i) 
942         {       
943                     mHit=(AliMUONhit*) (*fHits2)[i];
944                     Int_t   nch   = mHit->fChamber-1;  // chamber number
945                     if (nch >9) continue;
946                     iChamber = &(MUON->Chamber(nch));
947                     Int_t rmin = (Int_t)iChamber->RInner();
948                     Int_t rmax = (Int_t)iChamber->ROuter();
949                     Float_t xbgr=mHit->fX;
950                     Float_t ybgr=mHit->fY;
951                     Bool_t cond=kFALSE;
952                     
953                     for (Int_t imuon =0; imuon < nmuon[nch]; imuon++) {
954                         Float_t dist= (xbgr-xhit[nch][imuon])*(xbgr-xhit[nch][imuon])
955                             +(ybgr-yhit[nch][imuon])*(ybgr-yhit[nch][imuon]);
956                         if (dist<100) cond=kTRUE;
957                     }
958                     if (!cond) continue;
959                     
960 //
961 // Loop over pad hits
962                     for (AliMUONcluster* mPad=
963                              (AliMUONcluster*)MUON->FirstPad(mHit,fClusters2);
964                          mPad;
965                          mPad=(AliMUONcluster*)MUON->NextPad(fClusters2))
966                     {
967
968                         Int_t cathode  = mPad->fCathode;    // cathode number
969                         Int_t ipx      = mPad->fPadX;       // pad number on X
970                         Int_t ipy      = mPad->fPadY;       // pad number on Y
971                         Int_t iqpad    = Int_t(mPad->fQpad*kScale);// charge per pad
972 //                      Int_t iqpad    = mPad->fQpad;       // charge per pad
973
974                         if (cathode != (icat+1)) continue;
975                         //if (!HitMap[nch]->CheckBoundary()) continue;
976                         // fill the info array
977                         Float_t thex, they;
978                         segmentation=iChamber->GetSegmentationModel(cathode);
979                         segmentation->GetPadCxy(ipx,ipy,thex,they);
980                         Float_t rpad=TMath::Sqrt(thex*thex+they*they);
981                         if (rpad < rmin || iqpad ==0 || rpad > rmax) continue;
982
983                             new((*p_adr)[countadr++]) TVector(2);
984                             TVector &trinfo=*((TVector*) (*p_adr)[countadr-1]);
985                             trinfo(0)=-1;  // tag background
986                             trinfo(1)=-1;
987
988                         digits[0]=ipx;
989                         digits[1]=ipy;
990                         digits[2]=iqpad;
991                         digits[3]=0;
992                         digits[4]=-1;
993
994                         AliMUONlist* pdigit;
995                         // build the list of fired pads and update the info
996                         if (!HitMap[nch]->TestHit(ipx, ipy)) {
997                             list->AddAtAndExpand(new AliMUONlist(nch,digits),counter);
998                         
999                             HitMap[nch]->SetHit(ipx, ipy, counter);
1000                             counter++;
1001                             
1002                             pdigit=(AliMUONlist*)list->At(list->GetLast());
1003                             // list of tracks
1004                                 TObjArray *trlist=(TObjArray*)pdigit->
1005                                                    TrackList();
1006                                 trlist->Add(&trinfo);
1007                         } else {
1008                             pdigit=(AliMUONlist*) HitMap[nch]->GetHit(ipx, ipy);
1009                             // update charge
1010                             (*pdigit).fSignal+=iqpad;
1011
1012                             // update list of tracks
1013                                 TObjArray* trlist=(TObjArray*)pdigit->
1014                                                    TrackList();
1015                                 Int_t last_entry=trlist->GetLast();
1016                                 TVector *ptrk_p=(TVector*)trlist->
1017                                                  At(last_entry);
1018                                 TVector &ptrk=*ptrk_p;
1019                                 Int_t last_track=Int_t(ptrk(0));
1020                                 if (last_track==-1) {
1021                                     continue;
1022                                 } else {
1023                                     trlist->Add(&trinfo);
1024                                 }
1025                                 // check the track list
1026                                 Int_t nptracks=trlist->GetEntriesFast();
1027                                 if (nptracks > 0) {
1028                                     for (Int_t tr=0;tr<nptracks;tr++) {
1029                                         TVector *pptrk_p=(TVector*)trlist->At(tr);
1030                                         TVector &pptrk=*pptrk_p;
1031                                         trk[tr]=Int_t(pptrk(0));
1032                                         chtrk[tr]=Int_t(pptrk(1));
1033                                     }
1034                                 } // end if nptracks
1035                         } //  end if pdigit
1036                     } //end loop over clusters
1037                 } // hit loop
1038             } // track loop
1039             //Int_t nentr2=list->GetEntriesFast();
1040             //printf(" \n counter2, nentr2 %d %d \n",counter,nentr2);
1041             TTree *fAli=gAlice->TreeK();
1042             TFile *file;
1043             
1044             if (fAli) file =fAli->GetCurrentFile();
1045             file->cd();
1046         } // if Add     
1047
1048         Int_t tracks[10];
1049         Int_t charges[10];
1050         //cout<<"start filling digits \n "<<endl;
1051         //      const Float_t zero_supm =    6.;
1052         Int_t nentries=list->GetEntriesFast();
1053         //printf(" \n \n nentries %d \n",nentries);
1054         // start filling the digits
1055         
1056         for (Int_t nent=0;nent<nentries;nent++) {
1057             AliMUONlist *address=(AliMUONlist*)list->At(nent);
1058             if (address==0) continue; 
1059             Int_t ich=address->fChamber;
1060             Int_t q=address->fSignal; 
1061             iChamber=(AliMUONchamber*) (*fChambers)[ich];
1062             AliMUONresponse * response=iChamber->GetResponseModel();
1063             Int_t adcmax= (Int_t) response->MaxAdc();
1064             // add white noise and do zero-suppression and signal truncation
1065             Float_t MeanNoise = gRandom->Gaus(1, 0.2);
1066             Float_t Noise     = gRandom->Gaus(0, MeanNoise);
1067             q+=(Int_t)Noise; 
1068             if (address->fPhysics !=0 ) address->fPhysics+=(Int_t)Noise; 
1069             if ( q <= zero_supm ) continue;
1070             if ( q > adcmax)  q=adcmax;
1071             digits[0]=address->fPadX;
1072             digits[1]=address->fPadY;
1073             digits[2]=q;
1074             digits[3]=address->fPhysics;
1075             digits[4]=address->fHit;
1076             //printf("fSignal, fPhysics fTrack %d %d %d \n",digits[2],digits[3],digits[4]);
1077             
1078             TObjArray* trlist=(TObjArray*)address->TrackList();
1079             Int_t nptracks=trlist->GetEntriesFast();
1080             //printf("nptracks, trlist   %d  %p\n",nptracks,trlist);
1081
1082                 // this was changed to accomodate the real number of tracks
1083                 if (nptracks > 10) {
1084                     cout<<"Attention - nptracks > 10 "<<nptracks<<endl;
1085                     nptracks=10;
1086                 }
1087                 if (nptracks > 2) {
1088                     printf("Attention - nptracks > 2  %d \n",nptracks);
1089                     printf("cat,ich,ix,iy,q %d %d %d %d %d \n",icat,ich,digits[0],digits[1],q);
1090                 }
1091                 for (Int_t tr=0;tr<nptracks;tr++) {
1092                     TVector *pp_p=(TVector*)trlist->At(tr);
1093                     if(!pp_p ) printf("pp_p - %p\n",pp_p);
1094                     TVector &pp  =*pp_p;
1095                     tracks[tr]=Int_t(pp(0));
1096                     charges[tr]=Int_t(pp(1));
1097                 //printf("tracks, charges - %d %d\n",tracks[tr],charges[tr]);
1098                 }      //end loop over list of tracks for one pad
1099             // Sort list of tracks according to charge
1100                 if (nptracks > 1) {
1101                     SortTracks(tracks,charges,nptracks);
1102                 }
1103                 if (nptracks < 10 ) {
1104                     for (Int_t i=nptracks; i<10; i++) {
1105                         tracks[i]=0;
1106                         charges[i]=0;
1107                     }
1108                 }
1109
1110             // fill digits
1111             MUON->AddDigits(ich,tracks,charges,digits);
1112         }
1113         //cout<<"I'm out of the loops for digitisation"<<endl;
1114         gAlice->TreeD()->Fill();
1115         TTree *TD=gAlice->TreeD();
1116
1117         Stat_t ndig=TD->GetEntries();
1118         cout<<"number of digits  "<<ndig<<endl;
1119         TClonesArray *fDch;
1120         for (int k=0;k<10;k++) {
1121             fDch= MUON->DigitsAddress(k);
1122             int ndig=fDch->GetEntriesFast();
1123             printf (" i, ndig %d %d \n",k,ndig);
1124         }
1125
1126         MUON->ResetDigits();
1127         list->Delete();
1128         for(Int_t ii=0;ii<10;++ii) {
1129             if (HitMap[ii]) {
1130                 hm=HitMap[ii];
1131                 delete hm;
1132                 HitMap[ii]=0;
1133             }
1134         }
1135         
1136     } //end loop over cathodes
1137
1138        char hname[30];
1139        sprintf(hname,"TreeD%d",nev);
1140        gAlice->TreeD()->Write(hname);
1141        // reset tree
1142        gAlice->TreeD()->Reset();
1143        delete list;
1144        //Int_t nadr=p_adr->GetEntriesFast();
1145        // printf(" \n \n nadr %d \n",nadr);
1146
1147        p_adr->Clear();
1148        // gObjectTable->Print();
1149        
1150 }
1151
1152 void AliMUON::SortTracks(Int_t *tracks,Int_t *charges,Int_t ntr)
1153 {
1154   //
1155   // Sort the list of tracks contributing to a given digit
1156   // Only the 3 most significant tracks are acctually sorted
1157   //
1158   
1159   //
1160   //  Loop over signals, only 3 times
1161   //
1162   
1163   Int_t qmax;
1164   Int_t jmax;
1165   Int_t idx[3] = {-2,-2,-2};
1166   Int_t jch[3] = {-2,-2,-2};
1167   Int_t jtr[3] = {-2,-2,-2};
1168   Int_t i,j,imax;
1169   
1170   if (ntr<3) imax=ntr;
1171   else imax=3;
1172   for(i=0;i<imax;i++){
1173     qmax=0;
1174     jmax=0;
1175     
1176     for(j=0;j<ntr;j++){
1177       
1178       if((i == 1 && j == idx[i-1]) 
1179          ||(i == 2 && (j == idx[i-1] || j == idx[i-2]))) continue;
1180       
1181       if(charges[j] > qmax) {
1182         qmax = charges[j];
1183         jmax=j;
1184       }       
1185     } 
1186     
1187     if(qmax > 0) {
1188       idx[i]=jmax;
1189       jch[i]=charges[jmax]; 
1190       jtr[i]=tracks[jmax]; 
1191     }
1192     
1193   } 
1194   
1195   for(i=0;i<3;i++){
1196     if (jtr[i] == -2) {
1197          charges[i]=0;
1198          tracks[i]=0;
1199     } else {
1200          charges[i]=jch[i];
1201          tracks[i]=jtr[i];
1202     }
1203   }
1204
1205 }
1206
1207 void AliMUON::FindClusters(Int_t nev,Int_t last_entry)
1208 {
1209
1210 //
1211 // Loop on chambers and on cathode planes
1212 //
1213   for (Int_t icat=0;icat<2;icat++) {
1214             gAlice->ResetDigits();
1215             gAlice->TreeD()->GetEvent(last_entry+icat); // spurious +1 ...
1216             if (nev < 10) printf("last_entry , icat - %d %d \n",last_entry,icat);
1217             //gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
1218
1219       for (Int_t ich=0;ich<10;ich++) {
1220           AliMUONchamber* iChamber=(AliMUONchamber*) (*fChambers)[ich];
1221           TClonesArray *MUONdigits  = this->DigitsAddress(ich);
1222           if (MUONdigits == 0) continue;
1223           //
1224           // Get ready the current chamber stuff
1225           //
1226           AliMUONresponse* response = iChamber->GetResponseModel();
1227           AliMUONsegmentation*  seg = iChamber->GetSegmentationModel(icat+1);
1228           AliMUONClusterFinder* rec = iChamber->GetReconstructionModel();
1229           //printf("icat, ich, seg - %d %d %p\n",icat,ich,seg);
1230           if (seg) {      
1231               rec->SetSegmentation(seg);
1232               rec->SetResponse(response);
1233               rec->SetDigits(MUONdigits);
1234               rec->SetChamber(ich);
1235               if (nev==0) rec->CalibrateCOG(); 
1236               rec->FindRawClusters();
1237           }  
1238           //printf("Finish FindRawClusters for cathode %d in chamber %d\n",icat,ich);
1239           
1240           TClonesArray *fRch;
1241           fRch=RawClustAddress(ich);
1242           fRch->Sort();
1243           // it seems to work 
1244          
1245
1246       } // for ich
1247       // fill the tree
1248       TTree *TR=gAlice->TreeR();
1249
1250       gAlice->TreeR()->Fill();
1251
1252       Stat_t nent=TR->GetEntries();
1253       cout<<"number of entries  "<<nent<<endl;
1254       TClonesArray *fRch;
1255       for (int i=0;i<10;i++) {
1256           fRch=RawClustAddress(i);
1257           int nraw=fRch->GetEntriesFast();
1258           printf (" i, nraw %d %d \n",i,nraw);
1259       }
1260       ResetRawClusters();
1261
1262   } // for icat
1263
1264   char hname[30];
1265   sprintf(hname,"TreeR%d",nev);
1266   gAlice->TreeR()->Write(hname);
1267   gAlice->TreeR()->Reset();
1268
1269   //gObjectTable->Print();
1270
1271 }
1272  
1273 //______________________________________________________________________________
1274 //_____________________________________________________________________________ 
1275 void AliMUON::CathodeCorrelation(Int_t nev)
1276 {
1277
1278 // Correlates the clusters on the two cathode planes and build a list of
1279 // other possible combinations (potential ghosts) - for the moment use the
1280 // criteria of minimum distance between the CoGs of the two correlated
1281 // clusters
1282
1283
1284 //
1285 // Loop on chambers and on clusters on the cathode plane with the highest
1286 // number of clusters
1287
1288     static Bool_t first=kTRUE;
1289
1290      AliMUONRawCluster  *mRaw1;
1291      AliMUONRawCluster  *mRaw2;
1292      AliMUONchamber     *iChamber;
1293      AliMUONsegmentation *seg;
1294      TArrayF x1, y1, x2, y2, q1, q2;
1295      x1.Set(5000);
1296      x2.Set(5000);     
1297      y1.Set(5000);
1298      y2.Set(5000);     
1299      q1.Set(5000);
1300      q2.Set(5000);     
1301      
1302 // Get pointers to Alice detectors and Digits containers
1303      TTree *TR = gAlice->TreeR();
1304      Int_t nent=(Int_t)TR->GetEntries();
1305      if (nev < 10) printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
1306   
1307
1308      Int_t idx[4]; 
1309      Float_t xc2[4],yc2[4];
1310      Float_t xrec2, yrec2;
1311      Float_t xd0, xdif, ydif;
1312      Float_t ysrch,xd,xmax,ymax;
1313      Int_t ilow, iup, iraw1, i;
1314      //
1315      Float_t xarray[50];
1316      Float_t xdarray[50];
1317      Float_t yarray[50];
1318      Float_t qarray[50];
1319      Int_t idx2[50];
1320
1321      // Int_t nraw[2], entry,cathode;
1322
1323      for (i=0;i<50;i++) {
1324          xdarray[i]=1100.;
1325          xarray[i]=0.;
1326          yarray[i]=0.;
1327          qarray[i]=0.;
1328          idx2[i]=-1;
1329      }
1330      for (i=0;i<4;i++) {
1331           idx[i]=-1;
1332           xc2[i]=0.;
1333           yc2[i]=0.;
1334      }
1335
1336      // access to the Raw Clusters tree
1337      for (Int_t ich=0;ich<10;ich++) {
1338          iChamber = &(Chamber(ich));
1339          TClonesArray *MUONrawclust  = RawClustAddress(ich);
1340          ResetRawClusters();
1341          TR->GetEvent(nent-2);
1342          //TR->GetEvent(1);
1343          Int_t nrawcl1 = MUONrawclust->GetEntries();
1344          // printf("Found %d raw clusters for cathode 1 in chamber %d \n"
1345          //      ,nrawcl1,ich+1);
1346          if (!nrawcl1) continue;
1347
1348          seg = iChamber->GetSegmentationModel(1);
1349          // loop over raw clusters of first cathode
1350          for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1351                  mRaw1= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw1);
1352                  x1[iraw1]=mRaw1->fX;
1353                  y1[iraw1]=mRaw1->fY;
1354                  q1[iraw1]=(Float_t)mRaw1->fQ; //maybe better fPeakSignal
1355          } // rawclusters cathode 1
1356 //
1357          // Get information from 2nd cathode
1358          ResetRawClusters();
1359          TR->GetEvent(nent-1);
1360          //TR->GetEvent(2);
1361          Int_t nrawcl2 = MUONrawclust->GetEntries();
1362          if (!nrawcl2) {
1363              for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1364                  idx[3]=iraw1;
1365                  xc2[3]=x1[iraw1];
1366                  yc2[3]=y1[iraw1];
1367                  //printf("nrawcl2 is zero - idx[0] %d\n",idx[0]);
1368                  
1369                  AddCathCorrel(ich,idx,xc2,yc2);
1370                  // reset
1371                  idx[3]=-1;
1372                  xc2[3]=0.;
1373                  yc2[3]=0.;
1374                  
1375              } // store information from cathode 1 only 
1376          } else {
1377            //  printf("Found %d raw clusters for cathode 2 in chamber %d \n",
1378            // nrawcl2, ich+1);
1379
1380              for (Int_t iraw2=0; iraw2<nrawcl2; iraw2++) {
1381                  mRaw2= (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw2);
1382                  x2[iraw2]=mRaw2->fX;
1383                  y2[iraw2]=mRaw2->fY;   
1384                  q2[iraw2]=(Float_t)mRaw2->fQ;  
1385              } // rawclusters cathode 2
1386 //
1387 // Initalisation finished
1388              for (iraw1=0; iraw1<nrawcl1; iraw1++) {
1389              // find the sector
1390                  Int_t ix,iy;
1391                  seg->GetPadIxy(x1[iraw1],y1[iraw1],ix,iy);   
1392                  Int_t isec=seg->Sector(ix,iy);
1393                  // range to look for ghosts ?!
1394                  if (ich < 5) {
1395                      ymax = seg->Dpy(isec)*7/2;
1396                      xmax = seg->Dpx(isec)*7/2;
1397                  } else {
1398                      ymax = seg->Dpy(isec)*13/2;
1399                      xmax = seg->Dpx(isec)*3/2;
1400                  }
1401                  ysrch=ymax+y1[iraw1];
1402                  
1403                  ilow = AliMUONRawCluster::
1404                      BinarySearch(ysrch-2*ymax,y2,0,nrawcl2+1);
1405                  iup=   AliMUONRawCluster::
1406                      BinarySearch(ysrch,y2,ilow,nrawcl2+1);
1407                  if (ilow<0 || iup <0 || iup>nrawcl2) continue;
1408                  Int_t counter=0;
1409                  for (Int_t iraw2=ilow; iraw2<=iup; iraw2++) {
1410                      xrec2=x2[iraw2];
1411                      yrec2=y2[iraw2];   
1412                      xdif=x1[iraw1]-xrec2;
1413                      ydif=y1[iraw1]-yrec2;
1414                      xd=TMath::Sqrt(xdif*xdif+ydif*ydif);
1415                      if (iraw2==ilow) { 
1416                          if (ilow==iup) 
1417                              xd0=TMath::
1418                              Sqrt(2*xmax*2*xmax+2*ymax*2*ymax);
1419                          else xd0=101.; 
1420                      } 
1421                      Float_t qdif=TMath::Abs(q1[iraw1]-q2[iraw2])/q1[iraw1];
1422                      
1423                      if (x1[iraw1]*xrec2 > 0) {
1424                          if (xd <= xd0 )  {
1425 //                           printf("q1, q2 qdif % f %f %f \n",q1[iraw1],q2[iraw2],qdif);
1426 //                           printf("x1, x2 y1 y2 % f %f %f %f \n",x1[iraw1],xrec2,y1[iraw1],yrec2);
1427                            //if (qdif <0.3) { //check this number
1428                                  
1429                                  xd0=xd;
1430                                  idx2[counter]=iraw2;
1431                                  xdarray[counter]=xd;
1432                                  xarray[counter]=xdif;
1433                                  yarray[counter]=ydif;
1434                                  qarray[counter]=qdif;
1435                                  counter++;
1436                            // }
1437                              
1438                          }
1439                      } // check for same quadrant
1440                  } // loop over 2nd cathode range 
1441                  
1442                  
1443                  if (counter >=2) {
1444                      AliMUONRawCluster::
1445                          SortMin(idx2,xdarray,xarray,yarray,qarray,counter);
1446                      if (xdarray[0]<seg->Dpx(isec) && xdarray[1]<seg->Dpx(isec)) {
1447                          if (qarray[0]>qarray[1]){
1448                              Int_t swap=idx2[0];
1449                              idx2[0]=idx2[1];
1450                              idx2[1]=swap;
1451                          }
1452                      }
1453                  }
1454                  int imax;
1455                  if (counter <3) imax=counter;
1456                  else imax=3;
1457
1458                  for (int i=0;i<imax;i++) {
1459                      if (idx2[i] >= 0 && idx2[i] < nrawcl2) {
1460                          if (xarray[i] > xmax || yarray[i] > 2*ymax) 
1461                              continue;
1462                          idx[i]=idx2[i];
1463                          xc2[i]=x2[idx2[i]];
1464                          yc2[i]=y2[idx2[i]];
1465                      }
1466                  }
1467                  // add info about the cluster on the 'starting' cathode
1468
1469                  idx[3]=iraw1;
1470                  xc2[3]=x1[iraw1];
1471                  yc2[3]=y1[iraw1];
1472                  //if (idx[0] <0)  printf("iraw1 imax idx2[0] idx[0] %d %d %d %d\n",iraw1,imax,idx2[0],idx[0]);
1473                  AddCathCorrel(ich,idx,xc2,yc2);
1474                  // reset
1475                  for (Int_t ii=0;ii<counter;ii++) {
1476                      xdarray[ii]=1100.;
1477                      xarray[ii]=0.;
1478                      yarray[ii]=0.;
1479                      qarray[ii]=0.;
1480                      idx2[ii]=-1;
1481                  }
1482                  for (Int_t iii=0;iii<3;iii++) {
1483                      idx[iii]=-1;
1484                      xc2[iii]=0.;
1485                      yc2[iii]=0.;
1486                  }
1487              } // iraw1
1488          }
1489          x1.Reset();
1490          x2.Reset();     
1491          y1.Reset();
1492          y2.Reset();     
1493          q1.Reset();
1494          q2.Reset();     
1495      } //ich
1496 // 
1497      if (first) {
1498          MakeTreeC("C");
1499          first=kFALSE;
1500      }
1501      TTree *TC=TreeC();
1502      TC->Fill();
1503      //Int_t nentries=(Int_t)TC->GetEntries();
1504     //cout<<"number entries in tree of correlated clusters  "<<nentries<<endl;
1505      TClonesArray *fCch;
1506      static Int_t countev=0;
1507      Int_t countch=0;
1508
1509      for (Int_t ii=0;ii<10;ii++) {
1510            fCch= CathCorrelAddress(ii);
1511            Int_t ncor=fCch->GetEntriesFast();
1512            printf (" ii, ncor %d %d \n",ii,ncor);
1513            if (ncor>=2) countch++;
1514      }
1515
1516      // write
1517      char hname[30];
1518      sprintf(hname,"TreeC%d",nev);
1519      TC->Write(hname);
1520      // reset tree
1521      ResetCorrelation();
1522      TC->Reset();
1523
1524      if (countch==10) countev++;
1525      printf("countev - %d\n",countev);
1526     
1527 //     gObjectTable->Print();
1528      
1529      
1530 }
1531
1532
1533 //_____________________________________________________________________________
1534
1535 void AliMUON::MakeTreeC(Option_t *option)
1536 {
1537      char *C = strstr(option,"C");
1538      if (C && !fTreeC) fTreeC = new TTree("TC","CathodeCorrelation");
1539
1540 //  Create a branch for correlation 
1541
1542      const Int_t buffersize = 4000;
1543      char branchname[30];
1544
1545 // one branch for correlation per chamber
1546      for (int i=0; i<10 ;i++) {
1547          sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
1548       
1549          if (fCathCorrel   && fTreeC) {
1550             TreeC()->Branch(branchname,&((*fCathCorrel)[i]), buffersize);
1551             printf("Making Branch %s for correlation in chamber %d\n",branchname,i+1);
1552          }      
1553      }
1554 }
1555
1556 //_____________________________________________________________________________
1557 void AliMUON::GetTreeC(Int_t event)
1558 {
1559
1560     // set the branch address
1561     char treeName[20];
1562     char branchname[30];
1563
1564     ResetCorrelation();
1565     if (fTreeC) {
1566           delete fTreeC;
1567     }
1568
1569     sprintf(treeName,"TreeC%d",event);
1570     fTreeC = (TTree*)gDirectory->Get(treeName);
1571
1572
1573     TBranch *branch;
1574     if (fTreeC) {
1575         for (int i=0; i<10; i++) {
1576             sprintf(branchname,"%sCorrelation%d",GetName(),i+1);
1577             if (fCathCorrel) {
1578                 branch = fTreeC->GetBranch(branchname);
1579                 if (branch) branch->SetAddress(&((*fCathCorrel)[i]));
1580             }
1581         }
1582     } else {
1583         printf("ERROR: cannot find CathodeCorrelation Tree for event:%d\n",event);
1584     }
1585
1586     // gObjectTable->Print();
1587
1588 }
1589
1590
1591 void AliMUON::Streamer(TBuffer &R__b)
1592 {
1593    // Stream an object of class AliMUON.
1594       AliMUONchamber       *iChamber;
1595       AliMUONsegmentation  *segmentation;
1596       AliMUONresponse      *response;
1597       TClonesArray         *digitsaddress;
1598       TClonesArray         *rawcladdress;
1599       TClonesArray         *corcladdress;
1600       //      TObjArray            *clustaddress;
1601       
1602    if (R__b.IsReading()) {
1603       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1604       AliDetector::Streamer(R__b);
1605       R__b >> fNclusters;
1606       R__b >> fClusters; // diff
1607       R__b >> fDchambers;
1608       R__b >> fRawClusters;
1609       R__b >> fCathCorrel;
1610       R__b.ReadArray(fNdch);
1611       R__b.ReadArray(fNrawch);
1612       R__b.ReadArray(fNcorch);
1613       //
1614       R__b >> fAccCut;
1615       R__b >> fAccMin;
1616       R__b >> fAccMax; 
1617       //   
1618       // modifs perso  
1619       R__b >> fSPxzCut;  
1620       R__b >> fSSigmaCut;
1621       R__b >> fSXPrec; 
1622       R__b >> fSYPrec;
1623       //
1624       R__b >> fChambers;
1625 // Stream chamber related information
1626       for (Int_t i =0; i<10; i++) {
1627           iChamber=(AliMUONchamber*) (*fChambers)[i];
1628           iChamber->Streamer(R__b);
1629           if (iChamber->Nsec()==1) {
1630               segmentation=iChamber->GetSegmentationModel(1);
1631               segmentation->Streamer(R__b);
1632           } else {
1633               segmentation=iChamber->GetSegmentationModel(1);
1634               segmentation->Streamer(R__b);
1635               segmentation=iChamber->GetSegmentationModel(2);
1636               segmentation->Streamer(R__b);
1637           }
1638           response=iChamber->GetResponseModel();
1639           response->Streamer(R__b);       
1640           digitsaddress=(TClonesArray*) (*fDchambers)[i];
1641           digitsaddress->Streamer(R__b);
1642           rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1643           rawcladdress->Streamer(R__b);
1644           corcladdress=(TClonesArray*) (*fCathCorrel)[i];
1645           corcladdress->Streamer(R__b);
1646       }
1647       
1648    } else {
1649       R__b.WriteVersion(AliMUON::IsA());
1650       AliDetector::Streamer(R__b);
1651       R__b << fNclusters;
1652       R__b << fClusters; // diff
1653       R__b << fDchambers;
1654       R__b << fRawClusters;
1655       R__b << fCathCorrel;
1656       R__b.WriteArray(fNdch, 10);
1657       R__b.WriteArray(fNrawch, 10);
1658       R__b.WriteArray(fNcorch, 10);
1659       //
1660       R__b << fAccCut;
1661       R__b << fAccMin;
1662       R__b << fAccMax; 
1663       //   
1664       // modifs perso  
1665       R__b << fSPxzCut;  
1666       R__b << fSSigmaCut;
1667       R__b << fSXPrec; 
1668       R__b << fSYPrec;
1669       //
1670       R__b << fChambers;
1671 //  Stream chamber related information
1672       for (Int_t i =0; i<10; i++) {
1673           iChamber=(AliMUONchamber*) (*fChambers)[i];
1674           iChamber->Streamer(R__b);
1675           if (iChamber->Nsec()==1) {
1676               segmentation=iChamber->GetSegmentationModel(1);
1677               segmentation->Streamer(R__b);
1678           } else {
1679               segmentation=iChamber->GetSegmentationModel(1);
1680               segmentation->Streamer(R__b);
1681               segmentation=iChamber->GetSegmentationModel(2);
1682               segmentation->Streamer(R__b);
1683           }
1684           response=iChamber->GetResponseModel();
1685           response->Streamer(R__b);
1686           digitsaddress=(TClonesArray*) (*fDchambers)[i];
1687           digitsaddress->Streamer(R__b);
1688           rawcladdress=(TClonesArray*) (*fRawClusters)[i];
1689           rawcladdress->Streamer(R__b);
1690           corcladdress=(TClonesArray*) (*fCathCorrel)[i];
1691           corcladdress->Streamer(R__b);
1692       }
1693    }
1694 }
1695 AliMUONcluster* AliMUON::FirstPad(AliMUONhit*  hit, TClonesArray *clusters) 
1696 {
1697 //
1698     // Initialise the pad iterator
1699     // Return the address of the first padhit for hit
1700     TClonesArray *theClusters = clusters;
1701     Int_t nclust = theClusters->GetEntriesFast();
1702     if (nclust && hit->fPHlast > 0) {
1703         sMaxIterPad=hit->fPHlast;
1704         sCurIterPad=hit->fPHfirst;
1705         return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
1706     } else {
1707         return 0;
1708     }
1709 }
1710
1711 AliMUONcluster* AliMUON::NextPad(TClonesArray *clusters) 
1712 {
1713     sCurIterPad++;
1714     if (sCurIterPad <= sMaxIterPad) {
1715         return (AliMUONcluster*) clusters->UncheckedAt(sCurIterPad-1);
1716     } else {
1717         return 0;
1718     }
1719 }
1720
1721 //////////////////////////// modifs perso ///////////////
1722
1723 static TTree *ntuple_global;
1724 static TFile *hfile_global;
1725
1726 // variables of the tracking ntuple 
1727 struct { 
1728   Int_t ievr;           // number of event 
1729   Int_t ntrackr;        // number of tracks per event
1730   Int_t istatr[500];    // 1 = good muon, 2 = ghost, 0 = something else
1731   Int_t isignr[500];    // sign of the track
1732   Float_t pxr[500];     // x momentum of the reconstructed track
1733   Float_t pyr[500];     // y momentum of the reconstructed track
1734   Float_t pzr[500];     // z momentum of the reconstructed track
1735   Float_t zvr[500];     // z vertex 
1736   Float_t chi2r[500];   // chi2 of the fit of the track with the field map
1737   Float_t pxv[500];     // x momentum at vertex
1738   Float_t pyv[500];     // y momentum at vertex
1739   Float_t pzv[500];     // z momentum at vertex
1740 } ntuple_st;
1741
1742 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
1743 {
1744     TClonesArray *MUONrawclust  = RawClustAddress(ichamber);
1745     ResetRawClusters();
1746     TTree *TR = gAlice->TreeR();
1747     Int_t nent=(Int_t)TR->GetEntries();
1748     TR->GetEvent(nent-2+icathod-1);
1749     //TR->GetEvent(icathod);
1750     //Int_t nrawcl = (Int_t)MUONrawclust->GetEntriesFast();
1751
1752     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(icluster);
1753     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
1754     
1755     return  mRaw;
1756 }
1757
1758 void AliMUON::Reconst(Int_t &ifit, Int_t &idebug, Int_t bgd_ev, Int_t &nev, Int_t &idres, Int_t &ireadgeant, Option_t *option,Text_t *filename)
1759 {
1760   //
1761   // open kine and hits tree of background file for reconstruction of geant hits 
1762   // call tracking fortran program
1763   static Bool_t first=kTRUE;
1764   static TFile *File;
1765   char *Add = strstr(option,"Add");
1766   
1767   if (Add ) { // only in case of background with geant hits 
1768     if(first) {
1769       fFileName=filename;
1770       cout<<"filename  "<<fFileName<<endl;
1771       File=new TFile(fFileName);
1772       cout<<"I have opened "<<fFileName<<" file "<<endl;
1773       fHits2     = new TClonesArray("AliMUONhit",1000  );
1774       fParticles2 = new TClonesArray("TParticle",1000);
1775       first=kFALSE;
1776     }
1777     File->cd();
1778     if(fHits2) fHits2->Clear();
1779     if(fParticles2) fParticles2->Clear();
1780     if(TrH1) delete TrH1;
1781     TrH1=0;
1782     if(TK1) delete TK1;
1783     TK1=0;
1784     // Get Hits Tree header from file
1785     char treeName[20];
1786     sprintf(treeName,"TreeH%d",bgd_ev);
1787     TrH1 = (TTree*)gDirectory->Get(treeName);
1788     if (!TrH1) {
1789       printf("ERROR: cannot find Hits Tree for event:%d\n",bgd_ev);
1790     }
1791     // set branch addresses
1792     TBranch *branch;
1793     char branchname[30];
1794     sprintf(branchname,"%s",GetName());
1795     if (TrH1 && fHits2) {
1796       branch = TrH1->GetBranch(branchname);
1797       if (branch) branch->SetAddress(&fHits2);
1798     }
1799     TrH1->GetEntries();
1800     // get the Kine tree
1801     sprintf(treeName,"TreeK%d",bgd_ev);
1802     TK1 = (TTree*)gDirectory->Get(treeName);
1803     if (!TK1) {
1804       printf("ERROR: cannot find Kine Tree for event:%d\n",bgd_ev);
1805     }
1806     // set branch addresses
1807     if (TK1) 
1808       TK1->SetBranchAddress("Particles", &fParticles2);
1809     TK1->GetEvent(0);
1810     
1811     // get back to the first file
1812     TTree *TK = gAlice->TreeK();
1813     TFile *file1 = 0;
1814     if (TK) file1 = TK->GetCurrentFile();
1815     file1->cd();
1816     
1817   } // end if Add
1818   
1819   // call tracking fortran program
1820   reconstmuon(ifit,idebug,nev,idres,ireadgeant);
1821 }
1822
1823
1824 void AliMUON::InitTracking(Double_t &seff, Double_t &sb0, Double_t &sbl3)
1825 {
1826   //
1827   // introduce in fortran program somme parameters and cuts for tracking 
1828   // create output file "reconst.root" (histos + ntuple)
1829   cutpxz(fSPxzCut);          // Pxz cut (GeV/c) to begin the track finding
1830   sigmacut(fSSigmaCut);      // Number of sigmas delimiting the searching areas
1831   xpreci(fSXPrec);           // Chamber precision in X (cm) 
1832   ypreci(fSYPrec);           // Chamber precision in Y (cm)
1833   reco_init(seff,sb0,sbl3);
1834 }
1835
1836 void AliMUON::FinishEvent()
1837 {
1838     TTree *TK = gAlice->TreeK();
1839     if (TK) {
1840       TFile *file1 = TK->GetCurrentFile();
1841       if(file1) file1->cd();
1842     }
1843 }
1844
1845 void AliMUON::CloseTracking()
1846 {
1847   //
1848   // write histos and ntuple to "reconst.root" file
1849     reco_term();
1850 }
1851
1852 void chfill(Int_t &id, Float_t &x, Float_t &, Float_t &)
1853 {
1854   //
1855   // fill histo like hfill in fortran
1856     char name[5];
1857     sprintf(name,"h%d",id);
1858     TH1F *h1 = (TH1F*) gDirectory->Get(name);
1859     h1->Fill(x);
1860 }
1861
1862 void chfill2(Int_t &id, Float_t &x, Float_t &y, Float_t &w)
1863 {
1864   //
1865   // fill histo like hfill2 in fortran
1866     char name[5];
1867     sprintf(name,"h%d",id);
1868     TH2F *h2 = (TH2F*) gDirectory->Get(name);
1869     h2->Fill(x,y,w);
1870 }
1871
1872 void chf1(Int_t &id, Float_t &x, Float_t &w)
1873 {
1874   //
1875   // fill histo like hf1 in fortran
1876     char name[5];
1877     sprintf(name,"h%d",id);
1878     TH1F *h1 = (TH1F*) gDirectory->Get(name);
1879     h1->Fill(x,w);
1880 }
1881
1882 void hist_create()
1883 {
1884   //
1885   // Create an output file ("reconst.root")
1886   // Create some histograms and an ntuple
1887
1888     hfile_global = new TFile("reconst.root","RECREATE","Ntuple - reconstruction");
1889
1890    ntuple_global = new TTree("ntuple","Reconst ntuple");
1891    ntuple_global->Branch("ievr",&ntuple_st.ievr,"ievr/I");
1892    ntuple_global->Branch("ntrackr",&ntuple_st.ntrackr,"ntrackr/I");
1893    ntuple_global->Branch("istatr",&ntuple_st.istatr[0],"istatr[500]/I");
1894    ntuple_global->Branch("isignr",&ntuple_st.isignr[0],"isignr[500]/I");
1895    ntuple_global->Branch("pxr",&ntuple_st.pxr[0],"pxr[500]/F");
1896    ntuple_global->Branch("pyr",&ntuple_st.pyr[0],"pyr[500]/F");
1897    ntuple_global->Branch("pzr",&ntuple_st.pzr[0],"pzr[500]/F");
1898    ntuple_global->Branch("zvr",&ntuple_st.zvr[0],"zvr[500]/F");
1899    ntuple_global->Branch("chi2r",&ntuple_st.chi2r[0],"chi2r[500]/F");
1900    ntuple_global->Branch("pxv",&ntuple_st.pxv[0],"pxv[500]/F");
1901    ntuple_global->Branch("pyv",&ntuple_st.pyv[0],"pyv[500]/F");
1902    ntuple_global->Branch("pzv",&ntuple_st.pzv[0],"pzv[500]/F");
1903
1904    // test aliroot
1905
1906   new TH1F("h100","particule id du hit geant",20,0.,20.);
1907   new TH1F("h101","position en x du hit geant",100,-200.,200.);
1908   new TH1F("h102","position en y du hit geant",100,-200.,200.);
1909   new TH1F("h103","chambre de tracking concernee",15,0.,14.);
1910   new TH1F("h104","moment ptot du hit geant",50,0.,100.);
1911   new TH1F("h105","px au vertex",50,0.,20.);
1912   new TH1F("h106","py au vertex",50,0.,20.);
1913   new TH1F("h107","pz au vertex",50,0.,20.);
1914   new TH1F("h108","position zv",50,-15.,15.);
1915   new TH1F("h109","position en x du hit reconstruit",100,-300.,300.);
1916   new TH1F("h110","position en y du hit reconstruit",100,-300.,300.);
1917   new TH1F("h111","delta x ",100,-0.4,0.4);
1918   new TH1F("h112","delta y ",100,-0.4,0.4);
1919
1920   char hname[30];
1921   char hname1[30];
1922   for (int i=0;i<10;i++) {
1923     sprintf(hname,"deltax%d",i);
1924     sprintf(hname1,"h12%d",i);
1925     new TH1F(hname1,hname ,100,-0.4,0.4);
1926     sprintf(hname,"deltay%d",i);
1927     sprintf(hname1,"h13%d",i);
1928     new TH1F(hname1,hname ,100,-0.4,0.4);
1929   }
1930   new TH2F("h2000","VAR X st. 5",30,3.0,183.0,100,0.,25.);
1931   new TH2F("h2001","VAR Y st. 5",30,3.0,183.0,100,0.,25.);
1932
1933   new TH2F("h2500","P vs X HHIT",30,3.0,183.0,200,0.,200.);
1934   new TH2F("h2501","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
1935   new TH2F("h2502","P vs X EPH2 st. 5",30,3.0,183.0,100,0.,0.000005);
1936   new TH2F("h2503","P vs X EAL2 st. 5",30,3.0,183.0,100,0.,0.01);
1937   //new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,1.5);
1938   new TH2F("h2504","P vs X EXM2 st. 5",30,3.0,183.0,100,0.,0.1);
1939   new TH2F("h2505","P vs X EYM2 st. 5",30,3.0,183.0,100,0.,30.);
1940
1941   new TH2F("h2507","P vs X EPH st. 5",30,3.0,183.0,100,0.,0.003);
1942   new TH2F("h2508","P vs X EAL st. 5",30,3.0,183.0,100,0.,0.3);
1943   //new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,1.5);
1944   new TH2F("h2509","P vs X EXM st. 5",30,3.0,183.0,100,0.,0.4);
1945   new TH2F("h2510","P vs X EYM st. 5",30,3.0,183.0,100,0.,30.);
1946
1947   new TH2F("h2511","P vs X EPH cut st. 5",30,3.0,183.0,100,0.,0.01);
1948   new TH2F("h2512","P vs X EAL cut st. 5",30,3.0,183.0,100,0.,0.3);
1949   //new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,1.5);
1950   new TH2F("h2513","P vs X EXM cut st. 5",30,3.0,183.0,100,0.,0.4);
1951   new TH2F("h2514","P vs X EYM cut st. 5",30,3.0,183.0,100,0.,30.);
1952   // 4
1953   new TH2F("h2400","P vs X HHIT",30,3.0,183.0,200,0.,200.);
1954   new TH2F("h2401","P vs X HHIT**2",30,3.0,183.0,200,0.,5000.);
1955   new TH2F("h2402","P vs X EPH2 st. 4",30,3.0,183.0,100,0.,0.000005);
1956   new TH2F("h2403","P vs X EAL2 st. 4",30,3.0,183.0,100,0.,0.05);
1957   //new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,1.5);
1958   new TH2F("h2404","P vs X EXM2 st. 4",30,3.0,183.0,100,0.,0.1);
1959   new TH2F("h2405","P vs X EYM2 st. 4",30,3.0,183.0,100,0.,30.);
1960
1961   new TH2F("h2407","P vs X EPH st. 4",30,3.0,183.0,100,0.,0.003);
1962   new TH2F("h2408","P vs X EAL st. 4",30,3.0,183.0,100,0.,0.3);
1963   //new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,1.5);
1964   new TH2F("h2409","P vs X EXM st. 4",30,3.0,183.0,100,0.,0.1);
1965   new TH2F("h2410","P vs X EYM st. 4",30,3.0,183.0,100,0.,30.);
1966
1967   new TH2F("h2411","P vs X EPH cut st. 4",30,3.0,183.0,100,0.,0.01);
1968   new TH2F("h2412","P vs X EAL cut st. 4",30,3.0,183.0,100,0.,0.3);
1969   //new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,1.5);
1970   new TH2F("h2413","P vs X EXM cut st. 4",30,3.0,183.0,100,0.,0.1);
1971   new TH2F("h2414","P vs X EYM cut st. 4",30,3.0,183.0,100,0.,30.);
1972   // 3
1973   new TH1F("h2301","P2",30,3.0,183.0);
1974   new TH2F("h2302","P2 vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
1975   new TH2F("h2303","P2 vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.0005);
1976   //new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
1977   new TH2F("h2304","P2 vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
1978   new TH2F("h2305","P2 vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
1979
1980   new TH2F("h2307","P vs X EPH2 st. 3",30,3.0,183.0,100,0.,0.0006);
1981   new TH2F("h2308","P vs X EAL2 st. 3",30,3.0,183.0,100,0.,0.005);
1982   //new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,1.5);
1983   new TH2F("h2309","P vs X EXM2 st. 3",30,3.0,183.0,100,0.,2.);
1984   new TH2F("h2310","P vs X EYM2 st. 3",30,3.0,183.0,100,0.,3.);
1985
1986   new TH2F("h2311","P vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
1987   new TH2F("h2312","P vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
1988   //new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
1989   new TH2F("h2313","P vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
1990   new TH2F("h2314","P vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
1991
1992   new TH2F("h2315","P2 vs X EPH cut st. 3",30,3.0,183.0,100,0.,0.06);
1993   new TH2F("h2316","P2 vs X EAL cut st. 3",30,3.0,183.0,100,0.,0.05);
1994   //new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,1.5);
1995   new TH2F("h2317","P2 vs X EXM cut st. 3",30,3.0,183.0,100,0.,6.);
1996   new TH2F("h2318","P2 vs X EYM cut st. 3",30,3.0,183.0,100,0.,7.);
1997   
1998   // 2
1999   new TH1F("h2201","P2",30,3.0,183.0);
2000   new TH2F("h2202","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2001   new TH2F("h2203","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2002   //new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2003   new TH2F("h2204","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2004   new TH2F("h2205","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
2005
2006   new TH2F("h2207","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2007   new TH2F("h2208","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2008   //new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2009   new TH2F("h2209","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2010   new TH2F("h2210","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,5.);
2011
2012   new TH2F("h2211","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
2013   new TH2F("h2212","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2014   //new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2015   new TH2F("h2213","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2016   new TH2F("h2214","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
2017
2018   new TH2F("h2215","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.05);
2019   new TH2F("h2216","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2020   //new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2021   new TH2F("h2217","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2022   new TH2F("h2218","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,10.);
2023
2024   // 1
2025   new TH2F("h2102","P2 vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2026   new TH2F("h2103","P2 vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2027   //new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2028   new TH2F("h2104","P2 vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2029   new TH2F("h2105","P2 vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
2030
2031   new TH2F("h2107","P vs X EPH2 st. 2",30,3.0,183.0,100,0.,0.0006);
2032   new TH2F("h2108","P vs X EAL2 st. 2",30,3.0,183.0,100,0.,0.005);
2033   //new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,1.5);
2034   new TH2F("h2109","P vs X EXM2 st. 2",30,3.0,183.0,100,0.,7.);
2035   new TH2F("h2110","P vs X EYM2 st. 2",30,3.0,183.0,100,0.,7.);
2036
2037   new TH2F("h2111","P vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
2038   new TH2F("h2112","P vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2039   //new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2040   new TH2F("h2113","P vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2041   new TH2F("h2114","P vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
2042
2043   new TH2F("h2115","P2 vs X EPH cut st. 2",30,3.0,183.0,100,0.,0.1);
2044   new TH2F("h2116","P2 vs X EAL cut st. 2",30,3.0,183.0,100,0.,0.2);
2045   //new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,1.5);
2046   new TH2F("h2117","P2 vs X EXM cut st. 2",30,3.0,183.0,100,0.,11.);
2047   new TH2F("h2118","P2 vs X EYM cut st. 2",30,3.0,183.0,100,0.,11.);
2048
2049   // 2,3,4,5
2050   new TH1F("h2701","P2 fit 2",30,3.0,183.0);
2051   new TH2F("h2702","P2 vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
2052   new TH2F("h2703","P2 vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
2053   // new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2054   new TH2F("h2704","P2 vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
2055   new TH2F("h2705","P2 vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
2056
2057   new TH2F("h2707","P vs X EPH2 st. 1 fit 2",30,3.0,183.0,100,0.,0.0006);
2058   new TH2F("h2708","P vs X EAL2 st. 1 fit 2",30,3.0,183.0,100,0.,0.005);
2059   //new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2060   new TH2F("h2709","P vs X EXM2 st. 1 fit 2",30,3.0,183.0,100,0.,2.);
2061   new TH2F("h2710","P vs X EYM2 st. 1 fit 2",30,3.0,183.0,100,0.,3.);
2062
2063   new TH2F("h2711","P vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
2064   new TH2F("h2712","P vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
2065   //new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2066   new TH2F("h2713","P vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
2067   new TH2F("h2714","P vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
2068
2069   new TH2F("h2715","P2 vs X EPH cut st. 1 fit 2",30,3.0,183.0,100,0.,0.07);
2070   new TH2F("h2716","P2 vs X EAL cut st. 1 fit 2",30,3.0,183.0,100,0.,0.2);
2071   //new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,1.5);
2072   new TH2F("h2717","P2 vs X EXM cut st. 1 fit 2",30,3.0,183.0,100,0.,6.);
2073   new TH2F("h2718","P2 vs X EYM cut st. 1 fit 2",30,3.0,183.0,100,0.,7.);
2074
2075   // 1,3,4,5
2076   new TH1F("h2801","P2 fit 1",30,3.0,183.0);
2077   new TH2F("h2802","P2 vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
2078   new TH2F("h2803","P2 vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
2079   //new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2080   new TH2F("h2804","P2 vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
2081   new TH2F("h2805","P2 vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
2082
2083   new TH2F("h2807","P vs X EPH2 st. 2 fit 1",30,3.0,183.0,100,0.,0.0006);
2084   new TH2F("h2808","P vs X EAL2 st. 2 fit 1",30,3.0,183.0,100,0.,0.005);
2085   //new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2086   new TH2F("h2809","P vs X EXM2 st. 2 fit 1",30,3.0,183.0,100,0.,2.);
2087   new TH2F("h2810","P vs X EYM2 st. 2 fit 1",30,3.0,183.0,100,0.,3.);
2088
2089   new TH2F("h2811","P vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
2090   new TH2F("h2812","P vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
2091   //new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2092   new TH2F("h2813","P vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
2093   new TH2F("h2814","P vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
2094
2095   new TH2F("h2815","P2 vs X EPH cut st. 2 fit 1",30,3.0,183.0,100,0.,0.05);
2096   new TH2F("h2816","P2 vs X EAL cut st. 2 fit 1",30,3.0,183.0,100,0.,0.2);
2097   //new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,1.5);
2098   new TH2F("h2817","P2 vs X EXM cut st. 2 fit 1",30,3.0,183.0,100,0.,5.);
2099   new TH2F("h2818","P2 vs X EYM cut st. 2 fit 1",30,3.0,183.0,100,0.,7.);
2100   // fin de test
2101
2102   new TH1F("h500","Acceptance en H st. 4",500,0.,500.);
2103   new TH1F("h600","Acceptance en H st. 5",500,0.,500.);
2104   new TH1F("h700","X vertex track found",200,-10.,10.);
2105   new TH1F("h701","Y vertex track found",200,-10.,10.);
2106   new TH1F("h800","Rap. muon gen.",100,0.,5.);
2107   new TH1F("h801","Rap. muon gen. recons.",100,0.,5.);
2108   new TH1F("h802","Rap. muon gen. ghost ",100,0.,5.);
2109   new TH1F("h900","Pt muon gen.",100,0.,20.);
2110   new TH1F("h901","Pt muon gen. recons.",100,0.,20.);
2111   new TH1F("h902","Pt muon gen. ghost",100,0.,20.);
2112   new TH1F("h910","phi muon gen.",100,-10.,10.);
2113   new TH1F("h911","phi muon gen. recons.",100,-10.,10.);
2114   new TH1F("h912","phi muon gen. ghost",100,-10.,10.);
2115   new TH2F("h1001","Y VS X hit st. 1",300,-300.,300.,300,-300.,300.);
2116   new TH2F("h1002","Y VS X hit st. 2",300,-300.,300.,300,-300.,300.);
2117   new TH2F("h1003","Y VS X hit st. 3",300,-300.,300.,300,-300.,300.);
2118   new TH2F("h1004","Y VS X hit st. 4",300,-300.,300.,300,-300.,300.);
2119   new TH2F("h1005","Y VS X hit st. 5",300,-300.,300.,300,-300.,300.);
2120   //  Histos variance dans 4      
2121   new TH2F("h11","VAR X st. 4",30,3.0,183.0,100,0.,2.);
2122   new TH2F("h12","VAR Y st. 4",30,3.0,183.0,100,0.,600.);
2123   new TH2F("h13","VAR PHI st. 4",30,3.0,183.0,100,0.,0.0001);
2124   new TH2F("h14","VAR ALM st. 4",30,3.0,183.0,100,0.,0.05);
2125   new TH1F("h15","P",30,3.0,183.0);
2126   new TH1F("h411","VAR X st. 4",100,-1.42,1.42);
2127   new TH1F("h412","VAR Y st. 4",100,-25.,25.);
2128   new TH1F("h413","VAR PHI st. 4",100,-0.01,0.01);
2129   new TH1F("h414","VAR ALM st. 4",100,-0.23,0.23);
2130   // histo2
2131   new TH2F("h211","histo2-VAR X st. 4",30,3.0,183.0,100,0.,2.);
2132   new TH2F("h212","histo2-VAR Y st. 4",30,3.0,183.0,100,0.,600.);
2133   new TH1F("h213","histo2-VAR X st. 4",100,-1.42,1.42);
2134   new TH1F("h214","histo2-VAR Y st. 4",100,-25.,25.);
2135   new TH1F("h215","histo2-P",30,3.0,183.0);
2136
2137   //  Histos variance dans 2      
2138   new TH2F("h21","VAR X st. 2",30,3.0,183.0,100,0.,3.);
2139   new TH2F("h22","VAR Y st. 2",30,3.0,183.0,100,0.,7.);
2140   new TH2F("h23","VAR PHI st. 2",30,3.0,183.0,100,0.,0.006);
2141   new TH2F("h24","VAR ALM st. 2",30,3.0,183.0,100,0.,0.005);
2142   new TH1F("h25","P",30,3.0,183.0);
2143   new TH1F("h421","VAR X st. 2",100,-1.72,1.72);
2144   new TH1F("h422","VAR Y st. 2",100,-2.7,2.7);
2145   new TH1F("h423","VAR PHI st. 2",100,-0.08,0.08);
2146   new TH1F("h424","VAR ALM st. 2",100,-0.072,0.072);
2147   // histo2
2148   new TH2F("h221","histo2-VAR X st. 2",30,3.0,183.0,100,0.,3.);
2149   new TH2F("h222","histo2-VAR Y st. 2",30,3.0,183.0,100,0.,7.);
2150   new TH1F("h223","histo2-VAR X st. 2",100,-1.72,1.72);
2151   new TH1F("h224","histo2-VAR Y st. 2",100,-2.7,2.7);
2152   new TH1F("h225","histo2-P",30,3.0,183.0);
2153
2154   //  Histos variance dans 1      
2155   new TH2F("h31","VAR X st. 1",30,3.0,183.0,100,0.,2.);
2156   new TH2F("h32","VAR Y st. 1",30,3.0,183.0,100,0.,0.5);
2157   new TH2F("h33","VAR PHI st. 1",30,3.0,183.0,100,0.,0.006);
2158   new TH2F("h34","VAR ALM st. 1",30,3.0,183.0,100,0.,0.005);
2159   new TH1F("h35","P",30,3.0,183.0);
2160   new TH1F("h431","VAR X st. 1",100,-1.42,1.42);
2161   new TH1F("h432","VAR Y st. 1",100,-0.72,0.72);
2162   new TH1F("h433","VAR PHI st. 1",100,-0.08,0.08);
2163   new TH1F("h434","VAR ALM st. 1",100,-0.072,0.072);
2164   //  Histos variance dans 1      
2165   new TH2F("h41","VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
2166   new TH2F("h42","VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
2167   new TH2F("h43","VAR PHI st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
2168   new TH2F("h44","VAR ALM st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,0.005);
2169   new TH1F("h45","P",30,3.0,183.0);
2170   new TH1F("h441","VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
2171   new TH1F("h442","VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
2172   new TH1F("h443","VAR PHI st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
2173   new TH1F("h444","VAR ALM st. 1 fit 5,4,3,2,V",100,-0.072,0.072);
2174   // histo2
2175   new TH2F("h241","histo2-VAR X st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,4.);
2176   new TH2F("h242","histo2-VAR Y st. 1 fit 5,4,3,2,V",30,3.0,183.0,100,0.,20.);
2177   new TH1F("h243","histo2-VAR X st. 1 fit 5,4,3,2,V",100,-2.,2.);
2178   new TH1F("h244","histo2-VAR Y st. 1 fit 5,4,3,2,V",100,-4.5,4.5);
2179   new TH1F("h245","histo2-P",30,3.0,183.0);
2180
2181   //  Histos variance dans 2      
2182   new TH2F("h51","VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
2183   new TH2F("h52","VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
2184   new TH2F("h53","VAR PHI st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.005);
2185   new TH2F("h54","VAR ALM st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.01);
2186   new TH1F("h55","P",30,3.0,183.0);
2187   new TH1F("h451","VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
2188   new TH1F("h452","VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
2189   new TH1F("h453","VAR PHI st. 2 fit 5,4,3,1,V",100,-0.072,0.072);
2190   new TH1F("h454","VAR ALM st. 2 fit 5,4,3,1,V",100,-0.1,0.1);
2191   new TH1F("h999","PTOT",30,3.0,183.0);
2192   // histo2
2193   new TH2F("h251","histo2-VAR X st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,0.5);
2194   new TH2F("h252","histo2-VAR Y st. 2 fit 5,4,3,1,V",30,3.0,183.0,100,0.,2.);
2195   new TH1F("h253","histo2-VAR X st. 2 fit 5,4,3,1,V",100,-0.72,0.72);
2196   new TH1F("h254","histo2-VAR Y st. 2 fit 5,4,3,1,V",100,-1.42,1.42);
2197   new TH1F("h255","histo2-P",30,3.0,183.0);
2198   //  Histos variance dans 3      
2199   new TH2F("h61","VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
2200   new TH2F("h62","VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
2201   new TH2F("h63","VAR PHI st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
2202   new TH2F("h64","VAR ALM st. 3 fit 4,5,V",30,3.0,183.0,100,0.,0.0006);
2203   new TH1F("h65","P",30,3.0,183.0);
2204   new TH1F("h461","VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
2205   new TH1F("h462","VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
2206   new TH1F("h463","VAR PHI st. 3 fit 4,5,V",100,-0.024,0.024);
2207   new TH1F("h464","VAR ALM st. 3 fit 4,5,V",100,-0.024,0.024);
2208   // histo2
2209   new TH2F("h261","histo2-VAR X st. 3 fit 4,5,V",30,3.0,183.0,100,0.,5.);
2210   new TH2F("h262","histo2-VAR Y st. 3 fit 4,5,V",30,3.0,183.0,100,0.,2.);
2211   new TH1F("h263","histo2-VAR X st. 3 fit 4,5,V",100,-2.25,2.25);
2212   new TH1F("h264","histo2-VAR Y st. 3 fit 4,5,V",100,-1.42,1.42);
2213   new TH1F("h265","Phisto2-",30,3.0,183.0);
2214   // Histos dx,dy distribution between chambers inside stations
2215   new TH1F("h71","DX in st. ID-70",100,-5.,5.);
2216   new TH1F("h81","DY in st. ID-80",100,-5.,5.);
2217   new TH1F("h72","DX in st. ID-70",100,-5.,5.);
2218   new TH1F("h82","DY in st. ID-80",100,-5.,5.);
2219   new TH1F("h73","DX in st. ID-70",100,-5.,5.);
2220   new TH1F("h83","DY in st. ID-80",100,-5.,5.);
2221   new TH1F("h74","DX in st. ID-70",100,-5.,5.);
2222   new TH1F("h84","DY in st. ID-80",100,-5.,5.);
2223   new TH1F("h75","DX in st. ID-70",100,-5.,5.);
2224   new TH1F("h85","DY in st. ID-80",100,-5.,5.);
2225 }
2226
2227 void chfnt(Int_t &ievr, Int_t &ntrackr, Int_t *istatr, Int_t *isignr, Float_t *pxr, Float_t *pyr, Float_t *pzr, Float_t *zvr, Float_t *chi2r,  Float_t *pxv, Float_t *pyv, Float_t *pzv)
2228 {
2229   //
2230   // fill the ntuple 
2231     ntuple_st.ievr = ievr;
2232     ntuple_st.ntrackr = ntrackr;
2233     for (Int_t i=0; i<500; i++) {
2234         ntuple_st.istatr[i] = istatr[i];
2235         ntuple_st.isignr[i] = isignr[i]; 
2236         ntuple_st.pxr[i]    = pxr[i]; 
2237         ntuple_st.pyr[i]    = pyr[i];
2238         ntuple_st.pzr[i]    = pzr[i];
2239         ntuple_st.zvr[i]    = zvr[i];
2240         ntuple_st.chi2r[i]  = chi2r[i];
2241         ntuple_st.pxv[i]    = pxv[i]; 
2242         ntuple_st.pyv[i]    = pyv[i];
2243         ntuple_st.pzv[i]    = pzv[i];
2244     }
2245     ntuple_global->Fill();   
2246 }
2247
2248 void hist_closed()
2249 {
2250   //
2251   // write histos and ntuple to "reconst.root" file
2252   hfile_global->Write();
2253 }
2254
2255 void trackf_read_geant(Int_t *itypg, Double_t *xtrg, Double_t *ytrg, Double_t *ptotg, Int_t *idg, Int_t *izch, Double_t *pvert1g, Double_t *pvert2g, Double_t *pvert3g, Double_t *zvertg, Int_t &nhittot1, Double_t *cx, Double_t *cy, Double_t *cz, Int_t &ievr,Int_t &nev,Double_t *xgeant, Double_t *ygeant,Double_t *clsize1, Double_t *clsize2) 
2256 {
2257   //
2258   // introduce aliroot variables in fortran common 
2259   // tracking study from geant hits 
2260   //
2261
2262   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
2263   
2264   //  TTree *TK = gAlice->TreeK();
2265   TTree *TH = gAlice->TreeH();
2266   Int_t ntracks = (Int_t)TH->GetEntries();
2267   cout<<"ntrack="<<ntracks<<endl;
2268
2269   Int_t maxidg = 0;
2270   Int_t nres=0;
2271   
2272 //
2273 //  Loop over tracks
2274 //
2275
2276   for (Int_t track=0; track<ntracks;track++) {
2277       gAlice->ResetHits();
2278       TH->GetEvent(track);
2279       
2280       if (MUON)  {
2281 //
2282 //  Loop over hits
2283 //
2284           for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); 
2285               mHit;
2286               mHit=(AliMUONhit*)MUON->NextHit()) 
2287           {
2288               if (maxidg<=20000) {
2289                 
2290                 if (mHit->fChamber > 10) continue;
2291                 TClonesArray *fPartArray = gAlice->Particles();
2292                 TParticle *Part;
2293                 Int_t ftrack = mHit->fTrack;
2294                 Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
2295
2296                 if (id==kMuonPlus||id==kMuonMinus) {
2297                     
2298                     // inversion de x et y car le champ est inverse dans le programme tracking
2299                     xtrg[maxidg]   = 0;       
2300                     ytrg[maxidg]   = 0;       
2301                     xgeant[maxidg]   = mHit->fY;             // x-pos of hit
2302                     ygeant[maxidg]   = mHit->fX;             // y-pos of hit
2303                     clsize1[maxidg]   = 0;     // cluster size on 1-st cathode
2304                     clsize2[maxidg]   = 0;     // cluster size on 2-nd cathode
2305                     cx[maxidg]     = mHit->fCyHit;            // Px/P of hit
2306                     cy[maxidg]     = mHit->fCxHit;            // Py/P of hit
2307                     cz[maxidg]     = mHit->fCzHit;            // Pz/P of hit
2308                     izch[maxidg]   = mHit->fChamber; 
2309                     /*      
2310                     Int_t pdgtype  = Int_t(mHit->fParticle); // particle number
2311                     itypg[maxidg]  = gMC->IdFromPDG(pdgtype);
2312
2313                     */
2314                     if (id==kMuonPlus) itypg[maxidg]  = 5;
2315                     else  itypg[maxidg]  = 6;
2316
2317                     ptotg[maxidg]  = mHit->fPTot;          // P of hit 
2318                     
2319                     Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
2320                     Float_t thet = Part->Theta();
2321                     thet = thet*180./3.1416;
2322                     
2323                     Int_t iparent = Part->GetFirstMother();
2324                     if (iparent >= 0) {
2325                         Int_t ip;
2326                         while(1) {
2327                             ip=((TParticle*) fPartArray->UncheckedAt(iparent))->GetFirstMother();
2328                             if (ip < 0) {
2329                                 break;
2330                             } else {
2331                                 iparent = ip;
2332                             }
2333                         }
2334                     }
2335                     //printf("iparent - %d\n",iparent);
2336                     Int_t id1  = ftrack; // numero de la particule generee au vertex
2337                     Int_t idum = track+1;
2338                     Int_t id2 = ((TParticle*) fPartArray->UncheckedAt(iparent))->GetPdgCode();
2339
2340                     if (id2==443) id2=114;
2341                     else id2=116;
2342                    
2343                     if (id2==116) {
2344                       nres++;
2345                     }
2346                     //printf("id2 %d\n",id2);
2347                     idg[maxidg] = 30000*id1+10000*idum+id2;
2348                     
2349                     pvert1g[maxidg] = Part->Py();      // Px vertex
2350                     pvert2g[maxidg] = Part->Px();      // Py vertex  
2351                     pvert3g[maxidg] = Part->Pz();      // Pz vertex
2352                     zvertg[maxidg]  = Part->Vz();      // z vertex 
2353                     maxidg ++;
2354
2355                 }
2356               }
2357           } // hit loop
2358       } // if MUON
2359   } // track loop first file
2360
2361   if (TrH1 && fHits2 ) { // if background file
2362     ntracks =(Int_t)TrH1->GetEntries();
2363     printf("Trackf_read - 2-nd file - ntracks %d\n",ntracks);
2364
2365     //  Loop over tracks
2366     for (Int_t track=0; track<ntracks; track++) {
2367       
2368       if (fHits2) fHits2->Clear();
2369       TrH1->GetEvent(track);
2370
2371       //  Loop over hits
2372       for (int i=0;i<fHits2->GetEntriesFast();i++) 
2373         {
2374           AliMUONhit *mHit=(AliMUONhit*) (*fHits2)[i];
2375          
2376           if (mHit->fChamber > 10) continue;
2377
2378           if (maxidg<=20000) {
2379             
2380             // inversion de x et y car le champ est inverse dans le programme tracking !!!!
2381             xtrg[maxidg]   = 0;                    // only for reconstructed point
2382             ytrg[maxidg]   = 0;                    // only for reconstructed point
2383             xgeant[maxidg]   = mHit->fY;           // x-pos of hit
2384             ygeant[maxidg]   = mHit->fX;           // y-pos of hit
2385             clsize1[maxidg]   = 0;           // cluster size on 1-st cathode
2386             clsize2[maxidg]   = 0;           // cluster size on 2-nd cathode
2387             cx[maxidg]     = mHit->fCyHit;         // Px/P of hit
2388             cy[maxidg]     = mHit->fCxHit;         // Py/P of hit
2389             cz[maxidg]     = mHit->fCzHit;         // Pz/P of hit
2390             izch[maxidg]   = mHit->fChamber;       // chamber number
2391             ptotg[maxidg]  = mHit->fPTot;          // P of hit 
2392             
2393             Int_t ftrack = mHit->fTrack;
2394             Int_t id1  = ftrack;                   // track number 
2395             Int_t idum = track+1;
2396             
2397             TClonesArray *fPartArray = fParticles2;
2398             TParticle *Part;
2399             Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
2400             Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
2401             if (id==kMuonPlus||id==kMuonMinus) {
2402                 if (id==kMuonPlus) itypg[maxidg]  = 5;
2403                 else  itypg[maxidg]  = 6;
2404             } else itypg[maxidg]=0;
2405             
2406             Int_t id2=0; // set parent to 0 for background !!
2407             idg[maxidg] = 30000*id1+10000*idum+id2;
2408             
2409             pvert1g[maxidg] = Part->Py();      // Px vertex
2410             pvert2g[maxidg] = Part->Px();      // Py vertex  
2411             pvert3g[maxidg] = Part->Pz();      // Pz vertex
2412             zvertg[maxidg]  = Part->Vz();      // z vertex 
2413
2414             maxidg ++;
2415
2416           } // check limits (maxidg)
2417         } // hit loop 
2418     } // track loop
2419   } // if TrH1
2420
2421   ievr = nev;
2422   nhittot1 = maxidg ;
2423   cout<<"nhittot1="<<nhittot1<<endl;
2424
2425   static Int_t nbres=0;
2426   if (nres>=19) nbres++;
2427   printf("nres, nbres %d %d \n",nres,nbres);
2428
2429   hfile_global->cd();      
2430
2431 }
2432
2433
2434
2435 void trackf_read_spoint(Int_t *itypg, Double_t *xtrg, Double_t *ytrg, Double_t *ptotg, Int_t *idg, Int_t *izch, Double_t *pvert1g, Double_t *pvert2g, Double_t *pvert3g, Double_t *zvertg, Int_t &nhittot1, Double_t *cx, Double_t *cy, Double_t *cz, Int_t &ievr,Int_t &nev,Double_t *xgeant, Double_t *ygeant,Double_t *clsize1, Double_t *clsize2) 
2436
2437 {
2438   //
2439   // introduce aliroot variables in fortran common 
2440   // tracking study from reconstructed points 
2441   //
2442   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
2443
2444   cout<<"numero de l'evenement "<<nev<<endl;
2445   
2446   MUON->GetTreeC(nev);
2447   TTree *TC=MUON->TreeC();
2448   TC->GetEntries();
2449
2450   Int_t maxidg = 0;
2451   Int_t nres=0;
2452   Int_t nncor=0;
2453   static Int_t nuncor=0;
2454   static Int_t nbadcor=0;
2455   AliMUONRawCluster * mRaw;
2456   AliMUONRawCluster * mRaw1;
2457   TTree *TH = gAlice->TreeH();
2458
2459   Int_t ihit;
2460   Int_t mult1, mult2;
2461   if (MUON) {
2462       for (Int_t ich=0;ich<10;ich++) {
2463           TClonesArray *MUONcorrel  = MUON->CathCorrelAddress(ich);
2464           MUON->ResetCorrelation();
2465           TC->GetEvent();
2466           Int_t ncor = (Int_t)MUONcorrel->GetEntries();
2467           if (ncor>=2) nncor++;
2468           if (!ncor) continue;
2469
2470           //  Loop over correlated clusters
2471           for (Int_t icor=0;icor<ncor;icor++) {
2472               AliMUONcorrelation * mCor = (AliMUONcorrelation*)MUONcorrel->UncheckedAt(icor);
2473
2474               Int_t flag=0;   // = 1 if no information in the second cathode
2475               Int_t index = mCor->fCorrelIndex[0]; // for the second cathode
2476               if (index >= 0) {
2477                   Int_t index1 = mCor->fCorrelIndex[3]; // for the 1-st cathode
2478                   mRaw1 = MUON->RawCluster(ich,1,index1);
2479                   mult1=mRaw1->fMultiplicity;
2480                   mRaw = MUON->RawCluster(ich,2,index);
2481                   mult2=mRaw->fMultiplicity;
2482               } else {
2483                   index = mCor->fCorrelIndex[3];
2484                   mRaw = MUON->RawCluster(ich,1,index);
2485                   mult1=mRaw->fMultiplicity;
2486                   mult2=0;
2487                   flag=1;
2488                   nuncor++;
2489               }
2490               if (!mRaw) continue;
2491
2492               Int_t ftrack1 = mRaw->fTracks[1]; // qui doit etre le meme pour 
2493                                                 // la cathode 1 et 2
2494               ihit= mRaw->fTracks[0];
2495               //printf("icor, ftrack1 ihit %d %d %d\n",icor,ftrack1,ihit);
2496
2497               if (mRaw->fClusterType == 0 ) {
2498
2499                   if (maxidg<=20000) {
2500                       if (flag == 0) {
2501                           xtrg[maxidg]   = (Double_t) mCor->fY[3];
2502                           ytrg[maxidg]   = (Double_t) mCor->fX[0];
2503                           Int_t index1 = mCor->fCorrelIndex[3];
2504                           mRaw1 = MUON->RawCluster(ich,1,index1);
2505                           if (mRaw1->fClusterType==1 || mRaw1->fClusterType==2) {
2506                             Float_t xclust=mCor->fX[3];
2507                             Float_t yclust=mCor->fY[3];
2508                             AliMUONchamber *iChamber=&(MUON->Chamber(ich));
2509                             AliMUONsegmentation *seg = iChamber->GetSegmentationModel(1);
2510                             Int_t ix,iy;
2511                             seg->GetPadIxy(xclust,yclust,ix,iy);   
2512                             Int_t isec=seg->Sector(ix,iy);
2513                             printf("nev, CORRELATION with pure background in chamber sector %d  %d  %d !!!!!!!!!!!!!!!!!!!!!\n",nev,ich+1,isec);
2514                             nbadcor++;
2515                             
2516                           } // end if cluster type on cathode 1
2517                       }else {
2518                           xtrg[maxidg]   = (Double_t) mCor->fY[3];
2519                           ytrg[maxidg]   = (Double_t) mCor->fX[3];
2520                       } // if iflag
2521                       izch[maxidg]   = ich+1;
2522                       xgeant[maxidg] = 0;
2523                       ygeant[maxidg] = 0;
2524                       clsize1[maxidg] = mult1;
2525                       clsize2[maxidg] = mult2;
2526
2527                       cx[maxidg]     = 0; // Px/P of hit
2528                       cy[maxidg]     = 0; // Py/P of hit
2529                       cz[maxidg]     = 0; // Pz/P of hit
2530                       itypg[maxidg]  = 0; // particle number
2531                       ptotg[maxidg]  = 0; // P of hit
2532                       idg[maxidg]    = 0;
2533                       pvert1g[maxidg] = 0; // Px vertex
2534                       pvert2g[maxidg] = 0; // Py vertex  
2535                       pvert3g[maxidg] = 0; // Pz vertex
2536                       zvertg[maxidg]  = 0; // z vertex     
2537                       maxidg++;
2538                       
2539                   }// fin maxidg
2540                   
2541               } else if (mRaw->fClusterType ==1 && ftrack1 < 0) // background + resonance
2542                 {
2543                   nres++;
2544                   // get indexmap and loop over digits to find the signal
2545                   Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
2546                   gAlice->ResetDigits();
2547                   if (flag==0) {
2548                     //gAlice->TreeD()->GetEvent(2); // cathode 2
2549                     gAlice->TreeD()->GetEvent(nent-1); // cathode 2
2550                   } else {
2551                     //gAlice->TreeD()->GetEvent(1); // cathode 1
2552                     gAlice->TreeD()->GetEvent(nent-2); // cathode 1
2553                   }
2554
2555                    TClonesArray *MUONdigits  = MUON->DigitsAddress(ich);
2556                    Int_t mul=mRaw->fMultiplicity;
2557                    Int_t trsign;
2558                    for (int i=0;i<mul;i++) {
2559                      Int_t idx=mRaw->fIndexMap[i];
2560                      AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
2561                      trsign=dig->fTracks[0];
2562                      ihit=dig->fHit-1;
2563                      if (trsign > 0 && ihit >= 0) break;
2564
2565                    } // loop over indexmap
2566
2567                    //printf("trsign, ihit %d %d\n",trsign,ihit);
2568                    //printf("signal+background : trsign  %d\n",trsign);
2569                   
2570                    if (trsign < 0 || ihit < 0) { // no signal muon  was found
2571                      
2572                      if (maxidg<=20000) {
2573                        if (flag == 0) {
2574                          xtrg[maxidg]   = (Double_t) mCor->fY[3];
2575                          ytrg[maxidg]   = (Double_t) mCor->fX[0];
2576                        }else {
2577                          xtrg[maxidg]   = (Double_t) mCor->fY[3];
2578                          ytrg[maxidg]   = (Double_t) mCor->fX[3];
2579                        }
2580                        
2581                        izch[maxidg]   = ich+1;
2582
2583                       // initialisation of informations which 
2584                       // can't be reached for background
2585                        
2586                        xgeant[maxidg] = 0; // only for resonances
2587                        ygeant[maxidg] = 0; // only for resonances
2588                        clsize1[maxidg] = mult1;
2589                        clsize2[maxidg] = mult2;
2590
2591                        cx[maxidg]     = 0; // Px/P of hit
2592                        cy[maxidg]     = 0; // Py/P of hit
2593                        cz[maxidg]     = 0; // Pz/P of hit
2594                        itypg[maxidg]  = 0; // particle number
2595                        ptotg[maxidg]  = 0; // P of hit
2596                        idg[maxidg]    = 0;
2597                        pvert1g[maxidg] = 0; // Px vertex
2598                        pvert2g[maxidg] = 0; // Py vertex  
2599                        pvert3g[maxidg] = 0; // Pz vertex
2600                        zvertg[maxidg]  = 0;                
2601                        maxidg++;
2602                        
2603                      }// fin maxidg
2604                    } else { // signal muon - retrieve info
2605                      //printf("inside trsign, ihit %d %d\n",trsign,ihit);
2606                      if (maxidg<=20000) {
2607                        if (flag == 0) {
2608                          xtrg[maxidg]   = (Double_t) mCor->fY[3];
2609                          ytrg[maxidg]   = (Double_t) mCor->fX[0];
2610                        }else {
2611                          xtrg[maxidg]   = (Double_t) mCor->fY[3];
2612                          ytrg[maxidg]   = (Double_t) mCor->fX[3];
2613                        }
2614                        izch[maxidg]   = ich+1;
2615                        clsize1[maxidg] = mult1;
2616                        clsize2[maxidg] = mult2;
2617
2618                       // initialise and set to the correct values 
2619                       // if signal muons 
2620                        
2621                        xgeant[maxidg] = 0; // only for resonances
2622                        ygeant[maxidg] = 0; // only for resonances
2623                        
2624                        cx[maxidg]     = 0; // Px/P of hit
2625                        cy[maxidg]     = 0; // Py/P of hit
2626                        cz[maxidg]     = 0; // Pz/P of hit
2627                        itypg[maxidg]  = 0; // particle number
2628                        ptotg[maxidg]  = 0; // P of hit
2629                        idg[maxidg]    = 0;
2630                        pvert1g[maxidg] = 0; // Px vertex
2631                        pvert2g[maxidg] = 0; // Py vertex  
2632                        pvert3g[maxidg] = 0; // Pz vertex
2633                        zvertg[maxidg]  = 0;     
2634                        // try to retrieve info about signal muons          
2635                        gAlice->ResetHits();
2636                        TH->GetEvent(trsign);
2637
2638                        TClonesArray *MUONhits  = MUON->Hits();
2639                        AliMUONhit *mHit= (AliMUONhit*)MUONhits->
2640                                                         UncheckedAt(ihit);
2641                            TClonesArray *fPartArray = gAlice->Particles();
2642                            TParticle *Part;
2643                            Int_t nch=mHit->fChamber-1;
2644                            //printf("sig+bgr ich, nch %d %d \n",ich,nch);
2645                            if (nch==ich) {
2646                              Int_t ftrack = mHit->fTrack;
2647                              Int_t id = ((TParticle*) fPartArray->
2648                                         UncheckedAt(ftrack))->GetPdgCode();
2649                              if (id==kMuonPlus||id==kMuonMinus) {
2650                                  xgeant[maxidg] = (Double_t) mHit->fY;
2651                                  ygeant[maxidg] = (Double_t) mHit->fX;
2652                                  cx[maxidg]     = (Double_t) mHit->fCyHit; 
2653                                  cy[maxidg]     = (Double_t) mHit->fCxHit; 
2654                                  cz[maxidg]     = (Double_t) mHit->fCzHit; 
2655
2656                                  if (id==kMuonPlus) {
2657                                    itypg[maxidg]  = 5;
2658                                  } else if (id==kMuonMinus) {
2659                                    itypg[maxidg]  = 6;
2660                                  } else itypg[maxidg]  = 0;
2661                              
2662                                  ptotg[maxidg]  = (Double_t) mHit->fPTot;  
2663                                  Part = (TParticle*) fPartArray->
2664                                                      UncheckedAt(ftrack);
2665                                  Int_t iparent = Part->GetFirstMother();
2666                                  Int_t id2;
2667                                  id2 = ((TParticle*) fPartArray->
2668                                         UncheckedAt(ftrack))->GetPdgCode();
2669                              
2670                                  if (iparent >= 0) {
2671                                    Int_t ip;
2672                                    while(1) {
2673                                      ip=((TParticle*) fPartArray->
2674                                        UncheckedAt(iparent))->GetFirstMother();
2675                                      if (ip < 0) {
2676                                        id2 = ((TParticle*) fPartArray->
2677                                            UncheckedAt(iparent))->GetPdgCode();
2678                                        break;
2679                                      } else {
2680                                        iparent = ip;
2681                                        id2 = ((TParticle*) fPartArray->
2682                                            UncheckedAt(iparent))->GetPdgCode();
2683                                      } // ip<0
2684                                    } // while
2685                                  }// iparent
2686                                  Int_t id1  = ftrack; 
2687                                  Int_t idum = trsign+1;
2688                              
2689                                  if (id2==443 || id2==553) {
2690                                    nres++;
2691                                    if (id2==443) id2=114;
2692                                    else id2=116;
2693                                  }
2694                              
2695                                  idg[maxidg] = 30000*id1+10000*idum+id2;
2696                                  pvert1g[maxidg] = (Double_t) Part->Py(); 
2697                                  pvert2g[maxidg] = (Double_t) Part->Px();   
2698                                  pvert3g[maxidg] = (Double_t) Part->Pz(); 
2699                                  zvertg[maxidg]  = (Double_t) Part->Vz();  
2700                              } //if muon                             
2701                            } //if nch
2702                      maxidg++;
2703                      } // check limits
2704                    } // sign+bgr, highest bgr
2705               } 
2706               //pure resonance or mixed cluster with the highest 
2707               //contribution coming from resonance
2708               if (mRaw->fClusterType >= 1 && ftrack1>=0)  
2709                 {                             
2710                   if (maxidg<=20000) {
2711                     if (flag == 0) {
2712                       xtrg[maxidg]   = (Double_t) mCor->fY[3];
2713                       ytrg[maxidg]   = (Double_t) mCor->fX[0];
2714                     }else {
2715                       xtrg[maxidg]   = (Double_t) mCor->fY[3];
2716                       ytrg[maxidg]   = (Double_t) mCor->fX[3];
2717                     }
2718                     clsize1[maxidg] = mult1;
2719                     clsize2[maxidg] = mult2;
2720                     izch[maxidg]   = ich+1;
2721
2722                     Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
2723                     gAlice->ResetDigits();
2724                     if (flag==0) {
2725                       //gAlice->TreeD()->GetEvent(2); // cathode 2
2726                       gAlice->TreeD()->GetEvent(nent-1); // cathode 2
2727                     } else {
2728                       //gAlice->TreeD()->GetEvent(1);        // cathode 1
2729                       gAlice->TreeD()->GetEvent(nent-2); // cathode 1
2730                     }
2731
2732                     TClonesArray *MUONdigits  = MUON->DigitsAddress(ich);
2733                     Int_t mul=mRaw->fMultiplicity;
2734                     for (int i=0;i<mul;i++) {
2735                       Int_t idx=mRaw->fIndexMap[i];
2736                       AliMUONdigit *dig= (AliMUONdigit*)MUONdigits->UncheckedAt(idx);
2737                       ihit=dig->fHit-1;
2738                       if (ihit >= 0) break;
2739
2740                    } // loop over indexmap
2741                     //printf("fClusterType, ihit %d %d \n",mRaw->fClusterType,ihit);
2742                     if (ihit < 0) {
2743                        xgeant[maxidg] = 0; // only for resonances
2744                        ygeant[maxidg] = 0; // only for resonances
2745                        
2746                        cx[maxidg]     = 0; // Px/P of hit
2747                        cy[maxidg]     = 0; // Py/P of hit
2748                        cz[maxidg]     = 0; // Pz/P of hit
2749                        itypg[maxidg]  = 0; // particle number
2750                        ptotg[maxidg]  = 0; // P of hit
2751                        idg[maxidg]    = 0;
2752                        pvert1g[maxidg] = 0; // Px vertex
2753                        pvert2g[maxidg] = 0; // Py vertex  
2754                        pvert3g[maxidg] = 0; // Pz vertex
2755                        zvertg[maxidg]  = 0;     
2756                     } else {
2757                     gAlice->ResetHits();
2758                     TH->GetEvent(ftrack1);
2759                     TClonesArray *MUONhits  = MUON->Hits();
2760                     AliMUONhit *mHit= (AliMUONhit*)MUONhits->
2761                                                        UncheckedAt(ihit);
2762                            TClonesArray *fPartArray = gAlice->Particles();
2763                            TParticle *Part;
2764                            Int_t nch=mHit->fChamber-1;
2765                            //printf("signal ich, nch %d %d \n",ich,nch);
2766                            if (nch==ich) {
2767                              Int_t ftrack = mHit->fTrack;
2768                              Int_t id = ((TParticle*) fPartArray->
2769                                         UncheckedAt(ftrack))->GetPdgCode();
2770                              //printf("id %d \n",id);
2771                              if (id==kMuonPlus||id==kMuonMinus) {
2772                                  xgeant[maxidg] = (Double_t) mHit->fY;
2773                                  ygeant[maxidg] = (Double_t) mHit->fX;
2774                                  cx[maxidg]     = (Double_t) mHit->fCyHit; 
2775                                  cy[maxidg]     = (Double_t) mHit->fCxHit; 
2776                                  cz[maxidg]     = (Double_t) mHit->fCzHit; 
2777
2778                                  if (id==kMuonPlus) {
2779                                    itypg[maxidg]  = 5;
2780                                  } else if (id==kMuonMinus) {
2781                                    itypg[maxidg]  = 6;
2782                                  } else itypg[maxidg]  = 0;
2783                              
2784                                  ptotg[maxidg]  = (Double_t) mHit->fPTot;  
2785                                  Part = (TParticle*) fPartArray->
2786                                                      UncheckedAt(ftrack);
2787                                  Int_t iparent = Part->GetFirstMother();
2788                                  Int_t id2;
2789                                  id2 = ((TParticle*) fPartArray->
2790                                         UncheckedAt(ftrack))->GetPdgCode();
2791                              
2792                                  if (iparent >= 0) {
2793                                    Int_t ip;
2794                                    while(1) {
2795                                      ip=((TParticle*) fPartArray->
2796                                        UncheckedAt(iparent))->GetFirstMother();
2797                                      if (ip < 0) {
2798                                        id2 = ((TParticle*) fPartArray->
2799                                            UncheckedAt(iparent))->GetPdgCode();
2800                                        break;
2801                                      } else {
2802                                        iparent = ip;
2803                                        id2 = ((TParticle*) fPartArray->
2804                                            UncheckedAt(iparent))->GetPdgCode();
2805                                      } // ip<0
2806                                    } // while
2807                                  }// iparent
2808                                  Int_t id1  = ftrack; 
2809                                  Int_t idum = ftrack1+1;
2810                              
2811                                  if (id2==443 || id2==553) {
2812                                    nres++;
2813                                    if (id2==443) id2=114;
2814                                    else id2=116;
2815                                  }
2816                                  // printf("id2 %d\n",id2);
2817                                  idg[maxidg] = 30000*id1+10000*idum+id2;
2818                                  pvert1g[maxidg] = (Double_t) Part->Py(); 
2819                                  pvert2g[maxidg] = (Double_t) Part->Px();   
2820                                  pvert3g[maxidg] = (Double_t) Part->Pz(); 
2821                                  zvertg[maxidg]  = (Double_t) Part->Vz();  
2822                              } //if muon                             
2823                            } //if nch
2824                     } // ihit
2825                     maxidg++;
2826                   } // check limits
2827                 } // if cluster type
2828           } // icor loop
2829       } // ich loop
2830   }// if MUON
2831
2832
2833   ievr = nev;
2834   cout<<"evenement "<<ievr<<endl;
2835   nhittot1 = maxidg ;
2836   cout<<"nhittot1="<<nhittot1<<endl;
2837
2838   static Int_t nbres=0;
2839   static Int_t nbcor=0; 
2840   if (nres>=19) nbres++;
2841   printf("nres ,nncor - %d %d\n",nres,nncor);
2842   printf("nbres - %d\n",nbres);
2843   if (nncor>=20) nbcor++;
2844   printf("nbcor - %d\n",nbcor);
2845   printf("nuncor - %d\n",nuncor);
2846   printf("nbadcor - %d\n",nbadcor);
2847   
2848   TC->Reset();
2849
2850   hfile_global->cd();
2851   
2852 }
2853
2854 void trackf_fit(Int_t &ivertex, Double_t *pest, Double_t *pstep, Double_t &pxzinv, Double_t &tphi, Double_t &talam, Double_t &xvert, Double_t &yvert)
2855 {
2856   //
2857   //  Fit a track candidate with the following input parameters: 
2858   //  INPUT :  IVERTEX  : vertex flag, if IVERTEX=1 (XVERT,YVERT) are free paramaters
2859   //                                   if IVERTEX=1 (XVERT,YVERT)=(0.,0.) 
2860   //           PEST(5)  : starting value of parameters (minuit)
2861   //           PSTEP(5) : step size for parameters (minuit)
2862   //  OUTPUT : PXZINV,TPHI,TALAM,XVERT,YVERT : fitted value of the parameters
2863
2864   static Double_t arglist[10];
2865   static Double_t c[5] = {0.4, 0.45, 0.45, 90., 90.};
2866   static Double_t b1, b2, epxz, efi, exs, exvert, eyvert;
2867   TString chname;
2868   Int_t ierflg = 0;
2869   
2870   TMinuit *gMinuit = new TMinuit(5);
2871   gMinuit->mninit(5,10,7);
2872   gMinuit->SetFCN(fcnfwrap);  // constant m.f.
2873
2874   arglist[0] = -1;
2875   
2876   gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
2877   //      gMinuit->mnseti('track fitting');
2878   
2879   gMinuit->mnparm(0, "invmom",  pest[0], pstep[0], -c[0], c[0], ierflg);
2880   gMinuit->mnparm(1, "azimuth", pest[1], pstep[1], -c[1], c[1], ierflg);
2881   gMinuit->mnparm(2, "deep",    pest[2], pstep[2], -c[2], c[2], ierflg);
2882   if (ivertex==1) {
2883     gMinuit->mnparm(3, "x ", pest[3], pstep[3], -c[3], c[3], ierflg);
2884     gMinuit->mnparm(4, "y ", pest[4], pstep[4], -c[4], c[4], ierflg);  
2885   }   
2886   
2887   gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
2888   gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
2889   gMinuit->mnexcm("EXIT" , arglist, 0, ierflg);
2890   
2891   gMinuit->mnpout(0, chname, pxzinv, epxz, b1, b2, ierflg);
2892   gMinuit->mnpout(1, chname, tphi, efi, b1, b2, ierflg);
2893   gMinuit->mnpout(2, chname, talam, exs, b1, b2, ierflg);
2894   if (ivertex==1) {
2895     gMinuit->mnpout(3, chname, xvert, exvert, b1, b2, ierflg);
2896     gMinuit->mnpout(4, chname, yvert, eyvert, b1, b2, ierflg);
2897   }   
2898   
2899   delete gMinuit;
2900
2901 }
2902            
2903 void fcnf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *pest, Int_t iflag)
2904 {
2905   //
2906   // function called by trackf_fit
2907       Int_t futil = 0;
2908       fcn(npar,grad,fval,pest,iflag,futil);
2909 }
2910
2911 void prec_fit(Double_t &pxzinv, Double_t &fis, Double_t &alams, Double_t &xvert, Double_t &yvert, Double_t &pxzinvf, Double_t &fif, Double_t &alf, Double_t &xvertf, Double_t &yvertf, Double_t &epxzinv, Double_t &efi, Double_t &exs, Double_t &exvert, Double_t &eyvert)
2912 {
2913   // 
2914   // minuit fits for tracking finding 
2915                                                                             
2916       static Double_t arglist[10];
2917       static Double_t c1[5] = {0.001, 0.001, 0.001, 1., 1.};
2918       static Double_t c2[5] = {0.5, 0.5, 0.5, 120., 120.};
2919       static Double_t emat[9];
2920       static Double_t b1, b2;
2921       Double_t fmin, fedm, errdef; 
2922       Int_t npari, nparx, istat;
2923
2924       TString chname;
2925       Int_t ierflg = 0;
2926       
2927       TMinuit *gMinuit = new TMinuit(5);
2928       gMinuit->mninit(5,10,7);
2929       gMinuit->SetFCN(fcnfitfwrap);
2930
2931       arglist[0] = -1.;
2932       gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
2933       
2934       //      gMinuit->mnseti('track fitting');
2935
2936       gMinuit->mnparm(0,"invmom",   pxzinv, c1[0], -c2[0], c2[0], ierflg); // 0.003, 0.5
2937       gMinuit->mnparm(1,"azimuth ", fis,    c1[1], -c2[1], c2[1], ierflg);
2938       gMinuit->mnparm(2,"deep    ", alams,  c1[2], -c2[2], c2[2], ierflg);
2939       gMinuit->mnparm(3,"xvert",    xvert,  c1[3], -c2[3], c2[3], ierflg);
2940       gMinuit->mnparm(4,"yvert",    yvert,  c1[4], -c2[4], c2[4], ierflg);
2941
2942       gMinuit->mnexcm("SET NOGR", arglist, 0, ierflg);
2943       arglist[0] = 2.;
2944       gMinuit->mnexcm("MINIMIZE", arglist, 0, ierflg);
2945       gMinuit->mnexcm("EXIT", arglist, 0, ierflg);
2946  
2947       gMinuit->mnpout(0, chname, pxzinvf, epxzinv, b1, b2, ierflg);
2948       gMinuit->mnpout(1, chname, fif, efi, b1, b2, ierflg);
2949       gMinuit->mnpout(2, chname, alf, exs, b1, b2, ierflg);
2950       gMinuit->mnpout(3, chname, xvertf, exvert, b1, b2, ierflg);
2951       gMinuit->mnpout(4, chname, yvertf, eyvert, b1, b2, ierflg);
2952   
2953       gMinuit->mnemat(emat, 3);
2954       gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
2955
2956      delete gMinuit;
2957 }
2958
2959 void fcnfitf(Int_t &npar, Double_t *grad, Double_t &fval, Double_t *xval, Int_t iflag)
2960 {
2961   //
2962   // function called by prec_fit 
2963       Int_t futil = 0;
2964       fcnfit(npar,grad,fval,xval,iflag,futil);
2965 }
2966
2967 ///////////////////// fin modifs perso //////////////////////
2968
2969 ClassImp(AliMUONcluster)
2970  
2971 //___________________________________________
2972 AliMUONcluster::AliMUONcluster(Int_t *clhits)
2973 {
2974    fHitNumber=clhits[0];
2975    fCathode=clhits[1];
2976    fQ=clhits[2];
2977    fPadX=clhits[3];
2978    fPadY=clhits[4];
2979    fQpad=clhits[5];
2980    fRSec=clhits[6];
2981 }
2982 ClassImp(AliMUONdigit)
2983 //_____________________________________________________________________________
2984 AliMUONdigit::AliMUONdigit(Int_t *digits)
2985 {
2986   //
2987   // Creates a MUON digit object to be updated
2988   //
2989     fPadX        = digits[0];
2990     fPadY        = digits[1];
2991     fSignal      = digits[2];
2992     fPhysics     = digits[3];
2993     fHit       = digits[4];
2994
2995 }
2996 //_____________________________________________________________________________
2997 AliMUONdigit::AliMUONdigit(Int_t *tracks, Int_t *charges, Int_t *digits)
2998 {
2999     //
3000     // Creates a MUON digit object
3001     //
3002     fPadX        = digits[0];
3003     fPadY        = digits[1];
3004     fSignal      = digits[2];
3005     fPhysics     = digits[3];
3006     fHit       = digits[4];
3007     for(Int_t i=0; i<10; i++) {
3008         fTcharges[i]  = charges[i];
3009         fTracks[i]    = tracks[i];
3010     }
3011 }
3012
3013 AliMUONdigit::~AliMUONdigit()
3014 {
3015     
3016 }
3017
3018 ClassImp(AliMUONlist)
3019  
3020 //____________________________________________________________________________
3021     AliMUONlist::AliMUONlist(Int_t ich, Int_t *digits): 
3022         AliMUONdigit(digits)
3023 {
3024     //
3025     // Creates a MUON digit list object
3026     //
3027
3028     fChamber     = ich;
3029     fTrackList   = new TObjArray;
3030  
3031 }
3032
3033 ClassImp(AliMUONhit)
3034  
3035 //___________________________________________
3036     AliMUONhit::AliMUONhit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
3037         AliHit(shunt, track)
3038 {
3039     fChamber=vol[0];
3040     fParticle=hits[0];
3041     fX=hits[1];
3042     fY=hits[2];
3043     fZ=hits[3];
3044     fTheta=hits[4];
3045     fPhi=hits[5];
3046     fTlength=hits[6];
3047     fEloss=hits[7];
3048     fPHfirst=(Int_t) hits[8];
3049     fPHlast=(Int_t) hits[9];
3050
3051     // modifs perso
3052     fPTot=hits[10];
3053     fCxHit=hits[11];
3054     fCyHit=hits[12];
3055     fCzHit=hits[13];
3056 }
3057 ClassImp(AliMUONcorrelation)
3058 //___________________________________________
3059 //_____________________________________________________________________________
3060 AliMUONcorrelation::AliMUONcorrelation(Int_t *idx, Float_t *x, Float_t *y)
3061 {
3062     //
3063     // Creates a MUON correlation object
3064     //
3065     for(Int_t i=0; i<4; i++) {
3066         fCorrelIndex[i]  = idx[i];
3067         fX[i]    = x[i];
3068         fY[i]    = y[i];
3069     }
3070 }
3071 ClassImp(AliMUONRawCluster)
3072 Int_t AliMUONRawCluster::Compare(TObject *obj)
3073 {
3074   /*
3075          AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
3076          Float_t r=GetRadius();
3077          Float_t ro=raw->GetRadius();
3078          if (r>ro) return 1;
3079          else if (r<ro) return -1;
3080          else return 0;
3081   */
3082          AliMUONRawCluster *raw=(AliMUONRawCluster *)obj;
3083          Float_t y=fY;
3084          Float_t yo=raw->fY;
3085          if (y>yo) return 1;
3086          else if (y<yo) return -1;
3087          else return 0;
3088
3089 }
3090
3091 Int_t AliMUONRawCluster::
3092 BinarySearch(Float_t y, TArrayF coord, Int_t from, Int_t upto)
3093 {
3094    // Find object using a binary search. Array must first have been sorted.
3095    // Search can be limited by setting upto to desired index.
3096
3097    Int_t low=from, high=upto-1, half;
3098    while(high-low>1) {
3099         half=(high+low)/2;
3100         if(y>coord[half]) low=half;
3101         else high=half;
3102    }
3103    return low;
3104 }
3105
3106 void AliMUONRawCluster::SortMin(Int_t *idx,Float_t *xdarray,Float_t *xarray,Float_t *yarray,Float_t *qarray, Int_t ntr)
3107 {
3108   //
3109   // Get the 3 closest points(cog) one can find on the second cathode 
3110   // starting from a given cog on first cathode
3111   //
3112   
3113   //
3114   //  Loop over deltax, only 3 times
3115   //
3116   
3117     Float_t xmin;
3118     Int_t jmin;
3119     Int_t id[3] = {-2,-2,-2};
3120     Float_t jx[3] = {0.,0.,0.};
3121     Float_t jy[3] = {0.,0.,0.};
3122     Float_t jq[3] = {0.,0.,0.};
3123     Int_t jid[3] = {-2,-2,-2};
3124     Int_t i,j,imax;
3125   
3126     if (ntr<3) imax=ntr;
3127     else imax=3;
3128     for(i=0;i<imax;i++){
3129         xmin=1001.;
3130         jmin=0;
3131     
3132         for(j=0;j<ntr;j++){
3133             if ((i == 1 && j == id[i-1]) 
3134                   ||(i == 2 && (j == id[i-1] || j == id[i-2]))) continue;
3135            if (TMath::Abs(xdarray[j]) < xmin) {
3136               xmin = TMath::Abs(xdarray[j]);
3137               jmin=j;
3138            }       
3139         } // j
3140         if (xmin != 1001.) {    
3141            id[i]=jmin;
3142            jx[i]=xarray[jmin]; 
3143            jy[i]=yarray[jmin]; 
3144            jq[i]=qarray[jmin]; 
3145            jid[i]=idx[jmin];
3146         } 
3147     
3148     }  // i
3149   
3150     for (i=0;i<3;i++){
3151         if (jid[i] == -2) {
3152             xarray[i]=1001.;
3153             yarray[i]=1001.;
3154             qarray[i]=1001.;
3155             idx[i]=-1;
3156         } else {
3157             xarray[i]=jx[i];
3158             yarray[i]=jy[i];
3159             qarray[i]=jq[i];
3160             idx[i]=jid[i];
3161         }
3162     }
3163
3164 }
3165
3166
3167 Int_t AliMUONRawCluster::PhysicsContribution()
3168 {
3169   Int_t iPhys=0;
3170   Int_t iBg=0;
3171   Int_t iMixed=0;
3172   for (Int_t i=0; i<fMultiplicity; i++) {
3173     if (fPhysicsMap[i]==2) iPhys++;
3174     if (fPhysicsMap[i]==1) iMixed++;
3175     if (fPhysicsMap[i]==0) iBg++;
3176   }
3177   if (iMixed==0 && iBg==0) {
3178     return 2;
3179   } else if ((iPhys != 0 && iBg !=0) || iMixed != 0) {
3180     return 1;
3181   } else {
3182     return 0;
3183   }
3184 }
3185
3186    
3187 ClassImp(AliMUONreccluster) 
3188 ClassImp(AliMUONsegmentation)
3189 ClassImp(AliMUONresponse)       
3190
3191
3192
3193
3194
3195
3196
3197
3198
3199
3200
3201
3202
3203
3204
3205
3206
3207
3208