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