New version of the MUON ESD (C.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 "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       gAlice->GetMCApp()->AddHitList (fHits); // For purifyKine, only necessary when Hit list is created in AliMUONData
251   }
252   fHits = GetMUONData()->Hits(); // Added by Ivana to use the methods FisrtHit, NextHit of AliDetector
253             
254 }
255
256 //____________________________________________________________________
257 void AliMUON::SetPadSize(Int_t id, Int_t isec, Float_t p1, Float_t p2)
258 {
259 // Set the pad size for chamber id and cathode isec
260     Int_t i=2*(id-1);
261     ((AliMUONChamber*) fChambers->At(i))  ->SetPadSize(isec,p1,p2);
262     ((AliMUONChamber*) fChambers->At(i+1))->SetPadSize(isec,p1,p2);
263 }
264
265 //___________________________________________
266 void AliMUON::SetChambersZ(const Float_t *Z)
267 {
268   // Set Z values for all chambers (tracking and trigger)
269   // from the array pointed to by "Z"
270     for (Int_t ch = 0; ch < AliMUONConstants::NCh(); ch++)
271         ((AliMUONChamber*) fChambers->At(ch))->SetZ(Z[ch]);
272     return;
273 }
274 //_________________________________________________________________
275 void AliMUON::SetChambersZToDefault()
276 {
277   // Set Z values for all chambers (tracking and trigger)
278   // to default values
279   SetChambersZ(AliMUONConstants::DefaultChamberZ());
280   return;
281 }
282 //_________________________________________________________________
283 void AliMUON::SetChargeSlope(Int_t id, Float_t p1)
284 {
285 // Set the inverse charge slope for chamber id
286     Int_t i=2*(id-1);    //PH    ((AliMUONChamber*) (*fChambers)[i])->SetSigmaIntegration(p1);
287     //PH    ((AliMUONChamber*) (*fChambers)[i+1])->SetSigmaIntegration(p1);
288     ((AliMUONChamber*) fChambers->At(i))->SetChargeSlope(p1);
289     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSlope(p1);
290 }
291 //__________________________________________________________________
292 void AliMUON::SetChargeSpread(Int_t id, Float_t p1, Float_t p2)
293 {
294 // Set sigma of charge spread for chamber id
295     Int_t i=2*(id-1);
296     ((AliMUONChamber*) fChambers->At(i))->SetChargeSpread(p1,p2);
297     ((AliMUONChamber*) fChambers->At(i+1))->SetChargeSpread(p1,p2);
298 }
299 //___________________________________________________________________
300 void AliMUON::SetSigmaIntegration(Int_t id, Float_t p1)
301 {
302 // Set integration limits for charge spread
303     Int_t i=2*(id-1);
304     ((AliMUONChamber*) fChambers->At(i))->SetSigmaIntegration(p1);
305     ((AliMUONChamber*) fChambers->At(i+1))->SetSigmaIntegration(p1);
306 }
307
308 //__________________________________________________________________
309 void AliMUON::SetMaxAdc(Int_t id, Int_t p1)
310 {
311 // Set maximum number for ADCcounts (saturation)
312     Int_t i=2*(id-1);
313     ((AliMUONChamber*) fChambers->At(i))->SetMaxAdc(p1);
314     ((AliMUONChamber*) fChambers->At(i+1))->SetMaxAdc(p1);
315 }
316 //__________________________________________________________________
317 void AliMUON::SetMaxStepGas(Float_t p1)
318 {
319 // Set stepsize in gas
320   fMaxStepGas=p1;
321 }
322 //__________________________________________________________________
323 void AliMUON::SetMaxStepAlu(Float_t p1)
324 {
325 // Set step size in Alu
326     fMaxStepAlu=p1;
327 }
328 //__________________________________________________________________
329 void AliMUON::SetMaxDestepGas(Float_t p1)
330 {
331 // Set maximum step size in Gas
332     fMaxDestepGas=p1;
333 }
334 //__________________________________________________________________
335 void AliMUON::SetMaxDestepAlu(Float_t p1)
336 {
337 // Set maximum step size in Alu
338   fMaxDestepAlu=p1;
339 }
340 //___________________________________________________________________
341 void AliMUON::SetAcceptance(Bool_t acc, Float_t angmin, Float_t angmax)
342 {
343 // Set acceptance cuts 
344   fAccCut=acc;
345   fAccMin=angmin*TMath::Pi()/180;
346   fAccMax=angmax*TMath::Pi()/180;
347   Int_t ch;
348   if (acc) {
349     for (Int_t st = 0; st < AliMUONConstants::NCh() / 2; st++) {
350       // Loop over 2 chambers in the station
351       for (Int_t stCH = 0; stCH < 2; stCH++) {
352         ch = 2 * st + stCH;
353         //         Set chamber inner and outer radius according to acceptance cuts
354         Chamber(ch).SetRInner(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMin)));
355         Chamber(ch).SetROuter(TMath::Abs(AliMUONConstants::DefaultChamberZ(ch)*TMath::Tan(fAccMax)));
356       } // chamber loop
357     } // station loop
358   }
359 }
360
361 //____________________________________________________________________
362 Float_t  AliMUON::GetMaxStepGas() const
363 {
364 // Return stepsize in gas
365   
366   return fMaxStepGas;
367 }  
368
369 //____________________________________________________________________
370 Float_t  AliMUON::GetMaxStepAlu() const
371 {
372 // Return step size in Alu
373   
374   return fMaxStepAlu;
375 }
376   
377 //____________________________________________________________________
378 Float_t  AliMUON::GetMaxDestepGas() const
379 {
380 // Return maximum step size in Gas
381   
382   return fMaxDestepGas;
383 }
384   
385 //____________________________________________________________________
386 Float_t  AliMUON::GetMaxDestepAlu() const
387 {
388 // Return maximum step size in Gas
389   
390   return fMaxDestepAlu;
391 }
392
393 //____________________________________________________________________
394 void   AliMUON::SetSegmentationModel(Int_t id, Int_t isec, AliSegmentation *segmentation)
395 {
396 // Set the segmentation for chamber id cathode isec
397     ((AliMUONChamber*) fChambers->At(id))->SetSegmentationModel(isec, segmentation);
398
399 }
400 //____________________________________________________________________
401 void   AliMUON::SetResponseModel(Int_t id, AliMUONResponse *response)
402 {
403 // Set the response for chamber id
404     ((AliMUONChamber*) fChambers->At(id))->SetResponseModel(response);
405 }
406 //____________________________________________________________________
407 void   AliMUON::SetReconstructionModel(Int_t id, AliMUONClusterFinderVS *reconst)
408 {
409 // Set ClusterFinder for chamber id
410     ((AliMUONChamber*) fChambers->At(id))->SetReconstructionModel(reconst);
411 }
412 //____________________________________________________________________
413 void   AliMUON::SetNsec(Int_t id, Int_t nsec)
414 {
415 // Set number of segmented cathods for chamber id
416     ((AliMUONChamber*) fChambers->At(id))->SetNsec(nsec);
417 }
418 //____________________________________________________________________
419 AliDigitizer* AliMUON::CreateDigitizer(AliRunDigitizer* manager) const
420 {
421   return new AliMUONDigitizerv2(manager);
422 }
423 //_____________________________________________________________________
424 void AliMUON::SDigits2Digits()
425 {
426
427 // write TreeD here 
428
429     if (!fMerger) {
430       if (gAlice->GetDebug()>0) {
431         cerr<<"AliMUON::SDigits2Digits: create default AliMUONMerger "<<endl;
432         cerr<<" no merging, just digitization of 1 event will be done"<<endl;
433       }
434       fMerger = new AliMUONMerger();
435     }
436     fMerger->Init();
437     fMerger->Digitise();
438     char hname[30];
439     //    sprintf(hname,"TreeD%d",fLoader->GetHeader()->GetEvent());
440     fLoader->TreeD()->Write(hname,TObject::kOverwrite);
441     fLoader->TreeD()->Reset();
442 }
443
444 //_____________________________________________________________________
445 void AliMUON::Hits2SDigits()
446 {
447   // Adaption of AliMUONSDigitizerv1 to be excuted by the AliSimulation framework
448   AliRunLoader* runLoader = fLoader->GetRunLoader();
449   AliRunDigitizer   * manager = new AliRunDigitizer(1,1);
450   manager->SetInputStream(0,runLoader->GetFileName(),AliConfig::fgkDefaultEventFolderName);
451   AliMUONDigitizer * dMUON   = new AliMUONSDigitizerv1(manager);
452   fLoader->LoadHits("READ");
453   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
454     runLoader->GetEvent(iEvent);
455     dMUON->Exec("");
456   }
457   fLoader->UnloadHits();
458 }
459 //_______________________________________________________________________
460 AliLoader* AliMUON::MakeLoader(const char* topfoldername)
461
462 //builds standard getter (AliLoader type)
463 //if detector wants to use castomized getter, it must overload this method
464
465  if (GetDebug())
466    Info("MakeLoader",
467         "Creating standard getter for detector %s. Top folder is %s.",
468          GetName(),topfoldername);
469  fLoader   = new AliLoader(GetName(),topfoldername);
470  fMUONData = new AliMUONData(fLoader,GetName(),GetName()); 
471  fMUONData->SetSplitLevel(fSplitLevel);
472  return fLoader;
473 }
474
475 //_______________________________________________________________________
476 void AliMUON::Trigger(Int_t /*nev*/){
477 // call the Trigger Algorithm and fill TreeR
478
479   Int_t singlePlus[3]  = {0,0,0}; 
480   Int_t singleMinus[3] = {0,0,0}; 
481   Int_t singleUndef[3] = {0,0,0};
482   Int_t pairUnlike[3]  = {0,0,0}; 
483   Int_t pairLike[3]    = {0,0,0};
484   
485   ResetTrigger();
486   AliMUONTriggerDecision* decision= new AliMUONTriggerDecision(fLoader,1);
487   decision->Trigger();   
488   decision->GetGlobalTrigger(singlePlus, singleMinus, singleUndef,
489                              pairUnlike, pairLike);
490   
491   // add a local trigger in the list 
492   GetMUONData()->AddGlobalTrigger(singlePlus, singleMinus, singleUndef, pairUnlike, pairLike);
493   Int_t i;
494   
495   for (Int_t icirc=0; icirc<AliMUONConstants::NTriggerCircuit(); icirc++) { 
496     if(decision->GetITrigger(icirc)==1) {
497       Int_t localtr[7]={0,0,0,0,0,0,0};      
498       Int_t loLpt[2]={0,0}; Int_t loHpt[2]={0,0}; Int_t loApt[2]={0,0};
499       decision->GetLutOutput(icirc, loLpt, loHpt, loApt);
500       localtr[0] = icirc;
501       localtr[1] = decision->GetStripX11(icirc);
502       localtr[2] = decision->GetDev(icirc);
503       localtr[3] = decision->GetStripY11(icirc);
504       for (i=0; i<2; i++) {    // convert the Lut output in 1 digit 
505         localtr[4] = localtr[4]+Int_t(loLpt[i]*TMath::Power(2,i));
506         localtr[5] = localtr[5]+Int_t(loHpt[i]*TMath::Power(2,i));
507         localtr[6] = localtr[6]+Int_t(loApt[i]*TMath::Power(2,i));
508       }
509       GetMUONData()->AddLocalTrigger(localtr);  // add a local trigger in the list
510     }
511   }
512   
513   delete decision;
514
515   //  fLoader->TreeR()->Fill();
516   GetMUONData()->Fill("GLT"); //Filling Global and Local Trigger GLT
517   //  char hname[30];
518   //  sprintf(hname,"TreeR%d",nev);
519   //  fLoader->TreeR()->Write(hname,TObject::kOverwrite);
520     //  fLoader->TreeR()->Reset();
521   fLoader->WriteRecPoints("OVERWRITE");
522   
523   //  printf("\n End of trigger for event %d\n", nev);
524 }
525
526 //____________________________________________________________________
527 void AliMUON::Digits2Reco()
528 {
529   FindClusters();
530   Int_t nev = gAlice->GetHeader()->GetEvent();
531   GetMUONData()->Fill("RC"); //Filling Reconstructed Cluster
532   fLoader->WriteRecPoints("OVERWRITE");
533   GetMUONData()->ResetRawClusters();        
534   Info("Digits2Reco","End of cluster finding for event %d", nev);
535 }
536 //____________________________________________________________________
537 void AliMUON::FindClusters()
538 {
539 //
540 //  Perform cluster finding
541 //
542     TClonesArray *dig1, *dig2;
543     Int_t ndig, k;
544     dig1 = new TClonesArray("AliMUONDigit",1000);
545     dig2 = new TClonesArray("AliMUONDigit",1000);
546     AliMUONDigit *digit;
547 // Loop on chambers and on cathode planes
548 //
549     ResetRawClusters();        
550     TClonesArray * muonDigits;
551
552     for (Int_t ich = 0; ich < 10; ich++) {
553       //PH      AliMUONChamber* iChamber = (AliMUONChamber*) (*fChambers)[ich];
554         AliMUONChamber* iChamber = (AliMUONChamber*) fChambers->At(ich);
555         AliMUONClusterFinderVS* rec = iChamber->ReconstructionModel();
556     
557         ResetDigits();
558         GetMUONData()->GetCathode(0);
559         //TClonesArray *
560         muonDigits = GetMUONData()->Digits(ich); 
561         ndig=muonDigits->GetEntriesFast();
562         if(fDebug) 
563         printf("\n 1 Found %d digits in %p chamber %d", ndig, muonDigits,ich);
564         TClonesArray &lhits1 = *dig1;
565         Int_t n = 0;
566         for (k = 0; k < ndig; k++) {
567             digit = (AliMUONDigit*) muonDigits->UncheckedAt(k);
568             if (rec->TestTrack(digit->Track(0)))
569                 new(lhits1[n++]) AliMUONDigit(*digit);
570         }
571         GetMUONData()->ResetDigits();
572         GetMUONData()->GetCathode(1);
573         muonDigits =  GetMUONData()->Digits(ich);  
574         ndig=muonDigits->GetEntriesFast();
575         if(fDebug) 
576         printf("\n 2 Found %d digits in %p %d", ndig, muonDigits, ich);
577         TClonesArray &lhits2 = *dig2;
578         n=0;
579         
580         for (k=0; k<ndig; k++) {
581             digit= (AliMUONDigit*) muonDigits->UncheckedAt(k);
582             if (rec->TestTrack(digit->Track(0)))
583             new(lhits2[n++]) AliMUONDigit(*digit);
584         }
585
586         if (rec) {       
587             AliMUONClusterInput::Instance()->SetDigits(ich, dig1, dig2);
588             rec->FindRawClusters();
589         }
590         dig1->Delete();
591         dig2->Delete();
592     } // for ich
593     delete dig1;
594     delete dig2;
595 }
596 //______________________________________________________________________
597 #ifdef never
598 void AliMUON::Streamer(TBuffer &R__b)_
599 {
600    // Stream an object of class AliMUON.
601       AliMUONChamber        *iChamber;
602       AliMUONTriggerCircuit *iTriggerCircuit;
603       AliSegmentation       *segmentation;
604       AliMUONResponse       *response;
605       TClonesArray          *digitsaddress;
606       TClonesArray          *rawcladdress;
607       Int_t i;
608       if (R__b.IsReading()) {
609           Version_t R__v = R__b.ReadVersion(); if (R__v) { }
610           AliDetector::Streamer(R__b);
611           R__b >> fNPadHits;
612           R__b >> fPadHits; // diff
613           R__b >> fNLocalTrigger;       
614           R__b >> fLocalTrigger;       
615           R__b >> fNGlobalTrigger;       
616           R__b >> fGlobalTrigger;   
617           R__b >> fDchambers;
618           R__b >> fRawClusters;
619           R__b.ReadArray(fNdch);
620           R__b.ReadArray(fNrawch);
621           R__b >> fAccCut;
622           R__b >> fAccMin;
623           R__b >> fAccMax; 
624           R__b >> fChambers;
625           R__b >> fTriggerCircuits;
626           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
627               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
628               iTriggerCircuit->Streamer(R__b);
629           }
630 // Stream chamber related information
631           for (i =0; i<AliMUONConstants::NCh(); i++) {
632               iChamber=(AliMUONChamber*) (*fChambers)[i];
633               iChamber->Streamer(R__b);
634               if (iChamber->Nsec()==1) {
635                   segmentation=iChamber->SegmentationModel(1);
636                   if (segmentation)
637                       segmentation->Streamer(R__b);
638               } else {
639                   segmentation=iChamber->SegmentationModel(1);
640                   if (segmentation)
641                       segmentation->Streamer(R__b);
642                   if (segmentation)
643                       segmentation=iChamber->SegmentationModel(2);
644                   segmentation->Streamer(R__b);
645               }
646               response=iChamber->ResponseModel();
647               if (response)
648                   response->Streamer(R__b);       
649               digitsaddress=(TClonesArray*) (*fDchambers)[i];
650               digitsaddress->Streamer(R__b);
651               if (i < AliMUONConstants::NTrackingCh()) {
652                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
653                   rawcladdress->Streamer(R__b);
654               }
655           }
656           
657       } else {
658           R__b.WriteVersion(AliMUON::IsA());
659           AliDetector::Streamer(R__b);
660           R__b << fNPadHits;
661           R__b << fPadHits; // diff
662           R__b << fNLocalTrigger;       
663           R__b << fLocalTrigger;       
664           R__b << fNGlobalTrigger;       
665           R__b << fGlobalTrigger; 
666           R__b << fDchambers;
667           R__b << fRawClusters;
668           R__b.WriteArray(fNdch, AliMUONConstants::NCh());
669           R__b.WriteArray(fNrawch, AliMUONConstants::NTrackingCh());
670           
671           R__b << fAccCut;
672           R__b << fAccMin;
673           R__b << fAccMax; 
674           
675           R__b << fChambers;
676           R__b << fTriggerCircuits;
677           for (i =0; i<AliMUONConstants::NTriggerCircuit(); i++) {
678               iTriggerCircuit=(AliMUONTriggerCircuit*) (*fTriggerCircuits)[i];
679               iTriggerCircuit->Streamer(R__b);
680           }
681           for (i =0; i<AliMUONConstants::NCh(); i++) {
682               iChamber=(AliMUONChamber*) (*fChambers)[i];
683               iChamber->Streamer(R__b);
684               if (iChamber->Nsec()==1) {
685                   segmentation=iChamber->SegmentationModel(1);
686                   if (segmentation)
687                       segmentation->Streamer(R__b);
688               } else {
689                   segmentation=iChamber->SegmentationModel(1);
690                   if (segmentation)
691                       segmentation->Streamer(R__b);
692                   segmentation=iChamber->SegmentationModel(2);
693                   if (segmentation)
694                       segmentation->Streamer(R__b);
695               }
696               response=iChamber->ResponseModel();
697               if (response)
698                   response->Streamer(R__b);
699               digitsaddress=(TClonesArray*) (*fDchambers)[i];
700               digitsaddress->Streamer(R__b);
701               if (i < AliMUONConstants::NTrackingCh()) {
702                   rawcladdress=(TClonesArray*) (*fRawClusters)[i];
703                   rawcladdress->Streamer(R__b);
704               }
705           }
706       }
707 }
708 #endif
709 //_______________________________________________________________________
710 AliMUONPadHit* AliMUON::FirstPad(AliMUONHit*  hit, TClonesArray *clusters) 
711 {
712 // to be removed
713     // Initialise the pad iterator
714     // Return the address of the first padhit for hit
715     TClonesArray *theClusters = clusters;
716     Int_t nclust = theClusters->GetEntriesFast();
717     if (nclust && hit->PHlast() > 0) {
718         AliMUON::fMaxIterPad=hit->PHlast();
719         AliMUON::fCurIterPad=hit->PHfirst();
720         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
721     } else {
722         return 0;
723     }
724 }
725 //_______________________________________________________________________
726 AliMUONPadHit* AliMUON::NextPad(TClonesArray *clusters) 
727 {
728   // To be removed
729 // Get next pad (in iterator) 
730 //
731     AliMUON::fCurIterPad++;
732     if (AliMUON::fCurIterPad <= AliMUON::fMaxIterPad) {
733         return (AliMUONPadHit*) clusters->UncheckedAt(AliMUON::fCurIterPad-1);
734     } else {
735         return 0;
736     }
737 }
738 //_______________________________________________________________________
739
740 AliMUONRawCluster *AliMUON::RawCluster(Int_t ichamber, Int_t icathod, Int_t icluster)
741 {
742 //
743 //  Return rawcluster (icluster) for chamber ichamber and cathode icathod
744 //  Obsolete ??
745     TClonesArray *muonRawCluster  = GetMUONData()->RawClusters(ichamber);
746     ResetRawClusters();
747     TTree *treeR = fLoader->TreeR();
748     Int_t nent=(Int_t)treeR->GetEntries();
749     treeR->GetEvent(nent-2+icathod-1);
750     //treeR->GetEvent(icathod);
751     //Int_t nrawcl = (Int_t)muonRawCluster->GetEntriesFast();
752
753     AliMUONRawCluster * mRaw = (AliMUONRawCluster*)muonRawCluster->UncheckedAt(icluster);
754     //printf("RawCluster _ nent nrawcl icluster mRaw %d %d %d%p\n",nent,nrawcl,icluster,mRaw);
755     
756     return  mRaw;
757 }
758 //________________________________________________________________________
759 void   AliMUON::SetMerger(AliMUONMerger* merger)
760 {
761 // Set pointer to merger 
762     fMerger = merger;
763 }
764 //________________________________________________________________________
765 AliMUONMerger*  AliMUON::Merger()
766 {
767 // Return pointer to merger
768     return fMerger;
769 }
770 //________________________________________________________________________
771 AliMUON& AliMUON::operator = (const AliMUON& /*rhs*/)
772 {
773 // copy operator
774 // dummy version
775     return *this;
776 }
777 //________________________________________________________________________
778 void AliMUON::Reconstruct() const
779 {
780
781 //  AliLoader* loader = GetLoader();
782
783   AliRunLoader* runLoader = fLoader->GetRunLoader();
784   Int_t nEvents = runLoader->GetNumberOfEvents();
785
786 // used local container for each method
787 // passing fLoader as argument, could be avoided ???
788   AliMUONEventReconstructor* recoEvent = new AliMUONEventReconstructor(fLoader);
789   AliMUONData* dataEvent = recoEvent->GetMUONData();
790
791   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(fLoader);
792   AliMUONData* dataCluster = recoCluster->GetMUONData();
793
794   AliMUONTriggerDecision* trigDec = new AliMUONTriggerDecision(fLoader);
795   AliMUONData* dataTrig = trigDec->GetMUONData();
796
797
798   for (Int_t i = 0; i < 10; i++) {
799     AliMUONClusterFinderVS *RecModel = new AliMUONClusterFinderVS();
800     RecModel->SetGhostChi2Cut(10);
801     recoCluster->SetReconstructionModel(i,RecModel);
802   } 
803
804   fLoader->LoadDigits("READ");
805   fLoader->LoadRecPoints("RECREATE");
806   fLoader->LoadTracks("RECREATE");
807   
808   //   Loop over events              
809   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
810     printf("Event %d\n",ievent);
811     runLoader->GetEvent(ievent);
812
813     //----------------------- digit2cluster & Digits2Trigger -------------------
814     if (!fLoader->TreeR()) fLoader->MakeRecPointsContainer();
815      
816     // tracking branch
817     dataCluster->MakeBranch("RC");
818     dataCluster->SetTreeAddress("D,RC");
819     recoCluster->Digits2Clusters(); 
820     dataCluster->Fill("RC"); 
821
822     // trigger branch
823     dataTrig->MakeBranch("GLT");
824     dataTrig->SetTreeAddress("D,GLT");
825     trigDec->Digits2Trigger(); 
826     dataTrig->Fill("GLT");
827
828     fLoader->WriteRecPoints("OVERWRITE");
829
830     //---------------------------- Track & TriggerTrack ---------------------
831     if (!fLoader->TreeT()) fLoader->MakeTracksContainer();
832
833     // trigger branch
834     dataEvent->MakeBranch("RL"); //trigger track
835     dataEvent->SetTreeAddress("RL");
836     recoEvent->EventReconstructTrigger();
837     dataEvent->Fill("RL");
838
839     // tracking branch
840     dataEvent->MakeBranch("RT"); //track
841     dataEvent->SetTreeAddress("RT");
842     recoEvent->EventReconstruct();
843     dataEvent->Fill("RT");
844
845     fLoader->WriteTracks("OVERWRITE");  
846   
847     //--------------------------- Resetting branches -----------------------
848     dataCluster->ResetDigits();
849     dataCluster->ResetRawClusters();
850
851     dataTrig->ResetDigits();
852     dataTrig->ResetTrigger();
853
854     dataEvent->ResetRawClusters();
855     dataEvent->ResetTrigger();
856     dataEvent->ResetRecTracks();
857     dataEvent->ResetRecTriggerTracks();
858   
859   }
860   fLoader->UnloadDigits();
861   fLoader->UnloadRecPoints();
862   fLoader->UnloadTracks();
863
864   delete recoCluster;
865   delete recoEvent;
866   delete trigDec;
867 }
868 //________________________________________________________________________
869 void AliMUON::FillESD(AliESD* event) const
870 {
871
872   TClonesArray* recTracksArray;
873   TClonesArray* recTrigTracksArray;
874   
875   //YS AliLoader* loader = GetLoader();
876   AliRunLoader* runLoader = fLoader->GetRunLoader(); 
877   fLoader->LoadTracks("READ"); //YS
878
879
880   // declaration  
881   Int_t iEvent;
882   Int_t nTrackHits;
883   Double_t fitFmin;
884  
885
886   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
887   Double_t xRec, yRec, zRec, chi2MatchTrigger;
888   Bool_t matchTrigger;
889
890   //YS Int_t nEvents = runLoader->GetNumberOfEvents();
891
892   // setting pointer for tracks, triggertracks& trackparam at vertex
893   AliMUONTrack* recTrack;
894   AliMUONTrackParam* trackParam;
895   AliMUONTriggerTrack* recTriggerTrack;
896
897   iEvent = runLoader->GetEventNumber() ; //YS, seems not to be implemented yet (Ch. F)
898   runLoader->GetEvent(iEvent);
899
900   // setting ESD MUON class
901   AliESDMuonTrack* ESDTrack = new  AliESDMuonTrack() ;
902
903   //-------------------- trigger tracks-------------
904   Long_t trigPat = 0;
905   fMUONData->SetTreeAddress("RL");
906   fMUONData->GetRecTriggerTracks();
907   recTrigTracksArray = fMUONData->RecTriggerTracks();
908
909   // ready global trigger pattern from first track
910   recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
911   trigPat = recTriggerTrack->GetGTPattern();
912
913   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
914  
915   // -------------------- tracks-------------
916   fMUONData->SetTreeAddress("RT");
917   fMUONData->GetRecTracks();
918   recTracksArray = fMUONData->RecTracks();
919         
920   Int_t nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
921   
922   // loop over tracks
923   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
924
925     // reading info from tracks
926     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
927
928     trackParam = recTrack->GetTrackParamAtVertex();
929
930     bendingSlope            = trackParam->GetBendingSlope();
931     nonBendingSlope         = trackParam->GetNonBendingSlope();
932     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
933     xRec  = trackParam->GetNonBendingCoor();
934     yRec  = trackParam->GetBendingCoor();
935     zRec  = trackParam->GetZ();
936
937     nTrackHits       = recTrack->GetNTrackHits();
938     fitFmin          = recTrack->GetFitFMin();
939     matchTrigger     = recTrack->GetMatchTrigger();
940     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
941
942     // setting data member of ESD MUON
943     ESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
944     ESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
945     ESDTrack->SetThetaY(TMath::ATan(bendingSlope));
946     ESDTrack->SetZ(zRec);
947     ESDTrack->SetBendingCoor(yRec);
948     ESDTrack->SetNonBendingCoor(xRec);
949     ESDTrack->SetChi2(fitFmin);
950     ESDTrack->SetNHit(nTrackHits);
951     ESDTrack->SetMatchTrigger(matchTrigger);
952     ESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
953
954     // storing ESD MUON Track into ESD Event 
955     if (nRecTracks != 0)  
956       event->AddMuonTrack(ESDTrack);
957   } // end loop tracks
958
959   // add global trigger pattern
960   if (nRecTracks != 0)  
961     event->SetTrigger(trigPat);
962
963   // reset muondata
964   fMUONData->ResetRecTracks();
965   fMUONData->ResetRecTriggerTracks();
966
967   //} // end loop on event  
968   fLoader->UnloadTracks(); 
969 }
970