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