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