1fb4b547277b055facf882e54ff86a853615350f
[u/mrichter/AliRoot.git] / TRD / AliTRDtrigger.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 trigger class                                                        //
20 //                                                                           //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <TTree.h>
25 #include <TBranch.h>
26 #include <TMatrixD.h>
27
28 #include "AliLog.h"
29 #include "AliRun.h"
30 #include "AliLoader.h"
31
32 #include "AliTRDdigitsManager.h"
33 #include "AliTRDgeometry.h"
34 #include "AliTRDdataArrayI.h"
35 #include "AliTRDcalibDB.h"
36 #include "AliTRDCommonParam.h"
37 #include "Cal/AliTRDCalPIDLQ.h"
38 #include "AliTRDrawData.h"
39
40 #include "AliTRDtrigger.h"
41 #include "AliTRDmodule.h"
42 #include "AliTRDmcmTracklet.h"
43 #include "AliTRDtrigParam.h"
44 #include "AliTRDmcm.h"
45 #include "AliTRDzmaps.h"
46
47 //#include "AliHeader.h"
48
49 ClassImp(AliTRDtrigger)
50
51 //_____________________________________________________________________________
52 AliTRDtrigger::AliTRDtrigger():
53   TNamed(),
54   fTracks("AliTRDgtuTrack",0)
55 {
56   //
57   // AliTRDtrigger default constructor
58   //
59
60   fDigitsManager = NULL;
61   fTrackletTree = NULL;
62   fTracklets     = NULL;
63
64   fNROB = 0;
65   fTrigParam = NULL;
66   fMCM = NULL;
67   fTrk = NULL;
68   fGTUtrk = NULL;
69
70   fNtracklets = 0;
71
72   fDigits = NULL;
73   fTrack0 = NULL;
74   fTrack1 = NULL;
75   fTrack2 = NULL;
76
77   fModule = NULL;
78
79   fNPrimary = 0;
80
81 }
82
83 //_____________________________________________________________________________
84 AliTRDtrigger::AliTRDtrigger(const Text_t *name, const Text_t *title):
85   TNamed(name,title),
86   fTracks("AliTRDgtuTrack",1000)
87 {
88   //
89   // AliTRDtrigger constructor
90   //
91
92   fDigitsManager = new AliTRDdigitsManager();
93   fTrackletTree = NULL;
94   fTracklets = new TObjArray(400);
95
96   fNROB = 0;
97   fTrigParam = NULL;
98   fMCM = NULL;
99   fTrk = NULL;
100   fGTUtrk = NULL;
101
102   fNtracklets = 0;
103
104   fDigits = NULL;
105   fTrack0 = NULL;
106   fTrack1 = NULL;
107   fTrack2 = NULL;
108
109   fModule = NULL;
110
111   fNPrimary = 0;
112
113 }
114
115 //_____________________________________________________________________________
116 AliTRDtrigger::AliTRDtrigger(const AliTRDtrigger &p):TNamed(p)
117 {
118   //
119   // AliTRDtrigger copy constructor
120   //
121
122   ((AliTRDtrigger &) p).Copy(*this);
123
124 }
125
126 ///_____________________________________________________________________________
127 AliTRDtrigger::~AliTRDtrigger()
128 {
129   //
130   // AliTRDtrigger destructor
131   //
132
133   if (fTracklets) {
134     fTracklets->Delete();
135     delete fTracklets;
136   }
137
138   fTracks.Delete();
139
140 }
141
142 //_____________________________________________________________________________
143 AliTRDtrigger &AliTRDtrigger::operator=(const AliTRDtrigger &p)
144 {
145   //
146   // Assignment operator
147   //
148
149   if (this != &p) ((AliTRDtrigger &) p).Copy(*this);
150   return *this;
151
152 }
153
154 //_____________________________________________________________________________
155 void AliTRDtrigger::Copy(TObject &) const
156 {
157   //
158   // Copy function
159   //
160
161   AliFatal("Not implemented");
162
163 }
164
165 //_____________________________________________________________________________
166 void AliTRDtrigger::Init()
167 {
168
169   fModule = new AliTRDmodule(fTrigParam); 
170   /*
171   AliHeader *header = fRunLoader->GetHeader();
172   fNPrimary = header->GetNprimary();
173   */
174   fTracks.Clear();
175
176 }
177
178 //_____________________________________________________________________________
179 Bool_t AliTRDtrigger::Open(const Char_t *name, Int_t nEvent)
180 {
181   //
182   // Opens the AliROOT file.
183   //
184
185   TString evfoldname = AliConfig::GetDefaultEventFolderName();
186   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
187
188   if (!fRunLoader)
189     fRunLoader = AliRunLoader::Open(name);
190
191   if (!fRunLoader) {
192     Error("Open","Can not open session for file %s.",name);
193     return kFALSE;
194   }
195
196   // Open input
197
198   if (fRunLoader->GetAliRun() == 0x0) fRunLoader->LoadgAlice();
199   gAlice = fRunLoader->GetAliRun();
200
201   if (!(gAlice)) {
202     fRunLoader->LoadgAlice();
203     gAlice = fRunLoader->GetAliRun();
204     if (!(gAlice)) {
205       Error("Open","Could not find AliRun object.");
206       return kFALSE;
207     }
208   }
209
210   // Import the Trees for the event nEvent in the file
211   fRunLoader->GetEvent(nEvent);
212
213   // Open output
214
215   TObjArray *ioArray = 0;
216
217   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
218   loader->MakeTree("T");
219   fTrackletTree = loader->TreeT();
220   fTrackletTree->Branch("TRDmcmTracklet","TObjArray",&ioArray,32000,0);
221   /*
222   fRunLoader->LoadHeader();
223   */
224   Init();
225
226   return kTRUE;
227
228 }
229
230
231 //_____________________________________________________________________________
232 Bool_t AliTRDtrigger::ReadDigits() 
233 {
234   //
235   // Reads the digits arrays from the input aliroot file
236   //
237
238   if (!fRunLoader) {
239     Error("ReadDigits","Can not find the Run Loader");
240     return kFALSE;
241   }
242
243   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
244   if (!loader->TreeD()) loader->LoadDigits();
245
246   return (fDigitsManager->ReadDigits(loader->TreeD()));
247
248 }
249
250 //_____________________________________________________________________________
251 Bool_t AliTRDtrigger::ReadDigits(AliRawReader* rawReader)
252 {
253   //
254   // Reads the digits arrays from the ddl file
255   //
256
257   AliTRDrawData *raw = new AliTRDrawData();
258   raw->SetDebug(1);
259
260   fDigitsManager = raw->Raw2Digits(rawReader);
261
262   return kTRUE;
263
264 }
265
266 //_____________________________________________________________________________
267 Bool_t AliTRDtrigger::ReadTracklets(AliRunLoader *rl) 
268 {
269   //
270   // Reads the tracklets find the tracks
271   //
272
273   Int_t idet;
274
275   AliLoader *loader = rl->GetLoader("TRDLoader");
276   loader->LoadTracks();
277   fTrackletTree = loader->TreeT();
278
279   TBranch *branch = fTrackletTree->GetBranch("TRDmcmTracklet");
280   if (!branch) {
281     Error("ReadTracklets","Can't get the branch !");
282     return kFALSE;
283   }
284   TObjArray *tracklets = new TObjArray(400);
285   branch->SetAddress(&tracklets);
286
287   Int_t nEntries = (Int_t) fTrackletTree->GetEntries();
288   Int_t iEntry, itrk;
289   Int_t iStack, iStackPrev = -1;
290   
291   for (iEntry = 0; iEntry < nEntries; iEntry++) {    
292     fTrackletTree->GetEvent(iEntry);
293     
294     for (itrk = 0; itrk < tracklets->GetEntriesFast(); itrk++){
295
296       fTrk = (AliTRDmcmTracklet*)tracklets->UncheckedAt(itrk);
297       
298       idet = fTrk->GetDetector();
299
300       iStack = idet / (AliTRDgeometry::Nplan());
301       if (iStackPrev != iStack) {
302         if (iStackPrev == -1) {
303           iStackPrev = iStack;
304         } else {
305           MakeTracks(idet-AliTRDgeometry::Nplan());
306           ResetTracklets();
307           iStackPrev = iStack;
308         }
309       }
310       
311       Tracklets()->Add(fTrk);
312
313       if (iEntry == (nEntries-1) && itrk == (tracklets->GetEntriesFast()-1)) {
314         idet++;
315         MakeTracks(idet-AliTRDgeometry::Nplan());
316         ResetTracklets();
317       }
318
319     }
320
321   }
322
323   loader->UnloadTracks();
324
325   return kTRUE;
326
327 }
328
329 //_____________________________________________________________________________
330 Bool_t AliTRDtrigger::MakeTracklets(Bool_t makeTracks)
331 {
332
333   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
334   if (!calibration)
335   {
336     Error("MakeTracklets","No instance of AliTRDcalibDB.");
337     return kFALSE;  
338   }
339   
340   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
341   if (!commonParam)
342   {
343     Error("MakeTracklets","No common params.");
344     return kFALSE;
345   }
346     
347   AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
348
349   Int_t    chamBeg = 0;
350   Int_t    chamEnd = AliTRDgeometry::Ncham();
351   Int_t    planBeg = 0;
352   Int_t    planEnd = AliTRDgeometry::Nplan();
353   Int_t    sectBeg = 0;
354   Int_t    sectEnd = AliTRDgeometry::Nsect();
355
356   fMCM = new AliTRDmcm(fTrigParam,0);
357
358   Int_t time, col, row, col1, col2;
359   Float_t amp;
360   Int_t idet, iStack, iStackPrev;
361   idet = -1;
362
363   iStack     = -1;
364   iStackPrev = -1;
365   for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
366
367     for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
368
369       // number of ROBs in the chamber
370       if( icham == 2 ) {
371         fNROB = 6;
372       } else {
373         fNROB = 8;
374       }
375
376       for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
377
378         idet = geo->GetDetector(iplan,icham,isect);
379         ResetTracklets();
380         
381         if (makeTracks) {
382           iStack = idet / (AliTRDgeometry::Nplan());
383           if (iStackPrev != iStack) {
384             if (iStackPrev == -1) {
385               iStackPrev = iStack;
386             } else {
387               MakeTracks(idet-AliTRDgeometry::Nplan());
388               ResetTracklets();
389               iStackPrev = iStack;
390             }
391           }
392         }
393
394         Int_t    nRowMax     = commonParam->GetRowMax(iplan,icham,isect);
395         Int_t    nColMax     = commonParam->GetColMax(iplan);
396         Int_t    nTimeTotal  = calibration->GetNumberOfTimeBins();
397
398         // Get the digits
399         fDigits = fDigitsManager->GetDigits(idet);
400         fDigits->Expand();
401         fTrack0 = fDigitsManager->GetDictionary(idet,0);
402         fTrack0->Expand();
403         fTrack1 = fDigitsManager->GetDictionary(idet,1);
404         fTrack1->Expand();
405         fTrack2 = fDigitsManager->GetDictionary(idet,2); 
406         fTrack2->Expand();
407
408         for (Int_t iRob = 0; iRob < fNROB; iRob++) {
409
410           for (Int_t iMcm = 0; iMcm < kNMCM; iMcm++) {
411
412             fMCM->Reset();
413
414             fMCM->SetRobId(iRob);
415             fMCM->SetChaId(idet);
416
417             SetMCMcoordinates(iMcm);
418
419             row = fMCM->GetRow();
420
421             if (row < 0 || row >= nRowMax) {
422               Error("MakeTracklets","MCM row number out of range.");
423               continue;
424             }
425
426             fMCM->GetColRange(col1,col2);
427             
428             for (time = 0; time < nTimeTotal; time++) {
429               for (col = col1; col < col2; col++) {
430                 if (col >= 0 && col < nColMax) {
431                   amp = TMath::Abs(fDigits->GetDataUnchecked(row,col,time));
432                 } else {
433                   amp = 0.0;
434                 }
435                 fMCM->SetADC(col-col1,time,amp);
436
437               }
438             }
439
440             if (fTrigParam->GetTailCancelation()) {
441               fMCM->Filter(fTrigParam->GetNexponential(),fTrigParam->GetFilterType());
442             }
443             
444             if (fMCM->Run()) {
445
446               for (Int_t iSeed = 0; iSeed < kMaxTrackletsPerMCM; iSeed++) {
447                 
448                 if (fMCM->GetSeedCol()[iSeed] < 0) continue;
449
450                 if ( fTrigParam->GetDebugLevel() > 1 ) 
451                   printf("Add tracklet %d in col %02d \n",fNtracklets,fMCM->GetSeedCol()[iSeed]);
452
453                 if ( fTrigParam->GetDebugLevel() == -1 ) {
454                   printf("Add tracklet %d in col %02d \n",fNtracklets,fMCM->GetSeedCol()[iSeed]);
455                   for (time = 0; time < nTimeTotal; time++) {
456                     for (col = 0; col < kMcmCol; col++) {                   
457                       printf("%03.0f  ",fMCM->GetADC(col,time));
458                     }
459                     printf("\n");
460                   }
461                 }
462
463                 AddTracklet(idet,row,iSeed,fNtracklets++);
464
465               }
466
467             }
468
469           }
470
471       
472         }
473
474         // Compress the arrays
475         fDigits->Compress(1,0);
476         fTrack0->Compress(1,0);
477         fTrack1->Compress(1,0);
478         fTrack2->Compress(1,0);
479
480         WriteTracklets(idet);
481
482      }
483     }
484   }
485
486   if (makeTracks) {
487     idet++;
488     MakeTracks(idet-AliTRDgeometry::Nplan());
489     ResetTracklets();
490   }
491
492   return kTRUE;
493
494 }
495
496 //_____________________________________________________________________________
497 void AliTRDtrigger::SetMCMcoordinates(Int_t imcm)
498 {
499
500   Int_t robid = fMCM->GetRobId();
501
502   // setting the Row and Col range
503
504   const Int_t kNcolRob = 2;  // number of ROBs per chamber in column direction
505   const Int_t kNmcmRob = 4;  // number of MCMs per ROB in column/row direction
506
507   Int_t mcmid = imcm%(kNmcmRob*kNmcmRob);
508
509   if (robid%kNcolRob == 0) {
510
511     if ( mcmid%kNmcmRob == 0 ) {
512       fMCM->SetColRange(18*0-1,18*1-1+2+1);
513     }
514     if ( mcmid%kNmcmRob == 1 ) {
515       fMCM->SetColRange(18*1-1,18*2-1+2+1);
516     }
517     if ( mcmid%kNmcmRob == 2 ) {
518       fMCM->SetColRange(18*2-1,18*3-1+2+1);
519     }
520     if ( mcmid%kNmcmRob == 3 ) {
521       fMCM->SetColRange(18*3-1,18*4-1+2+1);
522     }
523
524   } else {
525
526     if ( mcmid%kNmcmRob == 0 ) {
527       fMCM->SetColRange(18*4-1,18*5-1+2+1);
528     }
529     if ( mcmid%kNmcmRob == 1 ) {
530       fMCM->SetColRange(18*5-1,18*6-1+2+1);
531     }
532     if ( mcmid%kNmcmRob == 2 ) {
533       fMCM->SetColRange(18*6-1,18*7-1+2+1);
534     }
535     if ( mcmid%kNmcmRob == 3 ) {
536       fMCM->SetColRange(18*7-1,18*8-1+2+1);
537     }
538
539   } 
540
541   fMCM->SetRow(kNmcmRob*(robid/kNcolRob)+mcmid/kNmcmRob);
542
543 }
544
545 //_____________________________________________________________________________
546 void AliTRDtrigger::AddTracklet(Int_t det, Int_t row, Int_t seed, Int_t n)
547 {
548
549   Float_t field = fTrigParam->GetField();
550   AliTRDgeometry *geo = (AliTRDgeometry*)AliTRDgeometry::GetGeometry(fRunLoader);
551
552   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
553   if (!calibration)
554   {
555     Error("AddTracklets","No instance of AliTRDcalibDB.");
556     return;  
557   }
558   
559   Int_t nTimeTotal  = calibration->GetNumberOfTimeBins();
560
561   fTrk = new AliTRDmcmTracklet(det,row,n);
562
563   Int_t iCol, iCol1, iCol2, track[3];
564   iCol = fMCM->GetSeedCol()[seed];  // 0....20 (MCM)
565   fMCM->GetColRange(iCol1,iCol2);   // range in the pad plane
566             
567   Float_t amp[3];
568   for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) {
569
570     amp[0] = fMCM->GetADC(iCol-1,iTime);
571     amp[1] = fMCM->GetADC(iCol  ,iTime);
572     amp[2] = fMCM->GetADC(iCol+1,iTime);
573
574     // extract track contribution only from the central pad
575     track[0] = fTrack0->GetDataUnchecked(row,iCol+iCol1,iTime);
576     track[1] = fTrack1->GetDataUnchecked(row,iCol+iCol1,iTime);
577     track[2] = fTrack2->GetDataUnchecked(row,iCol+iCol1,iTime);
578
579     if (fMCM->IsCluster(iCol,iTime)) {
580
581       fTrk->AddCluster(iCol+iCol1,iTime,amp,track);
582
583     } else if ((iCol+1+1) < kMcmCol) {
584
585       amp[0] = fMCM->GetADC(iCol-1+1,iTime);
586       amp[1] = fMCM->GetADC(iCol  +1,iTime);
587       amp[2] = fMCM->GetADC(iCol+1+1,iTime);
588
589       if (fMCM->IsCluster(iCol+1,iTime)) {
590
591         // extract track contribution only from the central pad
592         track[0] = fTrack0->GetDataUnchecked(row,iCol+1+iCol1,iTime);
593         track[1] = fTrack1->GetDataUnchecked(row,iCol+1+iCol1,iTime);
594         track[2] = fTrack2->GetDataUnchecked(row,iCol+1+iCol1,iTime);
595
596         fTrk->AddCluster(iCol+1+iCol1,iTime,amp,track);
597
598       }
599
600     } else {
601     }
602
603   }
604
605   fTrk->CookLabel(0.8);  
606   /*
607   if (fTrk->GetLabel() >= fNPrimary) {
608     Info("AddTracklet","Only primaries are stored!");
609     return;
610   }
611   */
612   // LTU Pt cut
613   fTrk->MakeTrackletGraph(geo,field);
614   fTrk->MakeClusAmpGraph();
615   if (TMath::Abs(fTrk->GetPt()) < fTrigParam->GetLtuPtCut()) {
616     return;
617   }
618       
619   Tracklets()->Add(fTrk);
620
621 }
622
623 //_____________________________________________________________________________
624 Bool_t AliTRDtrigger::WriteTracklets(Int_t det) 
625 {
626   //
627   // Fills TRDmcmTracklet branch in the tree with the Tracklets 
628   // found in detector = det. For det=-1 writes the tree. 
629   //
630
631   if ((det < -1) || (det >= AliTRDgeometry::Ndet())) {
632     Error("WriteTracklets","Unexpected detector index %d.",det);
633     return kFALSE;
634   }
635
636   TBranch *branch = fTrackletTree->GetBranch("TRDmcmTracklet");
637   if (!branch) {
638     TObjArray *ioArray = 0;
639     branch = fTrackletTree->Branch("TRDmcmTracklet","TObjArray",&ioArray,32000,0);
640   }
641
642   if ((det >= 0) && (det < AliTRDgeometry::Ndet())) {
643
644     Int_t nTracklets = Tracklets()->GetEntriesFast();
645     TObjArray *detTracklets = new TObjArray(400);
646
647     for (Int_t i = 0; i < nTracklets; i++) {
648       AliTRDmcmTracklet *trk = (AliTRDmcmTracklet *) Tracklets()->UncheckedAt(i);
649       
650       if (det == trk->GetDetector()) {
651         detTracklets->AddLast(trk);
652       }
653       else {
654       }
655     }
656
657     branch->SetAddress(&detTracklets);
658     fTrackletTree->Fill();
659
660     delete detTracklets;
661
662     return kTRUE;
663
664   }
665
666   if (det == -1) {
667
668     Info("WriteTracklets","Writing the Tracklet tree %s for event %d."
669          ,fTrackletTree->GetName(),fRunLoader->GetEventNumber());
670
671     AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
672     loader->WriteTracks("OVERWRITE");
673     
674     return kTRUE;  
675
676   }
677
678   return kFALSE;
679
680 }
681
682 //_____________________________________________________________________________
683 void AliTRDtrigger::MakeTracks(Int_t det)
684 {
685   //
686   // Create GTU tracks per module (stack of 6 chambers)
687   //
688   
689   fModule->Reset();
690
691   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
692   if (!commonParam)
693   {
694     Error("MakeTracks","No common params.");
695     return;
696   }
697     
698   Int_t nRowMax, iplan, icham, isect, row;
699
700   AliTRDgeometry *geo = (AliTRDgeometry*)AliTRDgeometry::GetGeometry(fRunLoader);
701
702   if ((det < 0) || (det >= AliTRDgeometry::Ndet())) {
703     Error("MakeTracks","Unexpected detector index %d.",det);
704     return;
705   }
706   
707   Int_t nTracklets = Tracklets()->GetEntriesFast();
708   
709   AliTRDmcmTracklet *trk;
710   for (Int_t i = 0; i < nTracklets; i++) {
711     
712     trk = (AliTRDmcmTracklet *) Tracklets()->UncheckedAt(i);
713     
714     iplan = geo->GetPlane(trk->GetDetector());
715     icham = geo->GetChamber(trk->GetDetector());
716     isect = geo->GetSector(trk->GetDetector());
717
718     nRowMax = commonParam->GetRowMax(iplan,icham,isect);
719     row = trk->GetRow();
720
721     fModule->AddTracklet(trk->GetDetector(),
722                          row,
723                          trk->GetRowz(),
724                          trk->GetSlope(),
725                          trk->GetOffset(),
726                          trk->GetTime0(),
727                          trk->GetNclusters(),
728                          trk->GetLabel(),
729                          trk->GetdQdl());
730     
731   }
732
733   fModule->SortTracklets();
734   fModule->RemoveMultipleTracklets();
735   fModule->SortZ((Int_t)geo->GetChamber(det));
736   fModule->FindTracks();
737   fModule->SortTracks();
738   fModule->RemoveMultipleTracks();
739
740   Int_t nModTracks = fModule->GetNtracks();
741   AliTRDgtuTrack *gtutrk;
742   for (Int_t i = 0; i < nModTracks; i++) {
743     gtutrk = (AliTRDgtuTrack*)fModule->GetTrack(i);
744     if (TMath::Abs(gtutrk->GetPt()) < fTrigParam->GetGtuPtCut()) continue;
745     gtutrk->CookLabel();
746     gtutrk->MakePID();
747     AddTrack(gtutrk,det);
748   }
749   
750 }
751
752