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