propagate cluster error parametrization
[u/mrichter/AliRoot.git] / TRD / AliTRDmodule.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 ///////////////////////////////////////////////////////////////////////////////
17 //                                                                           //
18 //                                                                           //
19 //  TRD 6-chambers stack                                                     //
20 //                                                                           //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <TObject.h>
25
26 #include "AliRun.h"
27 #include "AliLog.h"
28
29 #include "AliTRDgeometry.h"
30 #include "AliTRDmodule.h"
31 #include "AliTRDltuTracklet.h"
32 #include "AliTRDgtuTrack.h"
33 #include "AliTRDtrigParam.h"
34 #include "AliTRDzmaps.h"
35
36 ClassImp(AliTRDmodule)
37
38 //_____________________________________________________________________________
39 AliTRDmodule::AliTRDmodule()
40   :TObject()
41   ,fXprojPlane(0)
42   ,fField(0)
43   ,fTracklets(new TObjArray(400))
44   ,fTracks(new TObjArray(400))
45   ,fDeltaY(0)
46   ,fDeltaS(0)
47   ,fLTUtrk(0)
48   ,fGTUtrk(0)
49 {
50   //
51   // AliTRDmodule default constructor
52   //
53
54   fXprojPlane = AliTRDtrigParam::Instance()->GetXprojPlane();
55   fDeltaY     = AliTRDtrigParam::Instance()->GetDeltaY();
56   fDeltaS     = AliTRDtrigParam::Instance()->GetDeltaS();
57
58   // The magnetic field strength
59   Double_t x[3] = { 0.0, 0.0, 0.0 };
60   Double_t b[3];
61   gAlice->Field(x,b);  // b[] is in kilo Gauss
62   fField = b[2] * 0.1; // Tesla
63
64   for (Int_t iLayer = 0; iLayer < AliTRDgeometry::Nlayer(); iLayer++) {
65     for (Int_t i = 0; i < kNsubZchan; i++) {
66       fZnchan[iLayer][i] = 0;
67       for (Int_t j = 0; j < kNmaxZchan; j++) {
68         fZtrkid[iLayer][j][i] = -1;
69       }
70     }
71   }
72
73 }
74
75 //_____________________________________________________________________________
76 AliTRDmodule::AliTRDmodule(const AliTRDmodule &m)
77   :TObject(m)
78   ,fXprojPlane(m.fXprojPlane)
79   ,fField(m.fField)
80   ,fTracklets(NULL)
81   ,fTracks(NULL)
82   ,fDeltaY(m.fDeltaY)
83   ,fDeltaS(m.fDeltaS)
84   ,fLTUtrk(NULL)
85   ,fGTUtrk(NULL)
86 {
87   //
88   // AliTRDmodule copy constructor
89   //
90
91   for (Int_t iLayer = 0; iLayer < AliTRDgeometry::Nlayer(); iLayer++) {
92     for (Int_t i = 0; i < kNsubZchan; i++) {
93       ((AliTRDmodule &) m).fZnchan[iLayer][i] = 0;
94       for (Int_t j = 0; j < kNmaxZchan; j++) {
95         ((AliTRDmodule &) m).fZtrkid[iLayer][j][i] = -1;
96       }
97     }
98   }
99
100 }
101
102 //_____________________________________________________________________________
103 AliTRDmodule::~AliTRDmodule()
104 {
105   //
106   // Destructor
107   //
108
109 }
110
111 //_____________________________________________________________________________
112 AliTRDmodule &AliTRDmodule::operator=(const AliTRDmodule &m)
113 {
114   //
115   // Assignment operator
116   //
117
118   if (this != &m) ((AliTRDmodule &) m).Copy(*this); 
119   return *this;
120
121 }
122
123 //_____________________________________________________________________________
124 void AliTRDmodule::Copy(TObject &m) const
125 {
126   //
127   // copy function
128   //
129
130   ((AliTRDmodule &) m).fXprojPlane = fXprojPlane;
131   ((AliTRDmodule &) m).fField      = fField;
132   ((AliTRDmodule &) m).fTracklets  = NULL;
133   ((AliTRDmodule &) m).fTracks     = NULL;
134   ((AliTRDmodule &) m).fDeltaY     = fDeltaY;
135   ((AliTRDmodule &) m).fDeltaS     = fDeltaS;
136   ((AliTRDmodule &) m).fLTUtrk     = NULL;
137   ((AliTRDmodule &) m).fGTUtrk     = NULL;
138
139   for (Int_t iLayer = 0; iLayer < AliTRDgeometry::Nlayer(); iLayer++) {
140     for (Int_t i = 0; i < kNsubZchan; i++) {
141       ((AliTRDmodule &) m).fZnchan[iLayer][i] = 0;
142       for (Int_t j = 0; j < kNmaxZchan; j++) {
143         ((AliTRDmodule &) m).fZtrkid[iLayer][j][i] = -1;
144       }
145     }
146   }
147
148 }
149
150 //_____________________________________________________________________________
151 void AliTRDmodule::Reset() 
152 {
153   //
154   // Reset the tracks and tracklets in the module
155   //
156
157   ResetTracklets();
158   ResetTracks();
159
160   fLTUtrk    = 0;
161   fGTUtrk    = 0;
162   fTracklets = new TObjArray(400);
163   fTracks    = new TObjArray(400);
164
165 }
166
167 //_____________________________________________________________________________
168 void AliTRDmodule::ResetTracks() 
169 {
170   // 
171   // Reset the tracks in the module
172   //
173
174   if (fTracks) {
175
176     AliTRDgtuTrack *trk;
177     for (Int_t i = 0; i < GetNtracks(); i++) {
178       
179       trk = GetTrack(i);
180       trk->Reset();
181       
182     }
183
184     fTracks->Delete();
185
186   }
187
188 }
189
190 //_____________________________________________________________________________
191 AliTRDgtuTrack *AliTRDmodule::GetTrack(Int_t pos) const
192 {
193   //
194   // Return track at position "pos"
195   //
196
197   if (fTracks == 0) {
198     return 0;
199   }
200
201   void *trk = fTracks->UncheckedAt(pos);
202   if (trk == 0) {
203     return 0;
204   }
205
206   return (AliTRDgtuTrack *) trk;
207
208 }
209
210 //_____________________________________________________________________________
211 void AliTRDmodule::RemoveTrack(Int_t pos)
212 {
213   //
214   // Remove the track at position "pos"
215   //
216
217   if (fTracks == 0) {
218     return;
219   }
220
221   fTracks->RemoveAt(pos);
222   fTracks->Compress();
223
224 }
225
226 //_____________________________________________________________________________
227 void AliTRDmodule::AddTracklet(Int_t det, Int_t row, Float_t rowz, Float_t slope 
228                              , Float_t offset, Float_t time, Int_t ncl
229                              , Int_t label, Float_t q) 
230 {
231   // 
232   // Add a tracklet to this track
233   //
234   
235   fLTUtrk = new AliTRDltuTracklet(det,row,rowz,slope,offset,time,ncl,label,q);
236   Tracklets()->Add(fLTUtrk);
237
238 }
239
240 //_____________________________________________________________________________
241 AliTRDltuTracklet *AliTRDmodule::GetTracklet(Int_t pos) const
242 {
243   //
244   // Get the tracklet at position "pos"
245   //
246
247   if (fTracklets == 0) {
248     return 0;
249   }
250
251   void *trk = fTracklets->UncheckedAt(pos);
252   if (trk == 0) {
253     return 0;
254   }
255
256   return (AliTRDltuTracklet *) trk;
257
258 }
259
260 //_____________________________________________________________________________
261 void AliTRDmodule::RemoveTracklet(Int_t pos)
262 {
263   //
264   // Remove the tracklet at position "pos"
265   //
266
267   if (fTracklets == 0) {
268     return;
269   }
270
271   fTracklets->RemoveAt(pos);
272   fTracklets->Compress();
273
274 }
275
276 //_____________________________________________________________________________
277 void AliTRDmodule::RemoveMultipleTracklets()
278 {
279   //
280   // Remove multiple found tracklets
281   //
282
283   Float_t offDiffMin = 0.5;  // [cm]
284
285   AliTRDltuTracklet *trk;
286   Int_t   det1, det2, row1, row2, ncl1, ncl2, label1, label2;
287   Float_t off1, off2;
288   Int_t   itrk = 0;
289   while (itrk < (GetNtracklets() - 1)) {
290
291     trk    = GetTracklet(itrk);
292     det1   = trk->GetDetector();
293     row1   = trk->GetRow();
294     off1   = trk->GetOffset();
295     ncl1   = trk->GetNclusters();
296     label1 = trk->GetLabel();
297
298     trk    = GetTracklet(itrk+1);
299     det2   = trk->GetDetector();
300     row2   = trk->GetRow();
301     off2   = trk->GetOffset();
302     ncl2   = trk->GetNclusters();
303     label2 = trk->GetLabel();
304
305     if ((det1 == det2) && (row1 == row2)) {
306       if ((off2 - off1) < offDiffMin) {
307         if (ncl1 < ncl2) {
308           RemoveTracklet(itrk  );
309         }    
310         else {
311           RemoveTracklet(itrk+1);
312         }
313       }
314     }
315
316     itrk++;
317
318   }
319
320 }
321
322 //_____________________________________________________________________________
323 void AliTRDmodule::SortZ(Int_t cha)
324 {
325   //
326   // Match tracklets in the x-z plane (pad row sorting)
327   //
328
329   InitZLUT();
330
331   AliTRDltuTracklet *trk;
332   Int_t row, pla, det;
333
334   for (Int_t iTrk = 0; iTrk < GetNtracklets(); iTrk++) {
335
336     trk = GetTracklet(iTrk);
337     row = trk->GetRow();
338     det = trk->GetDetector();
339     pla = trk->GetPlane(det);
340
341     for (Int_t iZchan = 0; iZchan < kNsubZchan; iZchan++) {
342       if (fZChannelMap[cha][iZchan][pla][row] == 1) {
343         fZtrkid[pla][fZnchan[pla][iZchan]][iZchan] = iTrk;
344         fZnchan[pla][iZchan]++;
345       }
346     }
347
348   }
349
350 }
351
352 //_____________________________________________________________________________
353 void AliTRDmodule::InitZLUT()
354 {
355   //
356   // Initialize the pad row sorting look-up-table
357   //
358
359   for (Int_t iLayer = 0; iLayer < AliTRDgeometry::Nlayer(); iLayer++) {
360     for (Int_t i = 0; i < kNsubZchan; i++) {
361       fZnchan[iLayer][i] = 0;
362       for (Int_t j = 0; j < kNmaxZchan; j++) {
363         fZtrkid[iLayer][j][i] = -1;
364       }
365     }
366   }
367
368 }
369
370 //_____________________________________________________________________________
371 void AliTRDmodule::FindTracks() 
372 {
373   //
374   // Find tracks from tracklets
375   //
376
377   for (Int_t iZchan = 0; iZchan < kNsubZchan; iZchan++) {
378     FindTracksCombi(iZchan);
379   }
380
381 }
382
383 //_____________________________________________________________________________
384 void AliTRDmodule::FindTracksCombi(Int_t zchan) 
385 {
386   //
387   // Find tracks by pure combinatorics...
388   //
389   
390   static Int_t trkTrack[12];
391   
392   Int_t   nTracklets;
393   Int_t   nPlanes;
394   Int_t   ntrk1;
395   Int_t   trkId1;
396   Int_t   ntrk2;
397   Int_t   trkId2;
398
399   Float_t y1;
400   Float_t y1min;
401   Float_t y1max;
402   Float_t s1;
403   Float_t z1;
404   Float_t s1min;
405   Float_t s1max;
406   Float_t y2;
407   Float_t s2;
408   Float_t z2;
409
410   AliTRDltuTracklet *trk1;
411   AliTRDltuTracklet *trk2;
412   AliTRDltuTracklet *trk ;
413
414   Bool_t isPlane[kNplan];
415
416   for (Int_t iPlan1 = 0; iPlan1 < kNplan; iPlan1++) {
417
418     ntrk1 = fZnchan[iPlan1][zchan];
419
420     for (Int_t iTrk1 = 0; iTrk1 < ntrk1; iTrk1++) {
421
422       for (Int_t iPlan = 0; iPlan < kNplan; iPlan++) {
423         isPlane[iPlan] = kFALSE;
424       }
425
426       trkId1 = fZtrkid[iPlan1][iTrk1][zchan];
427
428       nTracklets = 0;
429       for (Int_t iList = 0; iList < kNmaxTrk; iList++) {
430         trkTrack[iList] = -1;
431       }
432       trkTrack[nTracklets++] = trkId1;
433
434       isPlane[iPlan1] = kTRUE;
435
436       trk1  = GetTracklet(trkId1);
437       y1    = trk1->GetYproj(fXprojPlane);
438       y1min = y1 - fDeltaY;
439       y1max = y1 + fDeltaY;
440       s1    = trk1->GetSlope();
441       s1min = s1 - fDeltaS;
442       s1max = s1 + fDeltaS;
443       z1    = trk1->GetZproj(fXprojPlane);      
444
445       for (Int_t iPlan2 = 0; iPlan2 < kNplan; iPlan2++) {
446
447         if (iPlan2 == iPlan1) continue;
448
449         ntrk2 = fZnchan[iPlan2][zchan];
450
451         for (Int_t iTrk2 = 0; iTrk2 < ntrk2; iTrk2++) {
452
453           trkId2 = fZtrkid[iPlan2][iTrk2][zchan];
454
455           if (trkId2 == trkId1) continue;
456
457           trk2 = GetTracklet(trkId2);
458           y2   = trk2->GetYproj(fXprojPlane);
459           s2   = trk2->GetSlope();
460           z2   = trk2->GetZproj(fXprojPlane);
461
462           if ((y1min < y2 && y2 < y1max) && 
463               (s1min < s2 && s2 < s1max)) {
464
465             if (nTracklets >= kNmaxTrk) {
466               AliWarning("Too many tracklets for this track.");
467             }    
468             else {
469               trkTrack[nTracklets++] = trkId2;
470               isPlane[iPlan2] = kTRUE;
471             }
472
473           }
474
475         }  // end trk 2
476
477       }  // end plan 2
478
479       nPlanes = 0;
480       for (Int_t iPlan = 0; iPlan < kNplan; iPlan++) {
481         nPlanes += (Int_t) isPlane[iPlan];
482       }
483       
484       if (nPlanes >= 4) {
485
486         Int_t cha1, cha2, npoints1, npoints2;
487         for (Int_t iList = 0; iList < (nTracklets - 1); iList++) {
488
489           if (trkTrack[iList] == -1 || trkTrack[iList+1] == -1) continue;
490           trk1 = GetTracklet(trkTrack[iList  ]);
491           trk2 = GetTracklet(trkTrack[iList+1]);
492
493           cha1 = trk1->GetDetector();
494           cha2 = trk2->GetDetector();
495           if (cha1 != cha2) continue;
496
497           npoints1 = trk1->GetNclusters();
498           npoints2 = trk2->GetNclusters();
499
500           if (npoints1 == npoints2) {
501             trkTrack[iList] = -1;
502           } 
503           else {
504             if (npoints1 > npoints2) trkTrack[iList+1] = -1;
505             if (npoints1 < npoints2) trkTrack[iList  ] = -1;
506           }
507
508         }
509
510         fGTUtrk = new AliTRDgtuTrack();
511         for (Int_t iList = 0; iList < nTracklets; iList++) {
512           if (trkTrack[iList] == -1) continue;
513           trk = GetTracklet(trkTrack[iList]);
514           fGTUtrk->AddTracklet(trk);
515         }
516         fGTUtrk->Track(fXprojPlane,fField);
517         AddTrack();
518
519       }
520            
521     }  // end trk 1
522
523   }  // end plan 1
524
525 }
526
527 //_____________________________________________________________________________
528 void AliTRDmodule::AddTrack() 
529 {
530   //
531   // Add a found track to the module
532   //
533   
534   Tracks()->Add(fGTUtrk);
535
536 }
537
538 //_____________________________________________________________________________
539 void AliTRDmodule::RemoveMultipleTracks()
540 {
541   //
542   // Remove multiple found tracks
543   //
544
545   AliTRDgtuTrack *trk1;
546   AliTRDgtuTrack *trk2;
547
548   Float_t yproj1;
549   Float_t yproj2;
550   Float_t alpha1;
551   Float_t alpha2;
552   Int_t   ntrk1;
553   Int_t   ntrk2;
554   Int_t   iTrack = 0;
555
556   while (iTrack < (GetNtracks()-1)) {
557
558     trk1   = GetTrack(iTrack  );
559     trk2   = GetTrack(iTrack+1);
560
561     ntrk1  = trk1->GetNtracklets();
562     yproj1 = trk1->GetYproj();
563     alpha1 = trk1->GetSlope();
564     ntrk2  = trk2->GetNtracklets();
565     yproj2 = trk2->GetYproj();
566     alpha2 = trk2->GetSlope();
567
568     if ((TMath::Abs(yproj1-yproj2) < fDeltaY) && 
569         (TMath::Abs(alpha1-alpha2) < fDeltaS)) {
570       if (ntrk1 < ntrk2) {
571         RemoveTrack(iTrack  );
572       } 
573       else {
574         RemoveTrack(iTrack+1);
575       }
576     } 
577     else {
578       iTrack++;
579     }
580     
581   }
582
583 }
584
585 //_____________________________________________________________________________
586 TObjArray *AliTRDmodule::Tracklets() 
587
588   //
589   // Returns the list of tracklets
590   //
591
592   if (!fTracklets) {
593     fTracklets = new TObjArray(400); 
594   }
595
596   return fTracklets; 
597
598 }
599
600 //_____________________________________________________________________________
601 void AliTRDmodule::ResetTracklets() 
602 {
603   //
604   // Resets the list of tracklets
605   //
606  
607   if (fTracklets) {
608     fTracklets->Delete();
609   } 
610
611 }
612
613 //_____________________________________________________________________________
614 void AliTRDmodule::SortTracklets()  
615 {
616   //
617   // Sorts the list of tracklets
618   //
619
620   if (fTracklets) {
621     fTracklets->Sort();
622   }
623  
624 }
625
626 //_____________________________________________________________________________
627 Int_t AliTRDmodule::GetNtracklets() const 
628 {
629   //
630   // Returns the number of tracklets
631   //
632
633   if (fTracklets) {
634     return fTracklets->GetEntriesFast();
635   }
636
637   return 0;
638
639 }
640
641 //_____________________________________________________________________________
642 TObjArray *AliTRDmodule::Tracks() 
643 {
644   //
645   // Returns the list of tracks
646   //
647  
648   if (!fTracks) {
649     fTracks = new TObjArray(400);
650   }
651  
652   return fTracks; 
653
654 }
655
656 //_____________________________________________________________________________
657 void AliTRDmodule::SortTracks()  
658
659   //
660   // Sort the list of tracks
661   //
662
663   if (fTracks) {
664     fTracks->Sort();
665   } 
666
667 }
668
669 //_____________________________________________________________________________
670 Int_t AliTRDmodule::GetNtracks() const 
671 {
672   //
673   // Returns the number of tracks
674   //
675
676   if (fTracks) {
677     return fTracks->GetEntriesFast();
678   }
679
680   return 0;
681
682 }