MUON/AliMUON.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18
  19 ///////////////////////////////////////////////
  20 //  Manager and hits classes for set:MUON     //
  21 ////////////////////////////////////////////////
  22
  23 #include "Riostream.h"
  24
  25 #include <AliPDG.h>
  26 #include <TBRIK.h>
  27 #include <TCanvas.h>
  28 #include <TDirectory.h>
  29 #include <TFile.h>
  30 #include <TGeometry.h>
  31 #include <TMinuit.h>
  32 #include <TNode.h>
  33 #include <TNtuple.h>
  34 #include <TObjArray.h>
  35 #include <TObject.h>
  36 #include <TObjectTable.h>
  37 #include <TPad.h>
  38 #include <TParticle.h>
  39 #include <TROOT.h>
  40 #include <TRandom.h>
  41 #include <TRotMatrix.h>
  42 #include <TTUBE.h>
  43 #include <TTUBE.h>
  44 #include <TTree.h>
  45 #include <TVector.h>
  46 #include <TVirtualMC.h>
  47
  48 #include "AliConst.h"
  49 #include "AliHeader.h"
  50 #include "AliHitMap.h"
  51 #include "AliLoader.h"
  52 #include "AliRunDigitizer.h"
  53 #include "AliESD.h"
  54 #include "AliESDMuonTrack.h"
  55 #include "AliMC.h"
  56 #include "AliMUONLoader.h"
  57 #include "AliMUON.h"
  58 #include "AliMUONTriggerTrack.h"
  59 #include "AliMUONEventReconstructor.h"
  60 #include "AliMUONClusterReconstructor.h"
  61 #include "AliMUONTrack.h"
  62 #include "AliMUONTrackParam.h"
  63 #include "AliMUONChamberTrigger.h"
  64 #include "AliMUONClusterFinderAZ.h"
  65 #include "AliMUONClusterInput.h"
  66 #include "AliMUONConstants.h"
  67 #include "AliMUONDigit.h"
  68 #include "AliMUONGlobalTrigger.h"
  69 #include "AliMUONHit.h"
  70 #include "AliMUONHitMapA1.h"
  71 #include "AliMUONLocalTrigger.h"
  72 #include "AliMUONMerger.h"
  73 #include "AliMUONPadHit.h"
  74 #include "AliMUONRawCluster.h"
  75 #include "AliMUONTransientDigit.h"
  76 #include "AliMUONTriggerCircuit.h"
  77 #include "AliMUONTriggerDecision.h"
  78 #include "AliMUONVGeometryBuilder.h"
  79 #include "AliRun.h"
  80 #include "AliMUONDigitizerv2.h"
  81 #include "AliMUONSDigitizerv1.h"
  82
  83
  84 // Defaults parameters for Z positions of chambers
  85 // taken from values for "stations" in AliMUON::AliMUON
  86 //     const Float_t zch[7]={528, 690., 975., 1249., 1449., 1610, 1710.};
  87 // and from array "dstation" in AliMUONv1::CreateGeometry
  88 //          Float_t dstation[5]={20., 20., 20, 20., 20.};
  89 //     for tracking chambers,
  90 //          according to (Z1 = zch - dstation) and  (Z2 = zch + dstation)
  91 //          for the first and second chambers in the station, respectively,
  92 // and from "DTPLANES" in AliMUONv1::CreateGeometry
  93 //           const Float_t DTPLANES = 15.;
  94 //     for trigger chambers,
  95 //          according to (Z1 = zch) and  (Z2 = zch + DTPLANES)
  96 //          for the first and second chambers in the station, respectively
  97
  98 ClassImp(AliMUON)
  99 //__________________________________________________________________
 100 AliMUON::AliMUON()
 101 {
 102 // Default Constructor
 103 //
 104     fNCh             = 0;
 105     fNTrackingCh     = 0;
 106     fIshunt          = 0;
 107     fChambers        = 0;
 108     fGeometryBuilders = 0;
 109     fTriggerCircuits = 0;
 110     fAccMin          = 0.;
 111     fAccMax          = 0.;
 112     fAccCut          = kFALSE;
 113     fMerger          = 0;
 114     fFileName        = 0;
 115     fMUONData        = 0;
 116     fSplitLevel      = 0;
 117 }
 118 //__________________________________________________________________
 119 AliMUON::AliMUON(const char *name, const char *title)
 120   : AliDetector(name,title)
 121 {
 122 //Begin_Html
 123 /*
 124 <img src="gif/alimuon.gif">
 125 */
 126 //End_Html
 127   fMUONData  = 0x0;
 128   fSplitLevel= 0;
 129   fIshunt     =  0;
 130
 131   fNCh             = AliMUONConstants::NCh();
 132   fNTrackingCh     = AliMUONConstants::NTrackingCh();
 133
 134   SetMarkerColor(kRed);//
 135 //
 136 // Creating List of Chambers
 137     Int_t ch;
 138     fChambers = new TObjArray(AliMUONConstants::NCh());
 139     fGeometryBuilders = new TObjArray(AliMUONConstants::NCh());
 140
 141     // Loop over stations
 142     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
 143       // Loop over 2 chambers in the station
 144       for (Int_t stCH = 0; stCH < 2; stCH++) {
 145         //
 146         //
 147         //    Default Parameters for Muon Tracking Stations
 148         ch = 2 * st + stCH;
 149         if (ch < AliMUONConstants::NTrackingCh()) {
 150           fChambers->AddAt(new AliMUONChamber(ch),ch);
 151         } else {
 152           fChambers->AddAt(new AliMUONChamberTrigger(ch),ch);
 153         }
 154         AliMUONChamber* chamber = (AliMUONChamber*) fChambers->At(ch);
 155         //chamber->SetGid(0);
 156         // Default values for Z of chambers
 157         chamber->SetZ(AliMUONConstants::DefaultChamberZ(ch));
 158         //
 159         chamber->InitGeo(AliMUONConstants::DefaultChamberZ(ch));
 160         //          Set chamber inner and outer radius to default
 161         chamber->SetRInner(AliMUONConstants::Dmin(st)/2);
 162         chamber->SetROuter(AliMUONConstants::Dmax(st)/2);
 163         //
 164       } // Chamber stCH (0, 1) in
 165     }     // Station st (0...)
 166
 167     // Negatives values are ignored by geant3 CONS200 in the calculation of the tracking parameters
 168     fMaxStepGas=0.1;
 169     fMaxStepAlu=0.1;
 170     fMaxDestepGas=-1;
 171     fMaxDestepAlu=-1;
 172
 173     fMaxIterPad   = 0;
 174     fCurIterPad   = 0;
 175
 176     fAccMin          = 0.;
 177     fAccMax          = 0.;
 178     fAccCut          = kFALSE;
 179
 180     // cp new design of AliMUONTriggerDecision
 181     fTriggerCircuits = new TObjArray(AliMUONConstants::NTriggerCircuit());
 182     for (Int_t circ=0; circ<AliMUONConstants::NTriggerCircuit(); circ++) {
 183       fTriggerCircuits->AddAt(new AliMUONTriggerCircuit(),circ);
 184     }
 185     fMerger = 0;
 186 }
 187 //____________________________________________________________________
 188 AliMUON::AliMUON(const AliMUON& rMUON):AliDetector(rMUON)
 189 {
 190 // Dummy copy constructor
 191     ;
 192
 193 }
 194 //____________________________________________________________________
 195 AliMUON::~AliMUON()
 196 {
 197 // Destructor
 198   if(fDebug) printf("%s: Calling AliMUON destructor !!!\n",ClassName());
 199   fIshunt  = 0;
 200   if (fMerger) delete fMerger;
 201
 202   if (fGeometryBuilders){
 203     fGeometryBuilders->Delete();
 204     delete fGeometryBuilders;
 205   }
 206   if (fChambers){
 207     fChambers->Delete();
 208     delete fChambers;
 209   }
 210   if (fTriggerCircuits){
 211     fTriggerCircuits->Delete();
 212     delete fTriggerCircuits;
 213   }
 214   delete fMUONData;
 215 }
 216 //_____________________________________________________________________________
 217 void AliMUON::AddGeometryBuilder(AliMUONVGeometryBuilder* geomBuilder)
 218 {
 219 // Adds the geometry builder to the list
 220 // ---
 221
 222   fGeometryBuilders->Add(geomBuilder);
 223 }
 224 //____________________________________________________________________
 225 void AliMUON::BuildGeometry()
 226 {
 227 // Geometry for event display
 228   for (Int_t i=0; i<7; i++) {
 229     for (Int_t j=0; j<2; j++) {
 230       Int_t id=2*i+j+1;
 231       this->Chamber(id-1).SegmentationModel(1)->Draw("eventdisplay");
 232     }
 233   }
 234 }
 235 //___________________________________________________________________
 236 Int_t AliMUON::DistancetoPrimitive(Int_t , Int_t )
 237 {
 238   return 9999;
 239 }
 240 //__________________________________________________________________
 241 void  AliMUON::SetTreeAddress()
 242 {
 243   GetMUONData()->SetLoader(fLoader);
 244   //  GetMUONData()->MakeBranch("D,S,RC");
 245   //  GetMUONData()->SetTreeAddress("H,D,S,RC");
 246   GetMUONData()->SetTreeAddress("H");
 247   if (fHits !=  GetMUONData()->Hits())  {
 248     if ( gAlice->GetMCApp() )
 249       if ( gAlice->GetMCApp()->GetHitLists() ) {
 250         fHits = GetMUONData()->Hits();
 251         gAlice->GetMCApp()->AddHitList(fHits); // For purifyKine, only necessary when Hit list is created in AliMUONData
 252       }
 253   }
 254   fHits = GetMUONData()->Hits(); // Added by Ivana to use the methods FisrtHit, NextHit of AliDetector
 255 }
 256
 257 //____________________________________________________________________
 258 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
 259 {
 260 // Set the pad size for chamber id and cathode isec
 261     Int_t i=2*(id-1);
 262     ((AliMUONChamber*) fChambers->At(i))  ->SetPadSize(isec,p1,p2);
 263     ((AliMUONChamber*) fChambers->At(i+1))->SetPadSize(isec,p1,p2);
 264 }
 265
 266 //___________________________________________
 267 void AliMUON::SetChambersZ(const Float_t *Z)
 268 {
 269   // Set Z values for all chambers (tracking and trigger)
 270   // from the array pointed to by "Z"
 271     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
 272         ((AliMUONChamber*) fChambers->At(ch))->SetZ(Z[ch]);
 273     return;
 274 }
 275 //_________________________________________________________________
 276 void AliMUON::SetChambersZToDefault()
 277 {
 278   // Set Z values for all chambers (tracking and trigger)
 279   // to default values
 280   SetChambersZ(AliMUONConstants::DefaultChamberZ());
 281   return;
 282 }
 283 //_________________________________________________________________
 284 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
 285 {
 286 // Set the inverse charge slope for chamber id
 287     Int_t i=2*(id-1);    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
 288     //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
 289     ((AliMUONChamber*) fChambers->At(i))->SetChargeSlope(p1);
 290     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSlope(p1);
 291 }
 292 //__________________________________________________________________
 293 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
 294 {
 295 // Set sigma of charge spread for chamber id
 296     Int_t i=2*(id-1);
 297     ((AliMUONChamber*) fChambers->At(i))->SetChargeSpread(p1,p2);
 298     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSpread(p1,p2);
 299 }
 300 //___________________________________________________________________
 301 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
 302 {
 303 // Set integration limits for charge spread
 304     Int_t i=2*(id-1);
 305     ((AliMUONChamber*) fChambers->At(i))->SetSigmaIntegration(p1);
 306     ((AliMUONChamber*) fChambers->At(i+1))->SetSigmaIntegration(p1);
 307 }
 308
 309 //__________________________________________________________________
 310 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
 311 {
 312 // Set maximum number for ADCcounts (saturation)
 313     Int_t i=2*(id-1);
 314     ((AliMUONChamber*) fChambers->At(i))->SetMaxAdc(p1);
 315     ((AliMUONChamber*) fChambers->At(i+1))->SetMaxAdc(p1);
 316 }
 317 //__________________________________________________________________
 318 void AliMUON::SetMaxStepGas(Float_t p1)
 319 {
 320 // Set stepsize in gas
 321   fMaxStepGas=p1;
 322 }
 323 //__________________________________________________________________
 324 void AliMUON::SetMaxStepAlu(Float_t p1)
 325 {
 326 // Set step size in Alu
 327     fMaxStepAlu=p1;
 328 }
 329 //__________________________________________________________________
 330 void AliMUON::SetMaxDestepGas(Float_t p1)
 331 {
 332 // Set maximum step size in Gas
 333     fMaxDestepGas=p1;
 334 }
 335 //__________________________________________________________________
 336 void AliMUON::SetMaxDestepAlu(Float_t p1)
 337 {
 338 // Set maximum step size in Alu
 339   fMaxDestepAlu=p1;
 340 }
 341 //___________________________________________________________________
 342 void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
 343 {
 344 // Set acceptance cuts
 345   fAccCut=acc;
 346   fAccMin=angmin*TMath::Pi()/180;
 347   fAccMax=angmax*TMath::Pi()/180;
 348   Int_t ch;
 349   if (acc) {
 350     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
 351       // Loop over 2 chambers in the station
 352       for (Int_t stCH = 0; stCH < 2; stCH++) {
 353         ch = 2 * st + stCH;
 354         //         Set chamber inner and outer radius according to acceptance cuts
 355         Chamber(ch).SetRInner(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMin)));
 356         Chamber(ch).SetROuter(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMax)));
 357       } // chamber loop
 358     } // station loop
 359   }
 360 }
 361
 362 //____________________________________________________________________
 363 Float_t  AliMUON::GetMaxStepGas() const
 364 {
 365 // Return stepsize in gas
 366
 367   return fMaxStepGas;
 368 }
 369
 370 //____________________________________________________________________
 371 Float_t  AliMUON::GetMaxStepAlu() const
 372 {
 373 // Return step size in Alu
 374
 375   return fMaxStepAlu;
 376 }
 377
 378 //____________________________________________________________________
 379 Float_t  AliMUON::GetMaxDestepGas() const
 380 {
 381 // Return maximum step size in Gas
 382
 383   return fMaxDestepGas;
 384 }
 385
 386 //____________________________________________________________________
 387 Float_t  AliMUON::GetMaxDestepAlu() const
 388 {
 389 // Return maximum step size in Gas
 390
 391   return fMaxDestepAlu;
 392 }
 393
 394 //____________________________________________________________________
 395 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
 396 {
 397 // Set the segmentation for chamber id cathode isec
 398     ((AliMUONChamber*) fChambers->At(id))->SetSegmentationModel(isec, segmentation);
 399
 400 }
 401 //____________________________________________________________________
 402 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
 403 {
 404 // Set the response for chamber id
 405     ((AliMUONChamber*) fChambers->At(id))->SetResponseModel(response);
 406 }
 407 //____________________________________________________________________
 408 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
 409 {
 410 // Set ClusterFinder for chamber id
 411     ((AliMUONChamber*) fChambers->At(id))->SetReconstructionModel(reconst);
 412 }
 413 //____________________________________________________________________
 414 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
 415 {
 416 // Set number of segmented cathods for chamber id
 417     ((AliMUONChamber*) fChambers->At(id))->SetNsec(nsec);
 418 }
 419 //____________________________________________________________________
 420 AliDigitizer* AliMUON::CreateDigitizer(AliRunDigitizer* manager) const
 421 {
 422   return new AliMUONDigitizerv2(manager);
 423 }
 424 //_____________________________________________________________________
 425 void AliMUON::SDigits2Digits()
 426 {
 427
 428 // write TreeD here
 429
 430     if (!fMerger) {
 431       if (gAlice->GetDebug()>0) {
 432         cerr<<"AliMUON::SDigits2Digits: create default AliMUONMerger "<<endl;
 433         cerr<<" no merging, just digitization of 1 event will be done"<<endl;
 434       }
 435       fMerger = new AliMUONMerger();
 436     }
 437     fMerger->Init();
 438     fMerger->Digitise();
 439     char hname[30];
 440     //    sprintf(hname,"TreeD%d",fLoader->GetHeader()->GetEvent());
 441     fLoader->TreeD()->Write(hname,TObject::kOverwrite);
 442     fLoader->TreeD()->Reset();
 443 }
 444
 445 //_____________________________________________________________________
 446 void AliMUON::Hits2SDigits()
 447 {
 448   // Adaption of AliMUONSDigitizerv1 to be excuted by the AliSimulation framework
 449   AliRunLoader* runLoader = fLoader->GetRunLoader();
 450   AliRunDigitizer   * manager = new AliRunDigitizer(1,1);
 451   manager->SetInputStream(0,runLoader->GetFileName(),AliConfig::fgkDefaultEventFolderName);
 452   AliMUONDigitizer * dMUON   = new AliMUONSDigitizerv1(manager);
 453   fLoader->LoadHits("READ");
 454   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
 455     runLoader->GetEvent(iEvent);
 456     dMUON->Exec("");
 457   }
 458   fLoader->UnloadHits();
 459 }
 460 //_______________________________________________________________________
 461 AliLoader* AliMUON::MakeLoader(const char* topfoldername)
 462 {
 463 //builds standard getter (AliLoader type)
 464 //if detector wants to use castomized getter, it must overload this method
 465
 466  if (GetDebug())
 467    Info("MakeLoader",
 468         "Creating standard getter for detector %s. Top folder is %s.",
 469          GetName(),topfoldername);
 470  fLoader   = new AliLoader(GetName(),topfoldername);
 471  fMUONData = new AliMUONData(fLoader,GetName(),GetName());
 472  fMUONData->SetSplitLevel(fSplitLevel);
 473  return fLoader;
 474 }
 475
 476 //_______________________________________________________________________
 477 void AliMUON::Trigger(Int_t /*nev*/){
 478 // call the Trigger Algorithm and fill TreeR
 479
 480   Int_t singlePlus[3]  = {0,0,0};
 481   Int_t singleMinus[3] = {0,0,0};
 482   Int_t singleUndef[3] = {0,0,0};
 483   Int_t pairUnlike[3]  = {0,0,0};
 484   Int_t pairLike[3]    = {0,0,0};
 485
 486   ResetTrigger();
 487   AliMUONTriggerDecision* decision= new AliMUONTriggerDecision(fLoader,1);
 488   decision->Trigger();
 489   decision->GetGlobalTrigger(singlePlus, singleMinus, singleUndef,
 490                              pairUnlike, pairLike);
 491
 492   // add a local trigger in the list
 493   GetMUONData()->AddGlobalTrigger(singlePlus, singleMinus, singleUndef, pairUnlike, pairLike);
 494   Int_t i;
 495
 496   for (Int_t icirc=0; icirc<AliMUONConstants::NTriggerCircuit(); icirc++) {
 497     if(decision->GetITrigger(icirc)==1) {
 498       Int_t localtr[7]={0,0,0,0,0,0,0};
 499       Int_t loLpt[2]={0,0}; Int_t loHpt[2]={0,0}; Int_t loApt[2]={0,0};
 500       decision->GetLutOutput(icirc, loLpt, loHpt, loApt);
 501       localtr[0] = icirc;
 502       localtr[1] = decision->GetStripX11(icirc);
 503       localtr[2] = decision->GetDev(icirc);
 504       localtr[3] = decision->GetStripY11(icirc);
 505       for (i=0; i<2; i++) {    // convert the Lut output in 1 digit
 506         localtr[4] = localtr[4]+Int_t(loLpt[i]*TMath::Power(2,i));
 507         localtr[5] = localtr[5]+Int_t(loHpt[i]*TMath::Power(2,i));
 508         localtr[6] = localtr[6]+Int_t(loApt[i]*TMath::Power(2,i));
 509       }
 510       GetMUONData()->AddLocalTrigger(localtr);  // add a local trigger in the list
 511     }
 512   }
 513
 514   delete decision;
 515
 516   //  fLoader->TreeR()->Fill();
 517   GetMUONData()->Fill("GLT"); //Filling Global and Local Trigger GLT
 518   //  char hname[30];
 519   //  sprintf(hname,"TreeR%d",nev);
 520   //  fLoader->TreeR()->Write(hname,TObject::kOverwrite);
 521     //  fLoader->TreeR()->Reset();
 522   fLoader->WriteRecPoints("OVERWRITE");
 523
 524   //  printf("\n End of trigger for event %d\n", nev);
 525 }
 526
 527 //____________________________________________________________________
 528 void AliMUON::Digits2Reco()
 529 {
 530   FindClusters();
 531   Int_t nev = gAlice->GetHeader()->GetEvent();
 532   GetMUONData()->Fill("RC"); //Filling Reconstructed Cluster
 533   fLoader->WriteRecPoints("OVERWRITE");
 534   GetMUONData()->ResetRawClusters();
 535   Info("Digits2Reco","End of cluster finding for event %d", nev);
 536 }
 537 //____________________________________________________________________
 538 void AliMUON::FindClusters()
 539 {
 540 //
 541 //  Perform cluster finding
 542 //
 543     TClonesArray *dig1, *dig2;
 544     Int_t ndig, k;
 545     dig1 = new TClonesArray("AliMUONDigit",1000);
 546     dig2 = new TClonesArray("AliMUONDigit",1000);
 547     AliMUONDigit *digit;
 548 // Loop on chambers and on cathode planes
 549 //
 550     ResetRawClusters();
 551     TClonesArray * muonDigits;
 552
 553     for (Int_t ich = 0; ich < 10; ich++) {
 554       //PH      AliMUONChamber* iChamber = (AliMUONChamber*) (*fChambers)[ich];
 555         AliMUONChamber* iChamber = (AliMUONChamber*) fChambers->At(ich);
 556         AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();
 557
 558         ResetDigits();
 559         GetMUONData()->GetCathode(0);
 560         //TClonesArray *
 561         muonDigits = GetMUONData()->Digits(ich);
 562         ndig=muonDigits->GetEntriesFast();
 563         if(fDebug)
 564         printf("\n 1 Found %d digits in %p chamber %d", ndig, muonDigits,ich);
 565         TClonesArray &lhits1 = *dig1;
 566         Int_t n = 0;
 567         for (k = 0; k < ndig; k++) {
 568             digit = (AliMUONDigit*) muonDigits->UncheckedAt(k);
 569             if (rec->TestTrack(digit->Track(0)))
 570                 new(lhits1[n++]) AliMUONDigit(*digit);
 571         }
 572         GetMUONData()->ResetDigits();
 573         GetMUONData()->GetCathode(1);
 574         muonDigits =  GetMUONData()->Digits(ich);
 575         ndig=muonDigits->GetEntriesFast();
 576         if(fDebug)
 577         printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
 578         TClonesArray &lhits2 = *dig2;
 579         n=0;
 580
 581         for (k=0; k<ndig; k++) {
 582             digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
 583             if (rec->TestTrack(digit->Track(0)))
 584             new(lhits2[n++]) AliMUONDigit(*digit);
 585         }
 586
 587         if (rec) {
 588             AliMUONClusterInput::Instance()->SetDigits(ich, dig1, dig2);
 589             rec->FindRawClusters();
 590         }
 591         dig1->Delete();
 592         dig2->Delete();
 593     } // for ich
 594     delete dig1;
 595     delete dig2;
 596 }
 597 //______________________________________________________________________
 598 #ifdef never
 599 void AliMUON::Streamer(TBuffer &R__b)_
 600 {
 601    // Stream an object of class AliMUON.
 602       AliMUONChamber        *iChamber;
 603       AliMUONTriggerCircuit *iTriggerCircuit;
 604       AliSegmentation       *segmentation;
 605       AliMUONResponse       *response;
 606       TClonesArray          *digitsaddress;
 607       TClonesArray          *rawcladdress;
 608       Int_t i;
 609       if (R__b.IsReading()) {
 610           Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 611           AliDetector::Streamer(R__b);
 612           R__b >> fNPadHits;
 613           R__b >> fPadHits; // diff
 614           R__b >> fNLocalTrigger;
 615           R__b >> fLocalTrigger;
 616           R__b >> fNGlobalTrigger;
 617           R__b >> fGlobalTrigger;
 618           R__b >> fDchambers;
 619           R__b >> fRawClusters;
 620           R__b.ReadArray(fNdch);
 621           R__b.ReadArray(fNrawch);
 622           R__b >> fAccCut;
 623           R__b >> fAccMin;
 624           R__b >> fAccMax;
 625           R__b >> fChambers;
 626           R__b >> fTriggerCircuits;
 627           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
 628               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
 629               iTriggerCircuit->Streamer(R__b);
 630           }
 631 // Stream chamber related information
 632           for (i =0; i<AliMUONConstants::NCh(); i++) {
 633               iChamber=(AliMUONChamber*) (*fChambers)[i];
 634               iChamber->Streamer(R__b);
 635               if (iChamber->Nsec()==1) {
 636                   segmentation=iChamber->SegmentationModel(1);
 637                   if (segmentation)
 638                       segmentation->Streamer(R__b);
 639               } else {
 640                   segmentation=iChamber->SegmentationModel(1);
 641                   if (segmentation)
 642                       segmentation->Streamer(R__b);
 643                   if (segmentation)
 644                       segmentation=iChamber->SegmentationModel(2);
 645                   segmentation->Streamer(R__b);
 646               }
 647               response=iChamber->ResponseModel();
 648               if (response)
 649                   response->Streamer(R__b);
 650               digitsaddress=(TClonesArray*) (*fDchambers)[i];
 651               digitsaddress->Streamer(R__b);
 652               if (i < AliMUONConstants::NTrackingCh()) {
 653                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
 654                   rawcladdress->Streamer(R__b);
 655               }
 656           }
 657
 658       } else {
 659           R__b.WriteVersion(AliMUON::IsA());
 660           AliDetector::Streamer(R__b);
 661           R__b << fNPadHits;
 662           R__b << fPadHits; // diff
 663           R__b << fNLocalTrigger;
 664           R__b << fLocalTrigger;
 665           R__b << fNGlobalTrigger;
 666           R__b << fGlobalTrigger;
 667           R__b << fDchambers;
 668           R__b << fRawClusters;
 669           R__b.WriteArray(fNdch, AliMUONConstants::NCh());
 670           R__b.WriteArray(fNrawch, AliMUONConstants::NTrackingCh());
 671
 672           R__b << fAccCut;
 673           R__b << fAccMin;
 674           R__b << fAccMax;
 675
 676           R__b << fChambers;
 677           R__b << fTriggerCircuits;
 678           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
 679               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
 680               iTriggerCircuit->Streamer(R__b);
 681           }
 682           for (i =0; i<AliMUONConstants::NCh(); i++) {
 683               iChamber=(AliMUONChamber*) (*fChambers)[i];
 684               iChamber->Streamer(R__b);
 685               if (iChamber->Nsec()==1) {
 686                   segmentation=iChamber->SegmentationModel(1);
 687                   if (segmentation)
 688                       segmentation->Streamer(R__b);
 689               } else {
 690                   segmentation=iChamber->SegmentationModel(1);
 691                   if (segmentation)
 692                       segmentation->Streamer(R__b);
 693                   segmentation=iChamber->SegmentationModel(2);
 694                   if (segmentation)
 695                       segmentation->Streamer(R__b);
 696               }
 697               response=iChamber->ResponseModel();
 698               if (response)
 699                   response->Streamer(R__b);
 700               digitsaddress=(TClonesArray*) (*fDchambers)[i];
 701               digitsaddress->Streamer(R__b);
 702               if (i < AliMUONConstants::NTrackingCh()) {
 703                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
 704                   rawcladdress->Streamer(R__b);
 705               }
 706           }
 707       }
 708 }
 709 #endif
 710 //_______________________________________________________________________
 711 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters)
 712 {
 713 // to be removed
 714     // Initialise the pad iterator
 715     // Return the address of the first padhit for hit
 716     TClonesArray *theClusters = clusters;
 717     Int_t nclust = theClusters->GetEntriesFast();
 718     if (nclust && hit->PHlast() > 0) {
 719         AliMUON::fMaxIterPad=hit->PHlast();
 720         AliMUON::fCurIterPad=hit->PHfirst();
 721         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
 722     } else {
 723         return 0;
 724     }
 725 }
 726 //_______________________________________________________________________
 727 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters)
 728 {
 729   // To be removed
 730 // Get next pad (in iterator)
 731 //
 732     AliMUON::fCurIterPad++;
 733     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
 734         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
 735     } else {
 736         return 0;
 737     }
 738 }
 739 //_______________________________________________________________________
 740
 741 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
 742 {
 743 //
 744 //  Return rawcluster (icluster) for chamber ichamber and cathode icathod
 745 //  Obsolete ??
 746     TClonesArray *muonRawCluster  = GetMUONData()->RawClusters(ichamber);
 747     ResetRawClusters();
 748     TTree *treeR = fLoader->TreeR();
 749     Int_t nent=(Int_t)treeR->GetEntries();
 750     treeR->GetEvent(nent-2+icathod-1);
 751     //treeR->GetEvent(icathod);
 752     //Int_t nrawcl = (Int_t)muonRawCluster->GetEntriesFast();
 753
 754     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)muonRawCluster->UncheckedAt(icluster);
 755     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
 756
 757     return  mRaw;
 758 }
 759 //________________________________________________________________________
 760 void   AliMUON::SetMerger(AliMUONMerger* merger)
 761 {
 762 // Set pointer to merger
 763     fMerger = merger;
 764 }
 765 //________________________________________________________________________
 766 AliMUONMerger*  AliMUON::Merger()
 767 {
 768 // Return pointer to merger
 769     return fMerger;
 770 }
 771 //________________________________________________________________________
 772 AliMUON& AliMUON::operator = (const AliMUON& /*rhs*/)
 773 {
 774 // copy operator
 775 // dummy version
 776     return *this;
 777 }
 778 //________________________________________________________________________
 779 void AliMUON::Reconstruct() const
 780 {
 781
 782 //  AliLoader* loader = GetLoader();
 783
 784   AliRunLoader* runLoader = fLoader->GetRunLoader();
 785   Int_t nEvents = runLoader->GetNumberOfEvents();
 786
 787 // used local container for each method
 788 // passing fLoader as argument, could be avoided ???
 789   AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(fLoader);
 790   AliMUONData* dataEvent = recoEvent->GetMUONData();
 791
 792   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(fLoader);
 793   AliMUONData* dataCluster = recoCluster->GetMUONData();
 794
 795   AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(fLoader);
 796   AliMUONData* dataTrig = trigDec->GetMUONData();
 797
 798
 799   for (Int_t i = 0; i < 10; i++) {
 800     AliMUONClusterFinderVS *RecModel = new AliMUONClusterFinderVS();
 801     RecModel->SetGhostChi2Cut(10);
 802     recoCluster->SetReconstructionModel(i,RecModel);
 803   }
 804
 805   fLoader->LoadDigits("READ");
 806   fLoader->LoadRecPoints("RECREATE");
 807   fLoader->LoadTracks("RECREATE");
 808
 809   //   Loop over events
 810   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
 811     printf("Event %d\n",ievent);
 812     runLoader->GetEvent(ievent);
 813
 814     //----------------------- digit2cluster & Digits2Trigger -------------------
 815     if (!fLoader->TreeR()) fLoader->MakeRecPointsContainer();
 816
 817     // tracking branch
 818     dataCluster->MakeBranch("RC");
 819     dataCluster->SetTreeAddress("D,RC");
 820     recoCluster->Digits2Clusters();
 821     dataCluster->Fill("RC");
 822
 823     // trigger branch
 824     dataTrig->MakeBranch("GLT");
 825     dataTrig->SetTreeAddress("D,GLT");
 826     trigDec->Digits2Trigger();
 827     dataTrig->Fill("GLT");
 828
 829     fLoader->WriteRecPoints("OVERWRITE");
 830
 831     //---------------------------- Track & TriggerTrack ---------------------
 832     if (!fLoader->TreeT()) fLoader->MakeTracksContainer();
 833
 834     // trigger branch
 835     dataEvent->MakeBranch("RL"); //trigger track
 836     dataEvent->SetTreeAddress("RL");
 837     recoEvent->EventReconstructTrigger();
 838     dataEvent->Fill("RL");
 839
 840     // tracking branch
 841     dataEvent->MakeBranch("RT"); //track
 842     dataEvent->SetTreeAddress("RT");
 843     recoEvent->EventReconstruct();
 844     dataEvent->Fill("RT");
 845
 846     fLoader->WriteTracks("OVERWRITE");
 847
 848     //--------------------------- Resetting branches -----------------------
 849     dataCluster->ResetDigits();
 850     dataCluster->ResetRawClusters();
 851
 852     dataTrig->ResetDigits();
 853     dataTrig->ResetTrigger();
 854
 855     dataEvent->ResetRawClusters();
 856     dataEvent->ResetTrigger();
 857     dataEvent->ResetRecTracks();
 858     dataEvent->ResetRecTriggerTracks();
 859
 860   }
 861   fLoader->UnloadDigits();
 862   fLoader->UnloadRecPoints();
 863   fLoader->UnloadTracks();
 864
 865   delete recoCluster;
 866   delete recoEvent;
 867   delete trigDec;
 868 }
 869 //________________________________________________________________________
 870 void AliMUON::FillESD(AliESD* event) const
 871 {
 872
 873   TClonesArray* recTracksArray;
 874   TClonesArray* recTrigTracksArray;
 875
 876   //YS AliLoader* loader = GetLoader();
 877   AliRunLoader* runLoader = fLoader->GetRunLoader();
 878   fLoader->LoadTracks("READ"); //YS
 879
 880
 881   // declaration
 882   Int_t iEvent;
 883   Int_t nTrackHits;
 884   Double_t fitFmin;
 885
 886
 887   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
 888   Double_t xRec, yRec, zRec, chi2MatchTrigger;
 889   Bool_t matchTrigger;
 890
 891   //YS Int_t nEvents = runLoader->GetNumberOfEvents();
 892
 893   // setting pointer for tracks, triggertracks& trackparam at vertex
 894   AliMUONTrack* recTrack;
 895   AliMUONTrackParam* trackParam;
 896   AliMUONTriggerTrack* recTriggerTrack;
 897
 898   iEvent = runLoader->GetEventNumber() ; //YS, seems not to be implemented yet (Ch. F)
 899   runLoader->GetEvent(iEvent);
 900
 901   // setting ESD MUON class
 902   AliESDMuonTrack* ESDTrack = new  AliESDMuonTrack() ;
 903
 904   //-------------------- trigger tracks-------------
 905   Long_t trigPat = 0;
 906   fMUONData->SetTreeAddress("RL");
 907   fMUONData->GetRecTriggerTracks();
 908   recTrigTracksArray = fMUONData->RecTriggerTracks();
 909
 910   // ready global trigger pattern from first track
 911   recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
 912   trigPat = recTriggerTrack->GetGTPattern();
 913
 914   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
 915
 916   // -------------------- tracks-------------
 917   fMUONData->SetTreeAddress("RT");
 918   fMUONData->GetRecTracks();
 919   recTracksArray = fMUONData->RecTracks();
 920
 921   Int_t nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
 922
 923   // loop over tracks
 924   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 925
 926     // reading info from tracks
 927     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 928
 929     trackParam = recTrack->GetTrackParamAtVertex();
 930
 931     bendingSlope            = trackParam->GetBendingSlope();
 932     nonBendingSlope         = trackParam->GetNonBendingSlope();
 933     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
 934     xRec  = trackParam->GetNonBendingCoor();
 935     yRec  = trackParam->GetBendingCoor();
 936     zRec  = trackParam->GetZ();
 937
 938     nTrackHits       = recTrack->GetNTrackHits();
 939     fitFmin          = recTrack->GetFitFMin();
 940     matchTrigger     = recTrack->GetMatchTrigger();
 941     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 942
 943     // setting data member of ESD MUON
 944     ESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
 945     ESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
 946     ESDTrack->SetThetaY(TMath::ATan(bendingSlope));
 947     ESDTrack->SetZ(zRec);
 948     ESDTrack->SetBendingCoor(yRec);
 949     ESDTrack->SetNonBendingCoor(xRec);
 950     ESDTrack->SetChi2(fitFmin);
 951     ESDTrack->SetNHit(nTrackHits);
 952     ESDTrack->SetMatchTrigger(matchTrigger);
 953     ESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 954
 955     // storing ESD MUON Track into ESD Event
 956     if (nRecTracks != 0)
 957       event->AddMuonTrack(ESDTrack);
 958   } // end loop tracks
 959
 960   // add global trigger pattern
 961   if (nRecTracks != 0)
 962     event->SetTrigger(trigPat);
 963
 964   // reset muondata
 965   fMUONData->ResetRecTracks();
 966   fMUONData->ResetRecTriggerTracks();
 967
 968   //} // end loop on event
 969   fLoader->UnloadTracks();
 970 }
 971