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