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