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