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