MUON/AliMUON.cxx

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