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