Presenting to the outside world a (x,y) reference located at the center of the slat...
[u/mrichter/AliRoot.git] / MUON / AliMUONReconstructor.cxx
1 /**************************************************************************
2
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4
5  *                                                                        *
6
7  * Author: The ALICE Off-line Project.                                    *
8
9  * Contributors are mentioned in the code where appropriate.              *
10
11  *                                                                        *
12
13  * Permission to use, copy, modify and distribute this software and its   *
14
15  * documentation strictly for non-commercial purposes is hereby granted   *
16
17  * without fee, provided that the above copyright notice appears in all   *
18
19  * copies and that both the copyright notice and this permission notice   *
20
21  * appear in the supporting documentation. The authors make no claims     *
22
23  * about the suitability of this software for any purpose. It is          *
24
25  * provided "as is" without express or implied warranty.                  *
26
27  **************************************************************************/
28
29 /* $Id$ */
30
31
32
33 ///////////////////////////////////////////////////////////////////////////////
34
35 //                                                                           //
36
37 // class for MUON reconstruction                                             //
38
39 //                                                                           //
40
41 ///////////////////////////////////////////////////////////////////////////////
42
43 #include <TParticle.h>
44
45 #include <TArrayF.h>
46
47
48
49 #include "AliRunLoader.h"
50
51 #include "AliHeader.h"
52
53 #include "AliGenEventHeader.h"
54
55 #include "AliESD.h"
56
57 #include "AliMUONReconstructor.h"
58
59  
60
61 #include "AliMUONData.h"
62
63 #include "AliMUONTrackReconstructor.h"
64
65 #include "AliMUONClusterReconstructor.h"
66
67 #include "AliMUONClusterFinderVS.h"
68
69 #include "AliMUONClusterFinderAZ.h"
70
71 #include "AliMUONTrack.h"
72
73 #include "AliMUONTrackParam.h"
74
75 #include "AliMUONTriggerTrack.h"
76
77 #include "AliESDMuonTrack.h"
78
79 #include "AliMUONRawData.h"
80
81
82
83 #include "AliRawReader.h"
84
85
86
87
88
89 ClassImp(AliMUONReconstructor)
90
91 //_____________________________________________________________________________
92
93 AliMUONReconstructor::AliMUONReconstructor()
94
95   : AliReconstructor()
96
97 {
98
99 }
100
101 //_____________________________________________________________________________
102
103 AliMUONReconstructor::~AliMUONReconstructor()
104
105 {
106
107 }
108
109 //_____________________________________________________________________________
110
111 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
112
113 {
114
115 //  AliLoader
116
117   AliLoader* loader = runLoader->GetLoader("MUONLoader");
118
119   Int_t nEvents = runLoader->GetNumberOfEvents();
120
121
122
123 // used local container for each method
124
125 // passing fLoader as argument, could be avoided ???
126
127   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
128
129   AliMUONData* dataEvent = recoEvent->GetMUONData();
130
131   if (strstr(GetOption(),"Kalman")) recoEvent->SetTrackMethod(2); // Kalman
132
133   else recoEvent->SetTrackMethod(1); // original
134
135
136
137   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
138
139   AliMUONData* dataCluster = recoCluster->GetMUONData();
140
141   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
142
143   if (strstr(GetOption(),"AZ")) {
144
145     recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ(0,1);
146
147     recoCluster->SetRecoModel(recModel);
148
149   }
150
151   recModel->SetGhostChi2Cut(10);
152
153
154
155   loader->LoadDigits("READ");
156
157   loader->LoadRecPoints("RECREATE");
158
159   loader->LoadTracks("RECREATE");
160
161   
162
163   //   Loop over events              
164
165   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
166
167     printf("Event %d\n",ievent);
168
169     runLoader->GetEvent(ievent);
170
171
172
173     //----------------------- digit2cluster & Trigger2Trigger -------------------
174
175     if (!loader->TreeR()) loader->MakeRecPointsContainer();
176
177      
178
179     // tracking branch
180
181     dataCluster->MakeBranch("RC");
182
183     dataCluster->SetTreeAddress("D,RC");
184
185     recoCluster->Digits2Clusters(); 
186
187     dataCluster->Fill("RC"); 
188
189
190
191     // trigger branch
192
193     dataCluster->MakeBranch("TC");
194
195     dataCluster->SetTreeAddress("TC");
196
197     recoCluster->Trigger2Trigger(); 
198
199     dataCluster->Fill("TC");
200
201
202
203     loader->WriteRecPoints("OVERWRITE");
204
205
206
207     //---------------------------- Track & TriggerTrack ---------------------
208
209     if (!loader->TreeT()) loader->MakeTracksContainer();
210
211
212
213     // trigger branch
214
215     dataEvent->MakeBranch("RL"); //trigger track
216
217     dataEvent->SetTreeAddress("RL");
218
219     recoEvent->EventReconstructTrigger();
220
221     dataEvent->Fill("RL");
222
223
224
225     // tracking branch
226
227     dataEvent->MakeBranch("RT"); //track
228
229     dataEvent->SetTreeAddress("RT");
230
231     recoEvent->EventReconstruct();
232
233     dataEvent->Fill("RT");
234
235
236
237     loader->WriteTracks("OVERWRITE"); 
238
239   
240
241     //--------------------------- Resetting branches -----------------------
242
243     dataCluster->ResetDigits();
244
245     dataCluster->ResetRawClusters();
246
247     dataCluster->ResetTrigger();
248
249
250
251     dataEvent->ResetRawClusters();
252
253     dataEvent->ResetTrigger();
254
255     dataEvent->ResetRecTracks();  
256
257     dataEvent->ResetRecTriggerTracks();
258
259
260
261   }
262
263   loader->UnloadDigits();
264
265   loader->UnloadRecPoints();
266
267   loader->UnloadTracks();
268
269
270
271   delete recoCluster;
272
273   delete recoEvent;
274
275 }
276
277
278
279 //_____________________________________________________________________________
280
281 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
282
283 {
284
285 //  AliLoader
286
287   AliLoader* loader = runLoader->GetLoader("MUONLoader");
288
289
290
291 // used local container for each method
292
293 // passing fLoader as argument, could be avoided ???
294
295   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
296
297   AliMUONData* dataEvent = recoEvent->GetMUONData();
298
299
300
301   AliMUONRawData* rawData = new AliMUONRawData(loader);
302
303   AliMUONData* dataCluster = rawData->GetMUONData();
304
305
306
307   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, dataCluster);
308
309   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
310
311   recModel->SetGhostChi2Cut(10);
312
313
314
315   loader->LoadRecPoints("RECREATE");
316
317   loader->LoadTracks("RECREATE");
318
319   loader->LoadDigits("RECREATE");
320
321
322
323
324
325   //   Loop over events  
326
327   Int_t iEvent = 0;
328
329             
330
331   while (rawReader->NextEvent()) {
332
333     printf("Event %d\n",iEvent);
334
335     runLoader->GetEvent(iEvent++);
336
337
338
339     //----------------------- raw2digits & raw2trigger-------------------
340
341     if (!loader->TreeD()) loader->MakeDigitsContainer();
342
343
344
345     // tracking branch
346
347     dataCluster->MakeBranch("D");
348
349     dataCluster->SetTreeAddress("D");
350
351     rawData->ReadTrackerDDL(rawReader);
352
353     dataCluster->Fill("D"); 
354
355
356
357     // trigger branch
358
359     dataCluster->MakeBranch("GLT");
360
361     dataCluster->SetTreeAddress("GLT");
362
363     rawData->ReadTriggerDDL(rawReader);
364
365     dataCluster->Fill("GLT"); 
366
367
368
369     loader->WriteDigits("OVERWRITE");
370
371
372
373     //----------------------- digit2cluster & Trigger2Trigger -------------------
374
375     if (!loader->TreeR()) loader->MakeRecPointsContainer();
376
377      
378
379     // tracking branch
380
381     dataCluster->MakeBranch("RC");
382
383     dataCluster->SetTreeAddress("RC");
384
385     recoCluster->Digits2Clusters(); 
386
387     dataCluster->Fill("RC"); 
388
389
390
391     // trigger branch
392
393     dataCluster->MakeBranch("TC");
394
395     dataCluster->SetTreeAddress("TC");
396
397     recoCluster->Trigger2Trigger(); 
398
399     dataCluster->Fill("TC");
400
401
402
403     loader->WriteRecPoints("OVERWRITE");
404
405
406
407     //---------------------------- Track & TriggerTrack ---------------------
408
409     if (!loader->TreeT()) loader->MakeTracksContainer();
410
411
412
413     // trigger branch
414
415     dataEvent->MakeBranch("RL"); //trigger track
416
417     dataEvent->SetTreeAddress("RL");
418
419     recoEvent->EventReconstructTrigger();
420
421     dataEvent->Fill("RL");
422
423
424
425     // tracking branch
426
427     dataEvent->MakeBranch("RT"); //track
428
429     dataEvent->SetTreeAddress("RT");
430
431     recoEvent->EventReconstruct();
432
433     dataEvent->Fill("RT");
434
435
436
437     loader->WriteTracks("OVERWRITE");  
438
439   
440
441     //--------------------------- Resetting branches -----------------------
442
443     dataCluster->ResetDigits();
444
445     dataCluster->ResetRawClusters();
446
447     dataCluster->ResetTrigger();
448
449
450
451     dataEvent->ResetRawClusters();
452
453     dataEvent->ResetTrigger();
454
455     dataEvent->ResetRecTracks();
456
457     dataEvent->ResetRecTriggerTracks();
458
459   
460
461   }
462
463   loader->UnloadRecPoints();
464
465   loader->UnloadTracks();
466
467   loader->UnloadDigits();
468
469
470
471   delete recoCluster;
472
473   delete recoEvent;
474
475 }
476
477
478
479 //_____________________________________________________________________________
480
481 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
482
483 {
484
485   TClonesArray* recTracksArray = 0;
486
487   TClonesArray* recTrigTracksArray = 0;
488
489   
490
491   AliLoader* loader = runLoader->GetLoader("MUONLoader");
492
493   loader->LoadTracks("READ");
494
495   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
496
497
498
499    // declaration  
500
501   Int_t iEvent;// nPart;
502
503   Int_t nTrackHits;// nPrimary;
504
505   Double_t fitFmin;
506
507   TArrayF vertex(3);
508
509
510
511   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
512
513   Double_t xRec, yRec, zRec, chi2MatchTrigger;
514
515   Bool_t matchTrigger;
516
517
518
519   // setting pointer for tracks, triggertracks & trackparam at vertex
520
521   AliMUONTrack* recTrack = 0;
522
523   AliMUONTrackParam* trackParam = 0;
524
525   AliMUONTriggerTrack* recTriggerTrack = 0;
526
527 //   TParticle* particle = new TParticle();
528
529 //   AliGenEventHeader* header = 0;
530
531   iEvent = runLoader->GetEventNumber(); 
532
533   runLoader->GetEvent(iEvent);
534
535
536
537   // vertex calculation (maybe it exists already somewhere else)
538
539   vertex[0] = vertex[1] = vertex[2] = 0.;
540
541  //  nPrimary = 0;
542
543 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
544
545 //     header->PrimaryVertex(vertex);
546
547 //   } else {
548
549 //     runLoader->LoadKinematics("READ");
550
551 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
552
553 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
554
555 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
556
557 //       runLoader->TreeK()->GetEvent(iPart);
558
559 //       if (particle->GetFirstMother() == -1) {
560
561 //      vertex[0] += particle->Vx();
562
563 //      vertex[1] += particle->Vy();
564
565 //      vertex[2] += particle->Vz();
566
567 //      nPrimary++;
568
569 //       }
570
571 //       if (nPrimary) {
572
573 //      vertex[0] /= (double)nPrimary;
574
575 //      vertex[1] /= (double)nPrimary;
576
577 //      vertex[2] /= (double)nPrimary;
578
579 //       }
580
581 //     }
582
583 //   }
584
585   // setting ESD MUON class
586
587   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
588
589
590
591   //-------------------- trigger tracks-------------
592
593   Long_t trigPat = 0;
594
595   muonData->SetTreeAddress("RL");
596
597   muonData->GetRecTriggerTracks();
598
599   recTrigTracksArray = muonData->RecTriggerTracks();
600
601
602
603   // ready global trigger pattern from first track
604
605   if (recTrigTracksArray) 
606
607     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
608
609   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
610
611
612
613   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
614
615  
616
617   // -------------------- tracks-------------
618
619   muonData->SetTreeAddress("RT");
620
621   muonData->GetRecTracks();
622
623   recTracksArray = muonData->RecTracks();
624
625         
626
627   Int_t nRecTracks = 0;
628
629   if (recTracksArray)
630
631     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
632
633   
634
635   // loop over tracks
636
637   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
638
639
640
641     // reading info from tracks
642
643     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
644
645
646
647     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
648
649     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
650
651
652
653     bendingSlope            = trackParam->GetBendingSlope();
654
655     nonBendingSlope         = trackParam->GetNonBendingSlope();
656
657     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
658
659     xRec  = trackParam->GetNonBendingCoor();
660
661     yRec  = trackParam->GetBendingCoor();
662
663     zRec  = trackParam->GetZ();
664
665
666
667     nTrackHits       = recTrack->GetNTrackHits();
668
669     fitFmin          = recTrack->GetFitFMin();
670
671     matchTrigger     = recTrack->GetMatchTrigger();
672
673     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
674
675
676
677     // setting data member of ESD MUON
678
679     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
680
681     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
682
683     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
684
685     theESDTrack->SetZ(zRec);
686
687     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
688
689     theESDTrack->SetNonBendingCoor(xRec);
690
691     theESDTrack->SetChi2(fitFmin);
692
693     theESDTrack->SetNHit(nTrackHits);
694
695     theESDTrack->SetMatchTrigger(matchTrigger);
696
697     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
698
699
700
701     // storing ESD MUON Track into ESD Event 
702
703     if (nRecTracks != 0)  
704
705       esd->AddMuonTrack(theESDTrack);
706
707   } // end loop tracks
708
709
710
711   // add global trigger pattern
712
713   if (nRecTracks != 0)  
714
715     esd->SetTrigger(trigPat);
716
717
718
719   // reset muondata
720
721   muonData->ResetRecTracks();
722
723   muonData->ResetRecTriggerTracks();
724
725
726
727   //} // end loop on event  
728
729   loader->UnloadTracks(); 
730
731  //  if (!header)
732
733 //     runLoader->UnloadKinematics();
734
735   delete theESDTrack;
736
737   delete muonData;
738
739   // delete particle;
740
741 }//_____________________________________________________________________________
742
743 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
744
745 {
746
747   // don't need rawReader ???
748
749   FillESD(runLoader, esd);
750
751 }
752