MUON/AliMUON.cxx

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