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