HLT/trigger/AliHLTTRDTriggerComponent.cxx

   1 // $Id$
   2 //**************************************************************************
   3 //* This file is property of and copyright by the ALICE HLT Project        *
   4 //* ALICE Experiment at CERN, All rights reserved.                         *
   5 //*                                                                        *
   6 //* Primary Authors: Matthias Richter <Matthias.Richter@ift.uib.no>        *
   7 //*                  Jochen Thaeder <jochen@thaeder.de>                    *
   8 //*                  for The ALICE HLT Project.                            *
   9 //*                                                                        *
  10 //* Permission to use, copy, modify and distribute this software and its   *
  11 //* documentation strictly for non-commercial purposes is hereby granted   *
  12 //* without fee, provided that the above copyright notice appears in all   *
  13 //* copies and that both the copyright notice and this permission notice   *
  14 //* appear in the supporting documentation. The authors make no claims     *
  15 //* about the suitability of this software for any purpose. It is          *
  16 //* provided "as is" without express or implied warranty.                  *
  17 //**************************************************************************
  18
  19 /// @file   AliHLTTRDTriggerComponent.cxx
  20 /// @author Felix Rettig, Stefan Kirsch
  21 /// @date   2012-08-16
  22 /// @brief
  23
  24 #include <cstdlib>
  25 #include "TSystem.h"
  26 #include "TObjArray.h"
  27 #include "TObjString.h"
  28 #include "TFile.h"
  29 #include "TH1I.h"
  30 #include "TH2I.h"
  31 #include "AliESDtrack.h"
  32 #include "AliESDEvent.h"
  33 #include "AliHLTTRDDefinitions.h"
  34 #include "AliHLTTriggerDecision.h"
  35 #include "AliHLTDomainEntry.h"
  36 #include "AliHLTLogging.h"
  37 #include "AliHLTErrorGuard.h"
  38 #include "AliHLTGlobalBarrelTrack.h"
  39 #include "AliESDtrackCuts.h"
  40 //#include "AliHLTESDTrackCuts.h"
  41 #include "AliGeomManager.h"
  42 #include "AliHLTComponentBenchmark.h"
  43 #include "AliHLTTRDTriggerComponent.h"
  44 #include "AliTRDonlineTrackingDataContainer.h"
  45 #include "AliTRDpadPlane.h"
  46
  47 ClassImp(AliHLTTRDTriggerComponent)
  48
  49 #define LogError( ... ) { HLTError(__VA_ARGS__); if (fDebugLevel >= 1) { DbgLog("ERROR", __VA_ARGS__); } }
  50 #define LogInfo( ... ) { HLTInfo(__VA_ARGS__); if (fDebugLevel >= 1) { DbgLog("INFO", __VA_ARGS__); } }
  51 #define LogInspect( ... ) { HLTDebug(__VA_ARGS__); if (fDebugLevel >= 1) { DbgLog("INSPECT", __VA_ARGS__); } }
  52 #define LogDebug( ... ) { if (fDebugLevel >= 1) { HLTInfo(__VA_ARGS__); DbgLog("DEBUG", __VA_ARGS__); } }
  53
  54 AliHLTTRDTriggerComponent::AliHLTTRDTriggerComponent() :
  55   AliHLTTrigger(),
  56   fName(),
  57   fRefTrackSelectionEtaLimit(0.9),
  58   fRefTrackSelectionVertexXYLimit(20.),
  59   fRefTrackSelectionVertexZLimit(30.),
  60   fRefTrackSelectionPtThreshold(0.7),
  61   fMatchRatingThreshold(0.25),
  62   fElectronTriggerPtThresholdHSE(3.),
  63   fElectronTriggerPIDThresholdHSE(144),
  64   fElectronTriggerPtThresholdHQU(2.),
  65   fElectronTriggerPIDThresholdHQU(164),
  66   fApplyRefTrackCuts(kFALSE),
  67   fElectronTriggerOnL1TrgOnly(kFALSE),
  68   fHistoMode(1),
  69   fDebugLevel(0),
  70   fExtendedHistos(kFALSE),
  71   fEventRendering(kFALSE),
  72   fPushHistos(kFALSE),
  73   fWriteHistos(kFALSE),
  74   fEventId(fgkInvalidEventId),
  75   fRunNumber(-1),
  76   fChunkId(NULL),
  77   fSectorsWithData(0),
  78   fIsMinBiasEvent(kFALSE),
  79   fIsTRDElectronEvent(kFALSE),
  80   fESDtracksPresent(kFALSE),
  81   fHLTtracksPresent(kFALSE),
  82   fTRDGeometry(NULL),
  83   fEsdEvent(NULL),
  84   fTrackingData(NULL),
  85   fHLTTracks(NULL),
  86   fRefTrackCuts(NULL),
  87 #ifdef __TRDHLTDEBUG
  88   fEventDisplay(NULL),
  89   fBenchmark(NULL),
  90 #endif
  91   fHistArray(NULL),
  92   fHistMatchRating(NULL),
  93   fHistMatchRatingByPt(NULL),
  94   fHistMatchRatingByPid(NULL),
  95   fHistTrackPt(NULL),
  96   fHistTrackPtMatched(NULL),
  97   fHistTrackPtCorr(NULL),
  98   fHistTrackPid(NULL),
  99   fHistTrackPidMatched(NULL),
 100   fHistElectronCandidatePt(NULL),
 101   fHistElectronCandidateMatchedPt(NULL),
 102   fHistElectronCandidatePid(NULL),
 103   fHistElectronCandidateMatchedPid(NULL),
 104   fHistRefTrackPid(NULL),
 105   fHistMatchedRefTrackPid(NULL),
 106   fHistPIDvsTruncPID(NULL),
 107   fHistElectronFalsePIDvsTruncPID(NULL),
 108   fHistElectronConfirmedPIDvsTruncPID(NULL),
 109   fHistTrackMatchingCombinations(NULL),
 110   fHistElectronTriggerBaseMinBias(NULL),
 111   fHistElectronTriggerBaseTrdL1(NULL)
 112 {
 113 }
 114
 115 const char* AliHLTTRDTriggerComponent::fgkDefaultOCDBEntry = "HLT/ConfigHLT/TRDTrigger";
 116 const char* AliHLTTRDTriggerComponent::fgkTriggerDecisionElectronHSE = "TRD-ELECTRON-HSE";
 117 const char* AliHLTTRDTriggerComponent::fgkTriggerDecisionElectronHQU = "TRD-ELECTRON-HQU";
 118
 119 AliHLTTRDTriggerComponent::~AliHLTTRDTriggerComponent()
 120 {
 121 }
 122
 123 const char* AliHLTTRDTriggerComponent::GetTriggerName() const
 124 {
 125   if (!fName.IsNull())
 126     return fName.Data();
 127   else
 128     return "TRDTriggerComponent";
 129 }
 130
 131 AliHLTComponent* AliHLTTRDTriggerComponent::Spawn()
 132 {
 133   return new AliHLTTRDTriggerComponent;
 134 }
 135
 136 int AliHLTTRDTriggerComponent::DoInit(int argc, const char** argv)
 137 {
 138   int iResult = 0;
 139
 140   do {
 141
 142     fChunkId = new TString("XXXXX");
 143 #ifdef __TRDHLTDEBUG
 144     if (fChunkId){
 145       // chunk identification for debug output
 146       *fChunkId = gSystem->WorkingDirectory();
 147       fChunkId->Remove(0, fChunkId->Last('/') + 1);
 148       if (fChunkId->Contains("hlt_trd_comp"))
 149         *fChunkId = "L";
 150     } else {
 151       iResult = -ENOMEM;
 152       break;
 153     }
 154 #endif
 155
 156     AliGeomManager::LoadGeometry();
 157
 158     fTRDGeometry = new AliTRDgeometry();
 159     if (!fTRDGeometry) {
 160       iResult = -ENOMEM;
 161       break;
 162     }
 163
 164     fTrackingData = new AliTRDonlineTrackingDataContainer();
 165     if (!fTrackingData) {
 166       iResult = -ENOMEM;
 167       break;
 168     }
 169
 170     fHLTTracks = new vector<AliHLTGlobalBarrelTrack>;
 171     if (!fHLTTracks) {
 172       iResult = -ENOMEM;
 173       break;
 174     }
 175
 176     fRefTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "No track cuts");
 177     if (fRefTrackCuts){
 178       // fRefTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE, 0);
 179       fRefTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 180       fRefTrackCuts->SetEtaRange(-0.8, 0.8);
 181     } else {
 182       iResult = -ENOMEM;
 183       break;
 184     }
 185
 186     fHistArray = new TObjArray(25);
 187     if(!fHistArray)
 188       return -ENOMEM;
 189     fHistArray->SetOwner(kTRUE);
 190
 191     fHistMatchRating = new TH1I("trdtrg_match_rating", "Match rating distribution HLT tracks vs. GTU track;match rating;abundance;",
 192                                 1./0.02, 0., 1.);
 193     fHistArray->AddLast(fHistMatchRating);
 194
 195     fHistMatchRatingByPt = new TH2I("trdtrg_match_rating_pt", "Match rating distribution HLT tracks vs. GTU track by p_{T};match rating m;p_{T}^{TRD,online} (GeV/c)",
 196                                     1./0.02, 0., 1., 20./0.02, 0., 20.);
 197     fHistArray->AddLast(fHistMatchRatingByPt);
 198
 199     fHistMatchRatingByPid = new TH2I("trdtrg_match_rating_pid", "Match rating distribution HLT tracks vs. GTU track by PID;match rating m;PID (a.u.)",
 200                                      1./0.02, 0., 1., 256, 0, 256);
 201     fHistArray->AddLast(fHistMatchRatingByPid);
 202
 203     fHistTrackPt = new TH1I("trdtrg_track_pt", "GTU track p_{T};GTU track p_{T} (GeV/c);abundance", 200, -20, 20.);
 204     fHistArray->AddLast(fHistTrackPt);
 205
 206     fHistTrackPtMatched = new TH1I("trdtrg_matched_track_pt", "matched GTU track p_{T};GTU track p_{T} (GeV/c);abundance", 200, -20, 20.);
 207     fHistArray->AddLast(fHistTrackPtMatched);
 208
 209     fHistTrackPtCorr = new TH2I("trdtrg_track_pt_corr", "HLT vs. GTU track p_{T};p_{T}^{HLT} (GeV/c);p_{T}^{GTU} (GeV/c)",
 210                400, -40., 40., 400, -40., 40.);
 211     fHistArray->AddLast(fHistTrackPtCorr);
 212
 213     fHistTrackPid = new TH1I("trdtrg_track_pid", "GTU track PID;GTU track PID (a.u.);abundance", 256, 0, 256);
 214     fHistArray->AddLast(fHistTrackPid);
 215
 216     fHistTrackPidMatched = new TH1I("trdtrg_matched_track_pid", "matched GTU track PID;GTU track PID (a.u.);abundance", 256, 0, 256);
 217     fHistArray->AddLast(fHistTrackPidMatched);
 218
 219     fHistTrackMatchingCombinations = new TH2I("trdtrg_matching_combinations", "HLT-GTU track matching combinations;set;combinations to check;",
 220                                               10, 0, 10, 2000, 0, 10000);
 221     fHistArray->AddLast(fHistTrackMatchingCombinations);
 222
 223     fHistElectronTriggerBaseMinBias = new TH1I("trdtrg_electron_trigger_base_mb", "min. bias base numbers for electron trigger analysis;set;abundance;",
 224                                               10, 0, 10);
 225     fHistElectronTriggerBaseMinBias->GetXaxis()->SetBinLabel(1, "min. bias events");
 226     fHistElectronTriggerBaseMinBias->GetXaxis()->SetBinLabel(2, "TRD L1 electron triggered");
 227     fHistElectronTriggerBaseMinBias->GetXaxis()->SetBinLabel(3, "TRD HLT electron triggered");
 228     fHistArray->AddLast(fHistElectronTriggerBaseMinBias);
 229
 230     fHistElectronTriggerBaseTrdL1 = new TH1I("trdtrg_electron_trigger_base_l1", "TRD L1 base numbers for electron trigger analysis;set;abundance;",
 231                                               10, 0, 10);
 232     fHistElectronTriggerBaseTrdL1->GetXaxis()->SetBinLabel(1, "TRD L1 electron triggered");
 233     fHistElectronTriggerBaseTrdL1->GetXaxis()->SetBinLabel(2, "TRD HLT electron triggered");
 234     fHistArray->AddLast(fHistElectronTriggerBaseTrdL1);
 235
 236     fHistElectronCandidatePt = new TH1I("trdtrg_electron_candidate_pt", "GTU electron candidate p_{T};GTU track p_{T} (GeV/c);abundance", 200, -20, 20.);
 237     fHistArray->AddLast(fHistElectronCandidatePt);
 238
 239     fHistElectronCandidateMatchedPt = new TH1I("trdtrg_electron_candidate_matched_pt", "matching GTU electron candidate p_{T};GTU track p_{T} (GeV/c);abundance", 200, -20, 20.);
 240     fHistArray->AddLast(fHistElectronCandidateMatchedPt);
 241
 242     fHistElectronCandidatePid = new TH1I("trdtrg_electron_candidate_pid", "GTU electron candidate PID;GTU track PID (a.u.);abundance", 256, 0, 256);
 243     fHistArray->AddLast(fHistElectronCandidatePid);
 244
 245     fHistElectronCandidateMatchedPid = new TH1I("trdtrg_electron_candidate_matched_pid", "matching GTU electron candidate PID;GTU track PID (a.u.);abundance", 256, 0, 256);
 246     fHistArray->AddLast(fHistElectronCandidateMatchedPid);
 247
 248   } while (0);
 249
 250   if (iResult < 0){
 251
 252     if (fHistArray) delete fHistArray;
 253     fHistArray = NULL;
 254
 255     if (fTRDGeometry) delete fTRDGeometry;
 256     fTRDGeometry = NULL;
 257
 258     if (fRefTrackCuts) delete fRefTrackCuts;
 259     fRefTrackCuts = NULL;
 260
 261     if (fHLTTracks) delete fHLTTracks;
 262     fHLTTracks = NULL;
 263
 264     if (fTrackingData) delete fTrackingData;
 265     fTrackingData = NULL;
 266
 267     if (fChunkId) delete fChunkId;
 268     fChunkId = NULL;
 269
 270   }
 271
 272   // check if the -triggername argument is used, the trigger name determines the following initialization
 273   vector<const char*> remainingArgs;
 274   for (int i = 0; i < argc; ++i) {
 275     if (strcmp(argv[i], "-triggername") == 0) {
 276       if (++i < argc){
 277         fName = argv[i];
 278       } else {
 279         LogError("invalid parameter for argument '-triggername', string expected");
 280         return -EINVAL;
 281       }
 282       continue;
 283     }
 284     remainingArgs.push_back(argv[i]);
 285   }
 286
 287   TString cdbPath;
 288   if (!fName.IsNull()) {
 289     cdbPath = "HLT/ConfigHLT/";
 290     cdbPath += fName;
 291   } else {
 292     cdbPath = fgkDefaultOCDBEntry;
 293   }
 294
 295   LogInfo("cdbPath = <%s>", cdbPath.Data());
 296   iResult = ConfigureFromCDBObject(cdbPath);
 297
 298   if (iResult >= 0 && argc > 0)
 299     iResult = ConfigureFromArgumentString(remainingArgs.size(), &(remainingArgs[0]));
 300
 301   return iResult;
 302 }
 303
 304 int AliHLTTRDTriggerComponent::DoDeinit()
 305 {
 306
 307 #ifdef __TRDHLTDEBUG
 308   if (fEventDisplay) delete fEventDisplay;
 309   fEventDisplay = NULL;
 310 #endif
 311
 312   if ((fHistoMode == 1) && (fWriteHistos)){
 313     TFile out("trg_out/trg_hists.root", "RECREATE");
 314     if (!out.IsZombie()) {
 315       out.cd();
 316       UInt_t numHists = fHistArray->GetEntries();
 317       for (UInt_t iHist = 0; iHist < numHists; ++iHist)
 318         if (fHistArray->At(iHist))
 319           fHistArray->At(iHist)->Write();
 320       out.Close();
 321     }
 322   }
 323
 324   if (fHistArray) delete fHistArray;
 325   fHistArray = NULL;
 326
 327   if (fRefTrackCuts) delete fRefTrackCuts;
 328   fRefTrackCuts = NULL;
 329
 330   if (!fHLTTracks) delete fHLTTracks;
 331   fHLTTracks = NULL;
 332
 333   if (fTrackingData) delete fTrackingData;
 334   fTrackingData = NULL;
 335
 336   if (fTRDGeometry) delete fTRDGeometry;
 337   fTRDGeometry = NULL;
 338
 339   if (fChunkId) delete fChunkId;
 340   fChunkId = NULL;
 341
 342   return 0;
 343 }
 344
 345 int AliHLTTRDTriggerComponent::Reconfigure(const char* cdbEntry, const char* /*chainId*/)
 346 {
 347   TString cdbPath;
 348   if (!cdbEntry || cdbEntry[0] == 0) {
 349     if (!fName.IsNull()) {
 350       cdbPath = "HLT/ConfigHLT/";
 351       cdbPath += fName;
 352     } else {
 353       cdbPath = fgkDefaultOCDBEntry;
 354     }
 355   } else
 356     cdbPath = cdbEntry;
 357
 358   return ConfigureFromCDBObject(cdbPath);
 359 }
 360
 361 int AliHLTTRDTriggerComponent::ReadPreprocessorValues(const char* /*modules*/)
 362 {
 363   return 0;
 364 }
 365
 366 Int_t AliHLTTRDTriggerComponent::ConfigureFromCDBObject(TString /*cdbPath*/)
 367 {
 368   return 0;
 369 }
 370
 371 int AliHLTTRDTriggerComponent::ScanConfigurationArgument(int argc, const char** argv)
 372 {
 373
 374   if (argc <= 0)
 375     return 0;
 376
 377   UShort_t iArg = 0;
 378   TString argument(argv[iArg]);
 379
 380   if (!argument.CompareTo("-match-sel-eta")){
 381     if (++iArg >= argc) return -EPROTO;
 382     argument = argv[iArg];
 383     fRefTrackSelectionEtaLimit = argument.Atof();
 384     LogInfo("ref track selection eta limit set to %.1f.", fRefTrackSelectionEtaLimit);
 385     return 2;
 386   }
 387
 388   if (!argument.CompareTo("-match-sel-vxy")){
 389     if (++iArg >= argc) return -EPROTO;
 390     argument = argv[iArg];
 391     fRefTrackSelectionVertexXYLimit = argument.Atof();
 392     LogInfo("ref track selection vertex xy limit set to %.1f.", fRefTrackSelectionVertexXYLimit);
 393     return 2;
 394   }
 395
 396   if (!argument.CompareTo("-match-sel-vz")){
 397     if (++iArg >= argc) return -EPROTO;
 398     argument = argv[iArg];
 399     fRefTrackSelectionVertexZLimit = argument.Atof();
 400     LogInfo("ref track selection vertex z limit set to %.1f.", fRefTrackSelectionVertexZLimit);
 401     return 2;
 402   }
 403
 404   if (!argument.CompareTo("-match-sel-pt")){
 405     if (++iArg >= argc) return -EPROTO;
 406     argument = argv[iArg];
 407     fRefTrackSelectionPtThreshold = argument.Atof();
 408     LogInfo("ref track selection pt threshold set to %.1f GeV/c.", fRefTrackSelectionPtThreshold);
 409     return 2;
 410   }
 411
 412   if (!argument.CompareTo("-match-rating")){
 413     if (++iArg >= argc) return -EPROTO;
 414     argument = argv[iArg];
 415     fMatchRatingThreshold = argument.Atof();
 416     LogInfo("match rating threshold set to %.2f GeV/c.", fMatchRatingThreshold);
 417     return 2;
 418   }
 419
 420   if (!argument.CompareTo("-histo-mode")){
 421     if (++iArg >= argc) return -EPROTO;
 422     argument = argv[iArg];
 423     fHistoMode = argument.Atoi();
 424     LogInfo("histo mode set to %d", fHistoMode);
 425     return 2;
 426   }
 427
 428   if (!argument.CompareTo("-trghse")){
 429     LogInfo("TRD HLT electron trigger HSE enabled.");
 430     return 1;
 431   }
 432
 433   if (!argument.CompareTo("-hse-pt")){
 434     if (++iArg >= argc) return -EPROTO;
 435     argument = argv[iArg];
 436     fElectronTriggerPtThresholdHSE = argument.Atof();
 437     LogInfo("pt threshold for HSE trigger set to %.1f GeV/c.", fElectronTriggerPtThresholdHSE);
 438     return 2;
 439   }
 440
 441   if (!argument.CompareTo("-hse-pid")){
 442     if (++iArg >= argc) return -EPROTO;
 443     argument = argv[iArg];
 444     fElectronTriggerPIDThresholdHSE = argument.Atof();
 445     LogInfo("PID threshold for HSE trigger set to %d.", fElectronTriggerPIDThresholdHSE);
 446     return 2;
 447   }
 448
 449   if (!argument.CompareTo("-trghqu")){
 450     LogInfo("TRD HLT electron trigger HQU enabled.");
 451     return 1;
 452   }
 453
 454   if (!argument.CompareTo("-hqu-pt")){
 455     if (++iArg >= argc) return -EPROTO;
 456     argument = argv[iArg];
 457     fElectronTriggerPtThresholdHQU = argument.Atof();
 458     LogInfo("pt threshold for HQU trigger set to %.1f GeV/c.", fElectronTriggerPtThresholdHQU);
 459     return 2;
 460   }
 461
 462   if (!argument.CompareTo("-hqu-pid")){
 463     if (++iArg >= argc) return -EPROTO;
 464     argument = argv[iArg];
 465     fElectronTriggerPIDThresholdHQU = argument.Atof();
 466     LogInfo("PID threshold for HQU trigger set to %.1f GeV/c.", fElectronTriggerPIDThresholdHQU);
 467     return 2;
 468   }
 469
 470   if (!argument.CompareTo("-l1-only")){
 471     LogInfo("evaluation of electron trigger only for events with TRD L1 electron trigger enabled.");
 472     fElectronTriggerOnL1TrgOnly = kTRUE;
 473     return 1;
 474   }
 475
 476   if (!argument.CompareTo("-histo-ext")){
 477     LogInfo("extended histogramming enabled.");
 478     fExtendedHistos = kTRUE; // enable extended histogramming, for debugging/tuning only!
 479
 480     if (!fHistRefTrackPid){
 481       fHistRefTrackPid = new TH2I("trdtrg_reftrack_pid", "TPC dE/dx by p_{T} for matching GTU tracks;p_{T} (GeV/c);dE/dx (a. u.)",
 482                                   500, 0., 10., 500, 0., 500);
 483       fHistArray->AddLast(fHistRefTrackPid);
 484     }
 485
 486     if (!fHistMatchedRefTrackPid){
 487       fHistMatchedRefTrackPid = new TH2I("trdtrg_matched_reftrack_pid", "TPC dE/dx by p_{T} for matching GTU tracks;p_{T} (GeV/c);dE/dx (a. u.)",
 488                                          500, 0., 10., 500, 0., 500);
 489       fHistArray->AddLast(fHistMatchedRefTrackPid);
 490     }
 491
 492     if (!fHistPIDvsTruncPID){
 493       fHistPIDvsTruncPID = new TH2I("trdtrg_pid_trunc_pid", "GTU track PID vs. truncated PID;PID (a.u.);truncated PID (a. u.)",
 494                                     256, 0., 256., 256, 0., 256.);
 495       fHistArray->AddLast(fHistPIDvsTruncPID);
 496     }
 497
 498     if (!fHistElectronFalsePIDvsTruncPID){
 499       fHistElectronFalsePIDvsTruncPID = new TH2I("trdtrg_electron_false_pid_trunc_pid", "false electron PID vs. truncated PID;PID (a.u.);truncated PID (a. u.)",
 500                                                  256, 0., 256., 256, 0., 256.);
 501       fHistArray->AddLast(fHistElectronFalsePIDvsTruncPID);
 502     }
 503
 504     if (!fHistElectronConfirmedPIDvsTruncPID){
 505       fHistElectronConfirmedPIDvsTruncPID = new TH2I("trdtrg_electron_confirmed_pid_trunc_pid", "confirmed electron PID vs. truncated PID;PID (a.u.);truncated PID (a. u.)",
 506                                                      256, 0., 256., 256, 0., 256.);
 507       fHistArray->AddLast(fHistElectronConfirmedPIDvsTruncPID);
 508     }
 509
 510     return 1;
 511   }
 512
 513   if (!argument.CompareTo("-ref-cuts")){
 514     LogInfo("ref track cuts for matching enabled.");
 515     fApplyRefTrackCuts = kTRUE;
 516     return 1;
 517   }
 518
 519   if (!argument.CompareTo("-push-histograms")){
 520     LogInfo("pushing of histograms enabled.");
 521     fPushHistos = kTRUE; // enable histogram pushing, for debugging/tuning only!
 522     return 1;
 523   }
 524
 525   if (!argument.CompareTo("-write-histograms")){
 526     LogInfo("writing of histograms enabled.");
 527     fWriteHistos = kTRUE; // enable histogram writing, for debugging/tuning only!
 528     return 1;
 529   }
 530
 531   if (!argument.CompareTo("-render")){
 532 #ifdef __TRDHLTDEBUG
 533     LogInfo("rendering of interesting events enabled.");
 534     if (!fEventDisplay)
 535       fEventDisplay = new AliTRDtrackingEventDisplay();
 536       LogInfo("event rendering activated. this is for debugging only!");
 537     fEventRendering = kTRUE; // enable event rendering, for debugging/tuning only!
 538 #endif
 539     return 1;
 540   }
 541
 542   if (!argument.CompareTo("-debug")){
 543     if (++iArg >= argc) return -EPROTO;
 544     argument = argv[iArg];
 545     fDebugLevel = argument.Atoi();
 546     LogInfo("debug level set to %d.", fDebugLevel);
 547     return 2;
 548   }
 549
 550   if (!argument.CompareTo("-ref-pt")){
 551     if (++iArg >= argc) return -EPROTO;
 552     argument = argv[iArg];
 553     Float_t minPt, maxPt;
 554     fRefTrackCuts->GetPtRange(minPt, maxPt);
 555     maxPt = argument.Atof();
 556     fRefTrackCuts->SetPtRange(minPt, maxPt);
 557     LogInfo("ref track pt range set to %.1f .. %.1f GeV/c.", minPt, maxPt);
 558     return 2;
 559   }
 560
 561   if (!argument.CompareTo("-ref-tdca")){
 562     if (++iArg >= argc) return -EPROTO;
 563     argument = argv[iArg];
 564     fRefTrackCuts->SetMaxDCAToVertexXY(argument.Atof());
 565     LogInfo("ref track DCA transverse threshold set to %.3f", argument.Atof());
 566     return 2;
 567   }
 568
 569   if (!argument.CompareTo("-ref-ldca")){
 570     if (++iArg >= argc) return -EPROTO;
 571     argument = argv[iArg];
 572     fRefTrackCuts->SetMaxDCAToVertexZ(argument.Atof());
 573     LogInfo("ref track longitudinal DCA threshold set to %.3f", argument.Atof());
 574     return 2;
 575   }
 576
 577   return 1;
 578 }
 579
 580 void AliHLTTRDTriggerComponent::ScanTriggerClasses(const char* firedTriggerClasses) {
 581
 582   fIsMinBiasEvent = kFALSE;
 583
 584   TString trg(firedTriggerClasses);
 585   if (trg.Index("CINT7WU-S-NOPF-ALL") >= 0)
 586     fIsMinBiasEvent = kTRUE;
 587
 588   if (trg.Index("CINT8WU-S-NOPF-ALL") >= 0)
 589     fIsMinBiasEvent = kTRUE;
 590
 591 }
 592
 593 Bool_t AliHLTTRDTriggerComponent::CheckRefTrackCuts(AliESDtrack* track){
 594
 595   // cuts by cut class
 596   if (fApplyRefTrackCuts)
 597     return fRefTrackCuts->AcceptTrack(track) ? kTRUE : kFALSE;
 598
 599   // simple custom cuts
 600   Float_t dcaToVertexXY, dcaToVertexZ;
 601   track->GetImpactParametersTPC(dcaToVertexXY, dcaToVertexZ);
 602   // LogDebug("IMPACT TPC %.4f  %.4f", dcaToVertexXY, dcaToVertexZ);
 603
 604   if ((dcaToVertexXY > fRefTrackSelectionVertexXYLimit) || (dcaToVertexZ >= fRefTrackSelectionVertexZLimit))
 605     return kFALSE;
 606
 607   if (TMath::Abs(track->Eta()) > fRefTrackSelectionEtaLimit)
 608     return kFALSE;
 609
 610   return kTRUE;
 611 }
 612
 613
 614 void AliHLTTRDTriggerComponent::DbgLog(const char* prefix, ...){
 615 #ifdef __TRDHLTDEBUG
 616   AliHLTEventID_t eventNumber = fEventId;
 617   printf("TRDHLTGM %s-%s-%s: [TRG] %s",
 618          (fRunNumber >= 0) ? Form("%06d", fRunNumber) : "XXXXXX",
 619          fChunkId->Data(),
 620          (eventNumber != fgkInvalidEventId) ? Form("%05llu", eventNumber) : "XXXXX",
 621          (strlen(prefix) > 0) ? Form("<%s> ", prefix) : "");
 622 #endif
 623   va_list args;
 624   va_start(args, prefix);
 625   char* fmt = va_arg(args, char*);
 626   vprintf(fmt, args);
 627   printf("\n");
 628   va_end(args);
 629 }
 630
 631 int AliHLTTRDTriggerComponent::PrepareESDData(){
 632
 633   int result = 0;
 634   fESDtracksPresent = kFALSE;
 635
 636   // check access to ESD event data
 637   const TObject* obj = GetFirstInputObject(kAliHLTAllDataTypes, "AliESDEvent");
 638   fEsdEvent = dynamic_cast<AliESDEvent*>(const_cast<TObject*>(obj));
 639
 640   if (fEsdEvent) {
 641     fEsdEvent->GetStdContent();
 642     fRunNumber = fEsdEvent->GetRunNumber();
 643     fESDtracksPresent = kTRUE;
 644     unsigned int numTracks = fEsdEvent->GetNumberOfTracks();
 645
 646     // process trigger classes
 647     ScanTriggerClasses(fEsdEvent->GetFiredTriggerClasses().Data());
 648
 649     LogDebug("ESD event data: %d ESD tracks [event num: %d] [minbias: %d (fired:%s)]",
 650                     numTracks, fEsdEvent->GetEventNumberInFile(), fIsMinBiasEvent, fEsdEvent->GetFiredTriggerClasses().Data());
 651
 652     // process ESD tracks
 653     if ((fDebugLevel >= 2) || fExtendedHistos){
 654       AliESDtrack* esdTrack;
 655       TString paramStr("");
 656       for (unsigned int iTrack = 0; iTrack < numTracks; ++iTrack) {
 657         esdTrack = fEsdEvent->GetTrack(iTrack);
 658
 659         if (fExtendedHistos){
 660           fHistRefTrackPid->Fill(esdTrack->Pt(), esdTrack->GetTPCsignal());
 661         }
 662
 663         if (fDebugLevel >= 2){
 664           paramStr = "";
 665           if (esdTrack){
 666
 667             if (esdTrack->GetInnerParam())
 668               paramStr += "I";
 669
 670             if (esdTrack->GetOuterParam())
 671               paramStr += "O";
 672
 673             LogDebug("ESD track %4d - pt: %+.2fGeV/c  [params: S%s]", iTrack, esdTrack->GetSignedPt(), paramStr.Data());
 674           } else
 675             LogError("ESD data for track %d invalid", iTrack);
 676         }
 677
 678       } // loop over ESD tracks
 679     }
 680
 681     fTrackingData->SetGtuPtMultiplierFromMagField(fEsdEvent->GetMagneticField()); // used for sign correction
 682
 683     result = 1;
 684
 685   }
 686
 687   return result;
 688
 689 }
 690
 691 int AliHLTTRDTriggerComponent::PrepareHLTData(){
 692
 693   int iResult = 0;
 694   UInt_t numHLTTracks = 0;
 695   fHLTtracksPresent = kFALSE;
 696
 697   if (!fHLTTracks){
 698     LogError("HLT track vector instance not available.");
 699     return 0;
 700   }
 701
 702   fHLTTracks->clear();
 703
 704   for (const AliHLTComponentBlockData *pBlock = GetFirstInputBlock(kAliHLTDataTypeTrack | kAliHLTDataOriginTPC);
 705        pBlock != NULL; pBlock = GetNextInputBlock()) {
 706     LogDebug("#hlttrk - data block with HLT raw tracks received");
 707     fHLTtracksPresent = kTRUE;
 708
 709     // vector<AliHLTGlobalBarrelTrack> hltTracks;
 710     if ((iResult = AliHLTGlobalBarrelTrack::ConvertTrackDataArray(reinterpret_cast<const AliHLTTracksData*>(pBlock->fPtr), pBlock->fSize, *fHLTTracks)) >= 0) {
 711       for (vector<AliHLTGlobalBarrelTrack>::iterator element = fHLTTracks->begin();
 712            element != fHLTTracks->end(); element++) {
 713
 714         numHLTTracks++;
 715         //## AliHLTGlobalBarrelTrack -> AliKalmanTrack -> AliExternalTrackParam
 716
 717       } // loop over HLT tracks
 718     }
 719   } // loop over data blocks
 720
 721   LogDebug("#hlttrk - %d HLT raw tracks found", numHLTTracks);
 722
 723   return iResult;
 724
 725 }
 726
 727 int AliHLTTRDTriggerComponent::PrepareTRDData() {
 728
 729   int result = 1;
 730
 731   fSectorsWithData = 0;
 732   fTrackingData->Clear();
 733   for (const AliHLTComponentBlockData* datablock = GetFirstInputBlock(AliHLTTRDDefinitions::fgkOnlineDataType);
 734        datablock != NULL;
 735        datablock = GetNextInputBlock())
 736     {
 737       fSectorsWithData |= datablock->fSpecification;
 738       fTrackingData->Decompress(datablock->fPtr, datablock->fSize, kTRUE);
 739     }
 740
 741   fTrackingData->SetGtuPtMultiplierFromMagField(GetBz()); // used for sign correction
 742
 743   fTrackingData->PrintSummary("trigger component");
 744
 745   return result;
 746
 747 }
 748
 749 int AliHLTTRDTriggerComponent::MatchTRDTracksESD(){
 750
 751   if (!fEsdEvent) {
 752     LogError("ESD event data not available in MatchTRDTracks().");
 753     return 0;
 754   }
 755
 756   if (fHistoMode == 0){
 757     UInt_t numHists = fHistArray->GetEntries();
 758     for (UInt_t iHist = 0; iHist < numHists; ++iHist)
 759       if (fHistArray->At(iHist))
 760         ((TH1*) (fHistArray->At(iHist)))->Reset();
 761   }
 762
 763   int result = 1;
 764   unsigned int numRefTracks = fEsdEvent->GetNumberOfTracks();
 765   unsigned int numGtuTracks;
 766   Int_t refTrackIndices[fkMaxRefTracksPerStack];
 767   UInt_t refTrackCount = 0;
 768   Bool_t isRelevant;
 769   Double_t distY, distZ;
 770   Double_t matchRating;
 771   Double_t bestMatchRating;
 772   Int_t bestMatchRefIndex;
 773   AliESDtrack* refTrack = NULL;
 774   UInt_t numComparisonsDone = 0;
 775   UInt_t numUnmatchedTracks = 0;
 776   UInt_t numMatchedTracks = 0;
 777   Double_t magField = fEsdEvent->GetMagneticField();
 778
 779   for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack) {
 780     numGtuTracks = fTrackingData->GetNumTracks(iStack);
 781     refTrackCount = 0;
 782
 783     // preselect ESD relevant ESD tracks
 784     for (UInt_t iRefTrack = 0; iRefTrack < numRefTracks; ++iRefTrack) {
 785       refTrack = fEsdEvent->GetTrack(iRefTrack);
 786       isRelevant = (refTrack->Pt() >= fRefTrackSelectionPtThreshold); //## implement properly
 787       if (isRelevant){
 788         if (refTrackCount < fkMaxRefTracksPerStack)
 789           refTrackIndices[refTrackCount++] = iRefTrack;
 790         else {
 791           LogError("number of HLT tracks exceeding limit of %d. Skipping some tracks.", fkMaxRefTracksPerStack);
 792           break;
 793         }
 794       }
 795     }
 796
 797     // try to match GTU track with ref tracks
 798     for (UInt_t iGtuTrack = 0; iGtuTrack < numGtuTracks; ++iGtuTrack) {
 799       bestMatchRating = 0.;
 800       bestMatchRefIndex = -1;
 801       Double_t gpt = fTrackingData->GetTrackPt(iStack, iGtuTrack);
 802       UShort_t gpid = fTrackingData->GetTrackPID(iStack, iGtuTrack);
 803       UShort_t layerMask = fTrackingData->GetTrackLayerMask(iStack, iGtuTrack);
 804
 805       Float_t trklLocalY[fkTRDLayers];
 806       Int_t trklBinZ[fkTRDLayers];
 807       for (UShort_t iLayer = 0; iLayer < fkTRDLayers; ++iLayer){
 808         if ((layerMask >> iLayer) & 1){
 809           trklLocalY[iLayer] = fTrackingData->GetTrackTrackletLocalY(iStack, iGtuTrack, iLayer);
 810           trklBinZ[iLayer] = fTrackingData->GetTrackTrackletBinZ(iStack, iGtuTrack, iLayer);
 811         }
 812       }
 813
 814       for (UInt_t iRefTrack = 0; iRefTrack < refTrackCount; ++iRefTrack) {
 815         refTrack = fEsdEvent->GetTrack(refTrackIndices[iRefTrack]);
 816
 817 //      if (!CheckRefTrackCuts(refTrack))
 818 //        continue;
 819
 820         numComparisonsDone++;
 821
 822 //      if ((!refTrack->GetOuterParam()) && (!refTrack->GetInnerParam()))  // use tracks with TPC outer param only (all HLT tracks have it anyways)
 823 //        continue;
 824
 825         Int_t distRes = EstimateTrackDistance(refTrack, iStack, layerMask, trklLocalY, trklBinZ,
 826                                               magField, &distY, &distZ);
 827
 828         if (fDebugLevel >= 3){
 829           printf("CHECKMATCH = %i   distY %.2f   distZ %.2f   pt: %.2f   %.2f\n",
 830                 distRes, distY, distZ, gpt, refTrack->GetSignedPt());
 831         }
 832
 833         if (distRes == 0){
 834           matchRating = RateTrackMatch(distY, distZ, gpt, refTrack->GetSignedPt());
 835         } else {
 836           matchRating = 0.;
 837         }
 838
 839         if ((matchRating >= fMatchRatingThreshold) && (matchRating > bestMatchRating)) {
 840           bestMatchRefIndex = refTrackIndices[iRefTrack];
 841           bestMatchRating = matchRating;
 842         } else if (matchRating > bestMatchRating)
 843           bestMatchRating = matchRating;
 844
 845 //      DbgLog("", Form("#match - comparing GTU track %d in S%02d-%d with ref track %d: [gpt: %+5.2f rpt: %+5.2f] dy: %.1f dz: %.1f rating: %.1f",
 846 //                      iGtuTrack, iStack/5, iStack%5, refTrackIndices[iRefTrack],
 847 //                      fTrackingData->GetTrackPt(iStack, iGtuTrack), refTrack->GetSignedPt(),
 848 //                      distY, distZ, matchRating));
 849
 850       } // loop over ref tracks in stack
 851
 852       fHistMatchRating->Fill(bestMatchRating);
 853       fHistMatchRatingByPt->Fill(bestMatchRating, TMath::Abs(gpt));
 854       fHistMatchRatingByPid->Fill(bestMatchRating, gpid);
 855       fHistTrackPt->Fill(gpt);
 856       fHistTrackPid->Fill(gpid);
 857
 858       if (fExtendedHistos){
 859         fHistMatchedRefTrackPid->Fill(refTrack->Pt(), refTrack->GetTPCsignal());
 860       }
 861
 862       if (bestMatchRefIndex >= 0){
 863         // GTU track has matching reference track
 864         refTrack = fEsdEvent->GetTrack(bestMatchRefIndex);
 865         Double_t rpt = refTrack->GetSignedPt();
 866         LogDebug("#match - match found: rating %.2f, gpt: %+5.2f, rpt: %+5.2f (%i) -> diff: %.2f%%",
 867                             bestMatchRating,gpt, rpt, bestMatchRefIndex, (gpt - rpt)/rpt*100.);
 868         fTrackingData->SetTrackAddInfo(iStack, iGtuTrack, bestMatchRefIndex);
 869         if (fDebugLevel >= 3){
 870           LogDebug("#match-info rating: %.2f gpt: %+.2f  matchref: %d  rpt: %+.2f       gpid: %d  S%02d-%d %d   rpid: %.3f",
 871                           bestMatchRating, gpt, bestMatchRefIndex, rpt, gpid,
 872                           iStack/5, iStack%5, iGtuTrack, refTrack->GetTPCsignal());
 873         }
 874         numMatchedTracks++;
 875
 876         fHistTrackPtMatched->Fill(gpt);
 877         fHistTrackPtCorr->Fill(rpt, gpt);
 878         fHistTrackPidMatched->Fill(gpid);
 879
 880       } else {
 881         if (fDebugLevel >= 3){
 882           LogDebug("#match-info rating: %.2f gpt: %+.2f  no ref matching  gpid: %d  S%02d-%d %d",
 883                           bestMatchRating, gpt, gpid,
 884                           iStack/5, iStack%5, iGtuTrack);
 885         }
 886         numUnmatchedTracks++;
 887       }
 888
 889     } // loop over gtu tracks in stack
 890
 891   } // loop over stacks
 892
 893   LogInfo("#match - %d matched GTU tracks, %d unmatched",
 894                       numMatchedTracks, numUnmatchedTracks);
 895
 896   fHistTrackMatchingCombinations->Fill(0., (Double_t)(numRefTracks * fTrackingData->GetNumTracks()));  // index 0: full combinatorics
 897   fHistTrackMatchingCombinations->Fill(1., numComparisonsDone);  // index 1: track matching comparisons actually done
 898
 899   return result;
 900
 901 }
 902
 903 int AliHLTTRDTriggerComponent::MatchTRDTracksHLT(){
 904
 905   if (!fHLTTracks){
 906     LogError("HLT track data available.");
 907     return 0;
 908   }
 909
 910   Double_t magField = GetBz();
 911   unsigned int numGtuTracks;
 912   unsigned int numHLTTracks = 0;
 913   Bool_t isRelevant;
 914   Double_t distY, distZ;
 915   Double_t matchRating;
 916   Double_t bestMatchRating;
 917   Int_t bestMatchRefIndex;
 918   UInt_t numComparisonsDone = 0;
 919   UInt_t numUnmatchedTracks = 0;
 920   UInt_t numMatchedTracks = 0;
 921
 922   Bool_t hltTrackPreSel[fkMaxRefTracks];
 923   UInt_t iHltTrack;
 924
 925   for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack) {
 926     numGtuTracks = fTrackingData->GetNumTracks(iStack);
 927
 928     // preselect relevant HLT tracks
 929     iHltTrack = 0;
 930     for (vector<AliHLTGlobalBarrelTrack>::iterator element = fHLTTracks->begin();
 931          element != fHLTTracks->end(); element++) {
 932       numHLTTracks++;
 933       isRelevant = (element->Pt() >= fRefTrackSelectionPtThreshold); //## implement properly
 934
 935       //## use cuts here
 936       hltTrackPreSel[iHltTrack++] = isRelevant;
 937       if (iHltTrack >= fkMaxRefTracks) {
 938         LogError("maximum number of HLT tracks exceeded.");
 939         break;
 940       }
 941     } // loop over HLT tracks;
 942
 943     // search for matching HLT track for each GTU track
 944     for (UInt_t iGtuTrack = 0; iGtuTrack < numGtuTracks; ++iGtuTrack) {
 945       bestMatchRating = 0.;
 946       bestMatchRefIndex = -1;
 947       Double_t gpt = fTrackingData->GetTrackPt(iStack, iGtuTrack);
 948       UShort_t gpid = fTrackingData->GetTrackPID(iStack, iGtuTrack);
 949       UShort_t layerMask = fTrackingData->GetTrackLayerMask(iStack, iGtuTrack);
 950
 951       Float_t trklLocalY[fkTRDLayers];
 952       Int_t trklBinZ[fkTRDLayers];
 953       for (UShort_t iLayer = 0; iLayer < fkTRDLayers; ++iLayer){
 954         if ((layerMask >> iLayer) & 1){
 955           trklLocalY[iLayer] = fTrackingData->GetTrackTrackletLocalY(iStack, iGtuTrack, iLayer);
 956           trklBinZ[iLayer] = fTrackingData->GetTrackTrackletBinZ(iStack, iGtuTrack, iLayer);
 957         }
 958       }
 959
 960       iHltTrack = 0;
 961       for (vector<AliHLTGlobalBarrelTrack>::iterator element = fHLTTracks->begin();
 962            element != fHLTTracks->end(); element++) {
 963         if (!hltTrackPreSel[iHltTrack]){
 964           iHltTrack++;
 965           continue;
 966         }
 967
 968         // compare GTU track and relevant HLT track
 969         numComparisonsDone++;
 970
 971         AliHLTGlobalBarrelTrack hltPar(*element);
 972         Int_t distRes = EstimateTrackDistance(&hltPar, iStack, layerMask, trklLocalY, trklBinZ,
 973                                                magField, &distY, &distZ);
 974
 975         if (fDebugLevel >= 3){
 976           printf("CHECKMATCH = %i   distY %.2f   distZ %.2f   pt: %.2f   %.2f\n",
 977                  distRes, distY, distZ, gpt, element->GetSignedPt());
 978         }
 979
 980         if (distRes == 0){
 981           matchRating = RateTrackMatch(distY, distZ, gpt, element->GetSignedPt());
 982         } else {
 983           matchRating = 0.;
 984         }
 985
 986         if ((matchRating >= fMatchRatingThreshold) && (matchRating > bestMatchRating)) {
 987           bestMatchRefIndex = iHltTrack;
 988           bestMatchRating = matchRating;
 989         } else if (matchRating > bestMatchRating)
 990           bestMatchRating = matchRating;
 991
 992
 993         iHltTrack++;
 994       } // loop over HLT tracks;
 995
 996       fHistMatchRating->Fill(bestMatchRating);
 997       fHistMatchRatingByPt->Fill(bestMatchRating, TMath::Abs(gpt));
 998       fHistMatchRatingByPid->Fill(bestMatchRating, gpid);
 999       fHistTrackPt->Fill(gpt);
1000       fHistTrackPid->Fill(gpid);
1001
1002       if (bestMatchRefIndex >= 0){
1003         // GTU track has matching reference track
1004         Double_t rpt = fHLTTracks->at(bestMatchRefIndex).GetSignedPt();
1005         LogDebug("#match - match found: rating %.2f, gpt: %+5.2f, rpt: %+5.2f (%i) -> diff: %.2f%%",
1006                             bestMatchRating,gpt, rpt, bestMatchRefIndex, (gpt - rpt)/rpt*100.);
1007         fTrackingData->SetTrackAddInfo(iStack, iGtuTrack, bestMatchRefIndex);
1008
1009 //      if (fExtendedHistos){
1010 //        fHistMatchedRefTrackPid->Fill(element->Pt(), element->GetTPCsignal());
1011 //      }
1012 //
1013 //      if (fDebugLevel >= 3){
1014 //        LogDebug("#match-info rating: %.2f gpt: %+.2f  matchref: %d  rpt: %+.2f       gpid: %d  S%02d-%d %d   rpid: %.3f",
1015 //                        bestMatchRating, gpt, bestMatchRefIndex, rpt, gpid,
1016 //                        iStack/5, iStack%5, iGtuTrack, refTrack->GetTPCsignal());
1017 //      }
1018         numMatchedTracks++;
1019
1020         fHistTrackPtMatched->Fill(gpt);
1021         fHistTrackPtCorr->Fill(rpt, gpt);
1022         fHistTrackPidMatched->Fill(gpid);
1023
1024       } else {
1025         if (fDebugLevel >= 3){
1026           LogDebug("#match-info rating: %.2f gpt: %+.2f  no ref matching  gpid: %d  S%02d-%d %d",
1027                           bestMatchRating, gpt, gpid,
1028                           iStack/5, iStack%5, iGtuTrack);
1029         }
1030         numUnmatchedTracks++;
1031       }
1032
1033     } // loop over gtu tracks
1034   } // loop over stacks
1035
1036   LogInfo("#match - %d matched GTU tracks, %d unmatched (%d comparisons)",
1037           numMatchedTracks, numUnmatchedTracks, numComparisonsDone);
1038
1039   fHistTrackMatchingCombinations->Fill(0., (Double_t)(numHLTTracks * fTrackingData->GetNumTracks()));  // index 0: full combinatorics
1040   fHistTrackMatchingCombinations->Fill(1., numComparisonsDone);  // index 1: track matching comparisons actually done
1041
1042   return kTRUE;
1043 }
1044
1045 int AliHLTTRDTriggerComponent::MatchTRDTracks(){
1046
1047   if (!fEsdEvent) {
1048     LogError("ESD event data not available in MatchTRDTracks().");
1049     return 0;
1050   }
1051
1052   if (fHistoMode == 0){
1053     UInt_t numHists = fHistArray->GetEntries();
1054     for (UInt_t iHist = 0; iHist < numHists; ++iHist)
1055       if (fHistArray->At(iHist))
1056         ((TH1*) (fHistArray->At(iHist)))->Reset();
1057   }
1058
1059   int result = 1;
1060   unsigned int numRefTracks = fEsdEvent->GetNumberOfTracks();
1061   unsigned int numGtuTracks;
1062   Int_t refTrackIndices[fkMaxRefTracksPerStack];
1063   UInt_t refTrackCount = 0;
1064   Bool_t isRelevant;
1065   Double_t distY, distZ;
1066   Double_t matchRating;
1067   Double_t bestMatchRating;
1068   Int_t bestMatchRefIndex;
1069   AliESDtrack* refTrack;
1070   UInt_t numComparisonsDone = 0;
1071   UInt_t numUnmatchedTracks = 0;
1072   UInt_t numMatchedTracks = 0;
1073
1074   for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack) {
1075     numGtuTracks = fTrackingData->GetNumTracks(iStack);
1076     refTrackCount = 0;
1077
1078     // preselect ESD relevant ESD tracks
1079     for (UInt_t iRefTrack = 0; iRefTrack < numRefTracks; ++iRefTrack) {
1080       refTrack = fEsdEvent->GetTrack(iRefTrack);
1081       isRelevant = (refTrack->Pt() >= fRefTrackSelectionPtThreshold);
1082       if (isRelevant){
1083         if (refTrackCount < fkMaxRefTracksPerStack)
1084           refTrackIndices[refTrackCount++] = iRefTrack;
1085         else {
1086           LogError("number of HLT tracks exceeding limit of %d. Skipping some tracks.", fkMaxRefTracksPerStack);
1087           break;
1088         }
1089       }
1090     }
1091
1092     // try to match GTU track with ref tracks
1093     for (UInt_t iGtuTrack = 0; iGtuTrack < numGtuTracks; ++iGtuTrack) {
1094       bestMatchRating = 0.;
1095       bestMatchRefIndex = -1;
1096       Double_t gpt = fTrackingData->GetTrackPt(iStack, iGtuTrack);
1097       UShort_t gpid = fTrackingData->GetTrackPID(iStack, iGtuTrack);
1098       UShort_t layerMask = fTrackingData->GetTrackLayerMask(iStack, iGtuTrack);
1099
1100       for (UInt_t iRefTrack = 0; iRefTrack < refTrackCount; ++iRefTrack) {
1101         refTrack = fEsdEvent->GetTrack(refTrackIndices[iRefTrack]);
1102
1103         if (!CheckRefTrackCuts(refTrack))
1104           continue;
1105
1106         numComparisonsDone++;
1107
1108 //      if ((!refTrack->GetOuterParam()) && (!refTrack->GetInnerParam()))  // use tracks with TPC outer param only (all HLT tracks have it anyways)
1109 //        continue;
1110
1111         Float_t trklLocalY[fkTRDLayers];
1112         Int_t trklBinZ[fkTRDLayers];
1113         for (UShort_t iLayer = 0; iLayer < fkTRDLayers; ++iLayer){
1114           if ((layerMask >> iLayer) & 1){
1115             trklLocalY[iLayer] = fTrackingData->GetTrackTrackletLocalY(iStack, iGtuTrack, iLayer);
1116             trklBinZ[iLayer] = fTrackingData->GetTrackTrackletBinZ(iStack, iGtuTrack, iLayer);
1117           }
1118         }
1119
1120         Int_t distRes = EstimateTrackDistance(refTrack, iStack, layerMask, trklLocalY, trklBinZ,
1121                                               fEsdEvent->GetMagneticField(), &distY, &distZ);
1122
1123         if (distRes == 0){
1124           matchRating = RateTrackMatch(distY, distZ, gpt, refTrack->GetSignedPt());
1125         } else {
1126           matchRating = 0.;
1127         }
1128
1129         if ((matchRating >= fMatchRatingThreshold) && (matchRating > bestMatchRating)) {
1130           bestMatchRefIndex = refTrackIndices[iRefTrack];
1131           bestMatchRating = matchRating;
1132         } else if (matchRating > bestMatchRating)
1133           bestMatchRating = matchRating;
1134
1135 //      DbgLog("", Form("#match - comparing GTU track %d in S%02d-%d with ref track %d: [gpt: %+5.2f rpt: %+5.2f] dy: %.1f dz: %.1f rating: %.1f",
1136 //                      iGtuTrack, iStack/5, iStack%5, refTrackIndices[iRefTrack],
1137 //                      fTrackingData->GetTrackPt(iStack, iGtuTrack), refTrack->GetSignedPt(),
1138 //                      distY, distZ, matchRating));
1139
1140       } // loop over ref tracks in stack
1141
1142       fHistMatchRating->Fill(bestMatchRating);
1143       fHistMatchRatingByPt->Fill(bestMatchRating, TMath::Abs(gpt));
1144       fHistMatchRatingByPid->Fill(bestMatchRating, gpid);
1145       fHistTrackPt->Fill(gpt);
1146       fHistTrackPid->Fill(gpid);
1147
1148       if (bestMatchRefIndex >= 0){
1149         // GTU track has matching reference track
1150         refTrack = fEsdEvent->GetTrack(bestMatchRefIndex);
1151         Double_t rpt = refTrack->GetSignedPt();
1152         LogDebug("#match - match found: rating %.2f, gpt: %+5.2f, rpt: %+5.2f (%i) -> diff: %.2f%%",
1153                             bestMatchRating,gpt, rpt, bestMatchRefIndex, (gpt - rpt)/rpt*100.);
1154         fTrackingData->SetTrackAddInfo(iStack, iGtuTrack, bestMatchRefIndex);
1155
1156         if (fExtendedHistos){
1157           fHistMatchedRefTrackPid->Fill(refTrack->Pt(), refTrack->GetTPCsignal());
1158         }
1159
1160         if (fDebugLevel >= 3){
1161           LogDebug("#match-info rating: %.2f gpt: %+.2f  matchref: %d  rpt: %+.2f       gpid: %d  S%02d-%d %d   rpid: %.3f",
1162                           bestMatchRating, gpt, bestMatchRefIndex, rpt, gpid,
1163                           iStack/5, iStack%5, iGtuTrack, refTrack->GetTPCsignal());
1164         }
1165         numMatchedTracks++;
1166
1167         fHistTrackPtMatched->Fill(gpt);
1168         fHistTrackPtCorr->Fill(rpt, gpt);
1169         fHistTrackPidMatched->Fill(gpid);
1170
1171       } else {
1172         if (fDebugLevel >= 3){
1173           LogDebug("#match-info rating: %.2f gpt: %+.2f  no ref matching  gpid: %d  S%02d-%d %d",
1174                           bestMatchRating, gpt, gpid,
1175                           iStack/5, iStack%5, iGtuTrack);
1176         }
1177         numUnmatchedTracks++;
1178       }
1179
1180     } // loop over gtu tracks in stack
1181
1182   } // loop over stacks
1183
1184   LogInfo("#match - %d matched GTU tracks, %d unmatched",
1185                       numMatchedTracks, numUnmatchedTracks);
1186
1187   fHistTrackMatchingCombinations->Fill(0., (Double_t)(numRefTracks * fTrackingData->GetNumTracks()));  // index 0: full combinatorics
1188   fHistTrackMatchingCombinations->Fill(1., numComparisonsDone);  // index 1: track matching comparisons actually done
1189
1190   return result;
1191 }
1192
1193 void AliHLTTRDTriggerComponent::DumpTrackingData(){
1194
1195   if (fTrackingData->GetNumTracklets() + fTrackingData->GetNumTracks() == 0)
1196     return;
1197
1198   TString trklStr("");
1199   TString matchStr("");
1200   UShort_t layerMask;
1201
1202 //  for (UShort_t iSector = 0; iSector < 18; ++iSector){
1203 //    if (fTrackingData->GetSectorTrgWord(iSector) != ((UInt_t)0x2345352 ^ ((UInt_t)iSector + 34)))
1204 //      LogError("invalid sector trigger word in sector %02d: trg word is 0x%08x, should be 0x%08x",
1205 //             iSector, fTrackingData->GetSectorTrgWord(iSector), ((UInt_t)0x2345352 ^ ((UInt_t)iSector + 34)));
1206 //    for (UShort_t iStack = 0; iStack < 5; ++iStack){
1207 //      ULong64_t dwl = ((ULong64_t)0xaffe << 16) | (ULong64_t)iSector << 8 | iStack;
1208 //      ULong64_t dw = (((ULong64_t)0xbead << 16) | (ULong64_t)iSector << 8 | iStack) | ((ULong64_t)dwl << 32);
1209 //      if (fTrackingData->GetStackTrgWord(iSector, iStack) != dw)
1210 //      LogError("stack %02d-%d trg word is 0x%016llx, should be 0x%016llx",
1211 //               iSector, iStack, fTrackingData->GetStackTrgWord(iSector, iStack), dw);
1212 //    }
1213 //  }
1214
1215   for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack){
1216     for (Int_t iTrk = 0; iTrk < fTrackingData->GetNumTracks(iStack); ++iTrk){
1217
1218       layerMask = fTrackingData->GetTrackLayerMask(iStack, iTrk);
1219
1220       trklStr = Form("trkl: ");
1221       for (Short_t iLayer = 5; iLayer >= 0; --iLayer){
1222         if ((layerMask >> iLayer) & 1)
1223           trklStr += Form("0x%08x (%+8.3f)  ",
1224                           fTrackingData->GetTrackTrackletWord(iStack, iTrk, iLayer),
1225                           fTrackingData->GetTrackTrackletLocalY(iStack, iTrk, iLayer));
1226         else
1227           trklStr += "---------------------  ";
1228       } // loop over layers
1229       trklStr.Remove(trklStr.Length() - 2, 2);
1230
1231       if (fTrackingData->GetTrackAddInfo(iStack, iTrk) >= 0){
1232         Double_t rpt = (fESDtracksPresent) ?  fEsdEvent->GetTrack(fTrackingData->GetTrackAddInfo(iStack, iTrk))->GetSignedPt() :
1233           fHLTTracks->at(fTrackingData->GetTrackAddInfo(iStack, iTrk)).GetSignedPt();
1234         matchStr = Form("mpt: %+7.2f", rpt);
1235       } else
1236         matchStr = "unmatched";
1237
1238       if (fDebugLevel >= 3){
1239
1240         printf("###DOTDA EV%04llu  GTU TRACK - S%02d-%d  pt: %+7.2f  pid: %3d  lm: 0x%02x %s  %s\n",
1241                fEventId,
1242                iStack/5, iStack%5, fTrackingData->GetTrackPt(iStack, iTrk), fTrackingData->GetTrackPID(iStack, iTrk),
1243                layerMask, trklStr.Data(),
1244                matchStr.Data());
1245
1246         printf("###DOTDB EV%04llu  GTU TRACK - S%02d-%d  pt: %+7.2f  pid: %3d  lm: 0x%02x %s\n",
1247                fEventId,
1248                iStack/5, iStack%5, fTrackingData->GetTrackPt(iStack, iTrk), fTrackingData->GetTrackPID(iStack, iTrk),
1249                layerMask, trklStr.Data());
1250       }
1251
1252
1253       // paranoia checks
1254       for (Short_t iLayer = 5; iLayer >= 0; --iLayer){
1255         if (((layerMask >> iLayer) & 1) && (fTrackingData->GetTrackTrackletWord(iStack, iTrk, iLayer) == 0x10001000))
1256           LogError("invalid layer mask / tracklet value combination A");
1257
1258         if ((((layerMask >> iLayer) & 1) == 0) && (fTrackingData->GetTrackTrackletWord(iStack, iTrk, iLayer) != 0x10001000))
1259           LogError("invalid layer mask / tracklet value combination B");
1260       }
1261
1262     } // loop over tracks in stack
1263   } // loop over stacks
1264
1265 }
1266
1267 void AliHLTTRDTriggerComponent::AssignTrackInfo(TString* infoStr, const UInt_t stack, const UInt_t trackIndex, const char* flagStr) {
1268
1269   TString flags(flagStr);
1270   Bool_t appendRefData = kFALSE;
1271
1272   if (fTrackingData->GetTrackAddInfo(stack, trackIndex) >= 0){
1273     appendRefData = kTRUE;
1274     if (flags.First('M') < 0)
1275       flags += "M";
1276   }
1277
1278   *infoStr = Form("EXCH-TRK-INFO DS<%s> CH<%s> EV%05llu SEC%02d-%d-%d  gpt: %+6.1f  gpid: %3d  flags: %s",
1279                   "TRDHLTTRG",
1280                   fChunkId->Data(), fEventId,
1281                   stack/5, stack%5, trackIndex,
1282                   fTrackingData->GetTrackPt(stack, trackIndex),
1283                   fTrackingData->GetTrackPID(stack, trackIndex), flags.Data());
1284
1285   if (appendRefData){
1286     Double_t rpt = (fESDtracksPresent) ? fEsdEvent->GetTrack(fTrackingData->GetTrackAddInfo(stack, trackIndex))->GetSignedPt() :
1287       fHLTTracks->at(fTrackingData->GetTrackAddInfo(stack, trackIndex)).GetSignedPt();
1288     *infoStr += Form("  rpt: %+6.1f", rpt);
1289   }
1290
1291 }
1292
1293 #ifdef __TRDHLTDEBUG
1294 void AliHLTTRDTriggerComponent::RenderEvent(const Bool_t showGtuTracks, const Bool_t showTracklets, const Bool_t showRefTracks) {
1295
1296   if ((!fEventDisplay) || (!fEsdEvent))
1297     return;
1298
1299   LogDebug("rendering event");
1300
1301   const Float_t refTrackPtDisplayThreshold = 0.7;
1302   const Float_t trackPtEmphasizeThreshold = 1.8;
1303
1304   fEventDisplay->Reset();
1305   fEventDisplay->SetMagField(fEsdEvent->GetMagneticField());
1306   for (UInt_t iDet = 0; iDet < fkTRDStacks*fkTRDLayers; ++iDet)
1307     fEventDisplay->SetChamberState(iDet/fkTRDLayers, iDet%6, ((fSectorsWithData >> (iDet/30)) & 1) );
1308
1309   if (showTracklets){
1310     for (UShort_t iDet = 0; iDet < fkTRDStacks*fkTRDLayers; ++iDet){
1311       UInt_t numTrkl = fTrackingData->GetNumTracklets(iDet);
1312       for (UInt_t iTrkl = 0; iTrkl < numTrkl; ++iTrkl){
1313         fEventDisplay->AddTracklet(fTrackingData->GetTracklet(iDet, iTrkl));
1314       }
1315     }
1316   }
1317
1318   if (showGtuTracks){
1319     for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack){
1320       for (Int_t iTrk = 0; iTrk < fTrackingData->GetNumTracks(iStack); ++iTrk){
1321         AliESDTrdTrack* track = fTrackingData->GetTrack(iStack, iTrk);
1322         Int_t trkIndex = fEventDisplay->AddTrack(track, (TMath::Abs(track->Pt()) >= trackPtEmphasizeThreshold) ? (kMagenta + 2) : kMagenta, 1, 1);
1323         if (fTrackingData->GetTrackAddInfo(iStack, iTrk) >= 0)
1324           fEventDisplay->SetTrackTrackletStyle(trkIndex, kRed, 12);
1325         else
1326           fEventDisplay->SetTrackTrackletStyle(trkIndex, kViolet - 5, 12);
1327       }
1328     }
1329   }
1330
1331   unsigned int numRefTracks = fEsdEvent->GetNumberOfTracks();
1332   unsigned int numRefTracksRendered = 0;
1333   if (showRefTracks){
1334     // check for some marks for rendering
1335     UShort_t marks[40000];
1336     memset(marks, 0, sizeof(UShort_t)*40000);
1337     for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack){
1338       for (Int_t iTrk = 0; iTrk < fTrackingData->GetNumTracks(iStack); ++iTrk){
1339         if (fTrackingData->GetTrackAddInfo(iStack, iTrk) >= 0){
1340           marks[fTrackingData->GetTrackAddInfo(iStack, iTrk)] |= 1;
1341           if (
1342               (fTrackingData->GetTrackPID(iStack, iTrk) >= fElectronTriggerPIDThresholdHQU) &&
1343               (TMath::Abs(fTrackingData->GetTrackPt(iStack, iTrk)) >= fElectronTriggerPtThresholdHQU))
1344             marks[fTrackingData->GetTrackAddInfo(iStack, iTrk)] |= 2;
1345         }
1346       }
1347     }
1348
1349     // add to rendering
1350     for (unsigned int iTrack = 0; iTrack < numRefTracks; ++iTrack){
1351       AliESDtrack* track = fEsdEvent->GetTrack(iTrack);
1352       if ((track->Pt() >= refTrackPtDisplayThreshold) || (marks[iTrack] != 0)){
1353         Color_t color = (track->Pt() >= trackPtEmphasizeThreshold) ? (kGray + 2) : kGray;
1354         UShort_t width = 1;
1355         UShort_t style = 1;
1356         if (marks[iTrack] & 0x2){
1357           color = kRed + 1;
1358           width = 6;
1359         }
1360         if (!CheckRefTrackCuts(track))
1361           style = 2;
1362
1363         fEventDisplay->AddTrack(track, color, width, style);
1364         numRefTracksRendered++;
1365       }
1366     }
1367   }
1368
1369   if (!fEventDisplay->IsEmpty()){
1370     fEventDisplay->SetTitle("");
1371     fEventDisplay->SetSetupText("HLT", Form("%.1f#scale[0.5]{ }/#scale[0.5]{ }%.1f",
1372                                           refTrackPtDisplayThreshold, trackPtEmphasizeThreshold));
1373     fEventDisplay->SetBottomText(Form("%05i-%s-%05llu",
1374                                       fRunNumber, fChunkId->Data(), fEventId));
1375     fEventDisplay->SetBottomTextRight(Form("%d (%d) HLT tracks, %d tracklets, %d GTU tracks",
1376                                            numRefTracks,
1377                                            numRefTracksRendered,
1378                                            fTrackingData->GetNumTracklets(),
1379                                            fTrackingData->GetNumTracks()));
1380     fEventDisplay->SetLook(AliTRDtrackingEventDisplay::dmMediumLight);
1381     fEventDisplay->SetDisplayMode(AliTRDtrackingEventDisplay::dmFullXY);
1382     fEventDisplay->SaveToFile(Form("display/event-%s-%05llu.eps",
1383                                    fChunkId->Data(), fEventId));
1384   }
1385
1386 }
1387 #endif
1388
1389 Bool_t AliHLTTRDTriggerComponent::TRDElectronTrigger(const char *ident, const Double_t minPt, const UShort_t minPID){
1390
1391   LogDebug("Electron trigger processing (%s: pt>=%.1f, pid>=%d)...", ident, minPt, minPID);
1392
1393   UInt_t numTracks;
1394   Bool_t highPtElectronSeenGTU = kFALSE;
1395   Bool_t highPtElectronSeen = kFALSE;
1396   UInt_t truncMeanPID = 0;
1397   TString trackExchangeInfo("");
1398   TString flags("");
1399
1400   UInt_t trdSectorTrgContribs = 0;
1401   for (UShort_t iSector = 0; iSector < fkTRDSectors; ++iSector)
1402     trdSectorTrgContribs |= fTrackingData->GetSectorTrgContribs(iSector);
1403   if ((trdSectorTrgContribs >> 5) & 0x3)
1404     fIsTRDElectronEvent = kTRUE;
1405   else
1406     fIsTRDElectronEvent = kFALSE;
1407
1408   if (fIsMinBiasEvent)
1409     fHistElectronTriggerBaseMinBias->Fill(0.);
1410
1411   if (fIsTRDElectronEvent)
1412     fHistElectronTriggerBaseTrdL1->Fill(0.);
1413
1414   if ((fElectronTriggerOnL1TrgOnly) && (!fIsTRDElectronEvent)){
1415     // evaluate trigger for events with TRD L1 electron trigger fired
1416     DbgLog("skipping %s electron trigger evaluation for event, because no TRD trigger flag set (ctbs: 0x%02x)",
1417            ident, trdSectorTrgContribs);
1418     return 0;
1419   }
1420
1421   for (UShort_t iStack = 0; iStack < fkTRDStacks; ++iStack){
1422     numTracks = fTrackingData->GetNumTracks(iStack);
1423     for (UInt_t iTrk = 0; iTrk < numTracks; ++iTrk){
1424       Double_t gpt = fTrackingData->GetTrackPt(iStack, iTrk);
1425       Bool_t trdElectronCandidate = kFALSE;
1426       Bool_t hltElectronCandidate = kFALSE;
1427
1428       // re-evaluate GTU only decision for comparison
1429       if (
1430           (TMath::Abs(gpt) >= minPt) &&
1431           (fTrackingData->GetTrackPID(iStack, iTrk) >= minPID)
1432           ) {
1433         trdElectronCandidate = kTRUE;
1434         highPtElectronSeenGTU = kTRUE;
1435         fHistElectronCandidatePt->Fill(gpt);
1436         fHistElectronCandidatePid->Fill(fTrackingData->GetTrackPID(iStack, iTrk));
1437
1438         if (fExtendedHistos){
1439
1440           // idea to be checked
1441           truncMeanPID = 0;
1442           Short_t minLyr = -1;
1443           UShort_t minValue = 255;
1444           Short_t maxLyr = -1;
1445           UShort_t maxValue = 0;
1446           UShort_t lyrCount = 0;
1447           UInt_t layerMask = fTrackingData->GetTrackLayerMask(iStack, iTrk);
1448
1449           // scan for min & max values
1450           for (UShort_t iLayer = 0; iLayer < fkTRDLayers; ++iLayer){
1451             if ((layerMask >> iLayer) & 1){
1452               if (fTrackingData->GetTrackTrackletPID(iStack, iTrk, iLayer) < minValue){
1453                 minValue = fTrackingData->GetTrackTrackletPID(iStack, iTrk, iLayer);
1454                 minLyr = iLayer;
1455               }
1456               if (fTrackingData->GetTrackTrackletPID(iStack, iTrk, iLayer) > maxValue){
1457                 maxValue = fTrackingData->GetTrackTrackletPID(iStack, iTrk, iLayer);
1458                 maxLyr = iLayer;
1459               }
1460             }
1461           } // loop over layers
1462
1463           // calculate trunc mean
1464           for (UShort_t iLayer = 0; iLayer < fkTRDLayers; ++iLayer){
1465             if (((layerMask >> iLayer) & 1) && (iLayer != minLyr) && (iLayer != maxLyr)){
1466               truncMeanPID += fTrackingData->GetTrackTrackletPID(iStack, iTrk, iLayer);
1467               lyrCount++;
1468             }
1469           } // loop over layers
1470           truncMeanPID = TMath::Nint((Double_t)(truncMeanPID)/(Double_t)(lyrCount));
1471
1472           fHistPIDvsTruncPID->Fill(fTrackingData->GetTrackPID(iStack, iTrk), truncMeanPID);
1473           if (fTrackingData->GetTrackAddInfo(iStack, iTrk) < 0)
1474             fHistElectronFalsePIDvsTruncPID->Fill(fTrackingData->GetTrackPID(iStack, iTrk), truncMeanPID);
1475
1476         }
1477
1478         LogInspect("#hlt-trd-trg - GTU flagged %s %s high-pt electron seen in S%02d-%d: pt: %+6.1f  pid: %d  [id: S%02d-%d-%d] [trunc-pid: %d]",
1479                    (fTrackingData->GetTrackAddInfo(iStack, iTrk) < 0) ? "unmatched" : "matched", ident,
1480                    iStack/5, iStack%5,
1481                    gpt, fTrackingData->GetTrackPID(iStack, iTrk),
1482                    iStack/5, iStack%5, iTrk,
1483                    truncMeanPID);
1484       }
1485
1486       // evaluate HLT-level trigger decision using match information
1487       if (
1488           (fTrackingData->GetTrackAddInfo(iStack, iTrk) >= 0) &&
1489           (TMath::Abs(gpt) >= minPt) &&
1490           (fTrackingData->GetTrackPID(iStack, iTrk) >= minPID)
1491           ) {
1492         hltElectronCandidate = kTRUE;
1493         highPtElectronSeen = kTRUE;
1494
1495         fHistElectronCandidateMatchedPt->Fill(gpt);
1496         fHistElectronCandidateMatchedPid->Fill(fTrackingData->GetTrackPID(iStack, iTrk));
1497
1498         if (fExtendedHistos){
1499           fHistElectronConfirmedPIDvsTruncPID->Fill(fTrackingData->GetTrackPID(iStack, iTrk), truncMeanPID);
1500         }
1501
1502         Double_t rpt = (fESDtracksPresent) ? fEsdEvent->GetTrack(fTrackingData->GetTrackAddInfo(iStack, iTrk))->GetSignedPt() :
1503           fHLTTracks->at(fTrackingData->GetTrackAddInfo(iStack, iTrk)).GetSignedPt();
1504         LogInspect("#hlt-trd-trg - HLT matched %s high-pt electron seen in S%02d-%d: gpt: %+6.1f  gpid: %d   rpt: %+6.1f  [id: S%02d-%d-%d]",
1505                    ident, iStack/5, iStack%5, fTrackingData->GetTrackPt(iStack, iTrk), fTrackingData->GetTrackPID(iStack, iTrk),
1506                    rpt, iStack/5, iStack%5, iTrk);
1507       }
1508
1509       // log output for subsequent offline analysis
1510       if ((fDebugLevel >= 3) && (trdElectronCandidate || hltElectronCandidate)){
1511         trackExchangeInfo = "";
1512         flags = "";
1513         if (trdElectronCandidate)
1514           flags += "G";
1515         if (hltElectronCandidate)
1516           flags += "H";
1517         AssignTrackInfo(&trackExchangeInfo, iStack, iTrk, flags.Data());
1518         LogDebug("%s\n", trackExchangeInfo.Data());
1519       }
1520
1521     } // loop over tracks in stack
1522   } // loop over stacks
1523
1524   if (highPtElectronSeenGTU || highPtElectronSeen){
1525     LogInspect("#hlt-trd-trg - event triggered by %s electron trigger (TRD L1: %d, TRD HLT: %d)",
1526                ident, highPtElectronSeenGTU, highPtElectronSeen);
1527   }
1528
1529   if (highPtElectronSeenGTU){
1530     if (fIsMinBiasEvent)
1531       fHistElectronTriggerBaseMinBias->Fill(1.);
1532   }
1533
1534   if (highPtElectronSeen){
1535     if (fIsMinBiasEvent)
1536       fHistElectronTriggerBaseMinBias->Fill(2.);
1537     if (fIsTRDElectronEvent)
1538       fHistElectronTriggerBaseTrdL1->Fill(1.);
1539   }
1540
1541   return (highPtElectronSeen) ? kTRUE : kFALSE;
1542 }
1543
1544
1545 int AliHLTTRDTriggerComponent::DoTrigger()
1546 {
1547
1548   fEsdEvent = NULL;
1549   fRunNumber = -1;
1550   int iResult = 0;
1551   UShort_t firedTriggers = 0;
1552
1553   const AliHLTComponentEventData* hltEventData = GetEventData();
1554   fEventId = hltEventData->fEventID;
1555   LogDebug("### START DoTrigger [event id: %llu, %d blocks, size: %d]",
1556                       fEventId, hltEventData->fBlockCnt, hltEventData->fStructSize);
1557
1558 //  LogDebug("### START DoTrigger [event id: %llu, %d blocks, size: %d]",
1559 //                    fEventId, hltEventData->fBlockCnt, hltEventData->fStructSize);
1560
1561   if (!IsDataEvent()) {  // process data events only
1562     IgnoreEvent();
1563     LogDebug("### END   DoTrigger [event id: %llu, %d blocks, size: %d] (skipped: no data event)",
1564              fEventId, hltEventData->fBlockCnt, hltEventData->fStructSize);
1565     return iResult;
1566   }
1567
1568   fTrackingData->SetLogPrefix(Form("TRDHLTGM XXXXXX-%05llu: [TRG] {TrkDat} ", fEventId));
1569
1570   do {
1571
1572     // access to TRD specific data from AliHLTTRDPreprocessorComponent
1573     if (!PrepareTRDData()){
1574       LogError("access to TRD data failed. Skipping event...");
1575       break;
1576     }
1577
1578     if (fTrackingData->GetNumTracks() + fTrackingData->GetNumTracklets() == 0) {
1579       LogDebug("no trigger-relevant TRD information, skipping further event processing");
1580       break;
1581     }
1582
1583     // access to ESD data
1584     if (!PrepareESDData()){
1585       LogInfo("access to ESD event data failed.");
1586     }
1587
1588     // access to alternative HLT data
1589     if (!fESDtracksPresent){
1590       if (!PrepareHLTData()){
1591         LogError("access to HLT event data failed.");
1592       }
1593     }
1594
1595     // match TRD and HLT tracks
1596     if (fESDtracksPresent){
1597       if (!MatchTRDTracksESD()){
1598         LogError("matching TRD tracks to ESD tracks failed. Skipping event...");
1599         break;
1600       }
1601     } else if (fHLTtracksPresent){
1602       if (!MatchTRDTracksHLT()){
1603         LogError("matching TRD tracks to HLT tracks failed. Skipping event...");
1604         break;
1605       }
1606     } else {
1607       LogError("No HLT track information available. Skipping event...");
1608       break;
1609     }
1610
1611 //    if (!MatchTRDTracks()){
1612 //      LogError("matching TRD tracks to TPC tracks failed. Skipping event...");
1613 //      break;
1614 //    }
1615
1616     if (fDebugLevel >= 1)
1617       DumpTrackingData();
1618
1619     // evaluate electron trigger conditions
1620     if (TRDElectronTrigger("HSE", fElectronTriggerPtThresholdHSE, fElectronTriggerPIDThresholdHSE))
1621       firedTriggers |= fkElectronTriggerHSE;
1622
1623     if (TRDElectronTrigger("HQU", fElectronTriggerPtThresholdHQU, fElectronTriggerPIDThresholdHQU))
1624       firedTriggers |= fkElectronTriggerHQU;
1625
1626     break;
1627
1628   } while (1);
1629
1630
1631   // trigger decision
1632   TString description("");
1633   if (firedTriggers & fkElectronTriggerHSE){
1634     if (description.Length() > 0)
1635       description += " ";
1636     description += fgkTriggerDecisionElectronHSE;
1637   }
1638
1639   if (firedTriggers & fkElectronTriggerHQU){
1640     if (description.Length() > 0)
1641       description += " ";
1642     description += fgkTriggerDecisionElectronHQU;
1643   }
1644
1645   SetDescription(description.Data());
1646   AliHLTTriggerDecision decision((firedTriggers) ? kTRUE : kFALSE,
1647                                  GetTriggerName(),
1648                                  GetReadoutList(),
1649                                  GetDescription()
1650                                  );
1651   TriggerEvent(&decision, kAliHLTDataTypeTObject | kAliHLTDataOriginOut);
1652
1653   if (firedTriggers){
1654     LogInspect("TRD HLT trigger fired for event: description: >%s<, flags: 0x%04x",
1655                description.Data(), firedTriggers);
1656 #ifdef __TRDHLTDEBUG
1657     if (fEventRendering)
1658       RenderEvent();
1659 #endif
1660   } else {
1661     LogDebug("TRD HLT trigger did not fire for event");
1662   }
1663
1664   if (fPushHistos){
1665     PushBack(fHistArray, (kAliHLTDataTypeTObjArray | kAliHLTDataOriginTRD), 0x3fffff);
1666   }
1667
1668   LogDebug("### END   DoTrigger [event id: %llu, %d blocks, size: %d]",
1669                       fEventId, hltEventData->fBlockCnt, hltEventData->fStructSize);
1670
1671   return iResult;
1672 }
1673
1674 Bool_t AliHLTTRDTriggerComponent::TrackPlaneIntersect(AliExternalTrackParam *trk, Double_t pnt[3], Double_t norm[3], Double_t mag){
1675
1676   UInt_t its = 0;
1677   Double_t r = 290.;
1678   Double_t dist = 99999, dist_prev = 99999;
1679   Double_t x[3] = {0., 0., 0.};
1680   Bool_t ret = kTRUE;
1681
1682   dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) *norm[1] + (x[2] - pnt[2]) * norm[2];
1683
1684   while(TMath::Abs(dist) > 0.1) {
1685
1686     trk->GetXYZAt(r, mag, x);
1687
1688     if ((x[0] * x[0] + x[1] * x[1]) < 100.) {  // extrapolation to radius failed
1689       ret = kFALSE;
1690       break;
1691     }
1692
1693     //distance between current track position and plane
1694     dist_prev = TMath::Abs(dist);
1695     dist = (x[0] - pnt[0]) * norm[0] + (x[1] - pnt[1]) * norm[1];
1696     r -= dist;
1697     its++;
1698     if(TMath::Abs(dist) >= dist_prev ||
1699        (r > 380.) || (r < 100.)){
1700       ret = kFALSE;
1701       break;
1702     }
1703   }
1704
1705   for (Int_t i=0; i<3; i++){
1706     if(ret)
1707       pnt[i] = x[i];
1708     else
1709       pnt[i] = 0.;
1710   }
1711
1712   return kTRUE;
1713 }
1714
1715
1716 Double_t AliHLTTRDTriggerComponent::RateTrackMatch(Double_t distY, Double_t distZ, Double_t rpt, Double_t gpt){
1717
1718   // maximum limits for spatial distance
1719   if ((distY > 5.) || (distZ > 20.))
1720     return 0.;
1721
1722   // same pt sign required
1723   if ((rpt * gpt) < 0.)
1724     return 0.;
1725
1726   Double_t rating_distY = -0.1 * distY + 1.;
1727   Double_t rating_distZ = -0.025 * distZ + 1.;
1728   Double_t rating_ptDiff = 1. - TMath::Abs((TMath::Abs(rpt) > 0.000001) ? ((gpt-rpt)/rpt) : 0.);
1729
1730   if (rating_ptDiff <  0.)
1731     rating_ptDiff = 0.2;
1732
1733   Double_t total = rating_distY * rating_distZ * rating_ptDiff;
1734
1735 //  DbgLog("", Form("#matching: rating:   dy: %.3f   dz: %.3f   dpt: %.3f     -> total: %.3f",
1736 //                  rating_distY, rating_distZ, rating_ptDiff, total));
1737
1738   if (total > 1.)
1739     LogError("track match rating exceeds limit of 1.0: %.3f", total);
1740
1741   return total;
1742 }
1743
1744 Int_t AliHLTTRDTriggerComponent::EstimateTrackDistance(AliESDtrack *esd_track,
1745                                                         const UShort_t stack,
1746                                                         const UShort_t layerMask,
1747                                                         const Float_t trklLocalY[6], const Int_t trklBinZ[6],
1748                                                         Double_t mag, Double_t *ydist, Double_t *zdist){
1749   if (!esd_track)
1750     return -3;
1751
1752   AliExternalTrackParam* refParam = NULL;
1753   if (esd_track->GetOuterParam())
1754     refParam = new AliExternalTrackParam(*(esd_track->GetOuterParam()));
1755   if (!refParam)
1756     refParam = new AliExternalTrackParam(*(esd_track));
1757
1758   Int_t res = EstimateTrackDistance(refParam, stack, layerMask, trklLocalY, trklBinZ, mag, ydist, zdist);
1759
1760   if (refParam)
1761     delete refParam;
1762
1763   return res;
1764
1765 }
1766
1767 Int_t AliHLTTRDTriggerComponent::EstimateTrackDistance(AliExternalTrackParam *refParam,
1768                                                         const UShort_t stack,
1769                                                         const UShort_t layerMask,
1770                                                         const Float_t trklLocalY[6], const Int_t trklBinZ[6],
1771                                                         Double_t mag, Double_t *ydist, Double_t *zdist){
1772
1773   Float_t diff_y = 0;
1774   Float_t diff_z = 0;
1775   Int_t nLayers = 0;
1776   Double_t xtrkl[3];
1777   Double_t ptrkl[3];
1778   Double_t ptrkl2[3];
1779   UInt_t trklDet;
1780   UShort_t trklLayer;
1781   // UInt_t stack_gtu;
1782   UShort_t stackInSector;
1783
1784   AliTRDpadPlane* padPlane;
1785
1786   for (UShort_t iLayer = 0; iLayer < 6; iLayer++){
1787     if ((layerMask >> iLayer) & 1){
1788       trklDet = stack*6 + iLayer;
1789       trklLayer = iLayer;
1790       // stack_gtu = stack;
1791       stackInSector = stack % 5;
1792
1793       // local coordinates of the outer end point of the tracklet
1794       xtrkl[0] = AliTRDgeometry::AnodePos();
1795       xtrkl[1] = trklLocalY[iLayer];
1796
1797       padPlane = fTRDGeometry->GetPadPlane(trklLayer, stackInSector);
1798       if(stackInSector == 2){ // corrected version by Felix Muecke
1799         xtrkl[2] = padPlane->GetRowPos(trklBinZ[iLayer]) - (padPlane->GetRowSize(trklBinZ[iLayer]))/2. - padPlane->GetRowPos(6);
1800       } else {
1801         xtrkl[2] = padPlane->GetRowPos(trklBinZ[iLayer]) - (padPlane->GetRowSize(trklBinZ[iLayer]))/2. - padPlane->GetRowPos(8);
1802       }
1803
1804       // transform to global coordinates
1805       TGeoHMatrix *matrix = fTRDGeometry->GetClusterMatrix(trklDet);
1806       if (!matrix){
1807         LogError("invalid TRD cluster matrix in EstimateTrackDistance for detector %i", trklDet);
1808         return -5;
1809       }
1810       matrix->LocalToMaster(xtrkl, ptrkl);
1811       fTRDGeometry->RotateBack((stack/5) * 30, ptrkl, ptrkl2);  // ptrkl2 now contains the global position of the outer end point of the tracklet
1812
1813       // calculate parameterization of plane representing the tracklets layer
1814       Double_t layer_zero_local[3] = {0., 0.,   0.};
1815       Double_t layer_zero_global[3], layer_zero_global2[3];
1816
1817       matrix->LocalToMaster(layer_zero_local, layer_zero_global);
1818       fTRDGeometry->RotateBack(trklDet, layer_zero_global, layer_zero_global2); // layer_zero_global2 points to chamber origin in global coords
1819
1820       Double_t layer_ref_local[3] = {AliTRDgeometry::AnodePos(), 0.,  0.};
1821       Double_t layer_ref_global[3], layer_ref_global2[3];
1822
1823       matrix->LocalToMaster(layer_ref_local, layer_ref_global);
1824       fTRDGeometry->RotateBack(trklDet, layer_ref_global, layer_ref_global2); // layer_ref_global2 points to center anode pos within plane in global coords
1825
1826       Double_t n0[3] = {layer_ref_global2[0]-layer_zero_global2[0],
1827                         layer_ref_global2[1]-layer_zero_global2[1],
1828                         layer_ref_global2[2]-layer_zero_global2[2]};
1829
1830       Double_t n_len = TMath::Sqrt(n0[0]*n0[0] + n0[1]*n0[1] + n0[2]*n0[2]);
1831       if (n_len == 0.){ // This should never happen
1832         //printf("<ERROR> divison by zero in estimate_track_distance!");
1833         n_len = 1.;
1834       }
1835       Double_t n[3] = {n0[0]/n_len, n0[1]/n_len, n0[2]/n_len}; // normal vector of plane
1836
1837       Bool_t isects = TrackPlaneIntersect(refParam, layer_ref_global2, n, mag); // find intersection point between track and TRD layer
1838
1839       if (isects == kFALSE){ // extrapolation fails, because track never reaches the TRD radius
1840         return -1;
1841       }
1842
1843       Double_t m[2] = {ptrkl2[0] - layer_ref_global2[0], ptrkl2[1] - layer_ref_global2[1]};
1844       Double_t len_m = TMath::Sqrt(m[0]*m[0] + m[1]*m[1]);
1845       diff_y += len_m;
1846       diff_z += TMath::Abs(ptrkl2[2] - layer_ref_global2[2]);
1847       nLayers++;
1848     }
1849   }
1850
1851   if (nLayers > 0){
1852     *ydist = diff_y / nLayers;
1853     *zdist = diff_z / nLayers;
1854     return 0;
1855   } else {
1856     LogError("invalid number of contributing layers (%d) in EstimateTrackDistance()", nLayers);
1857     return -4;
1858   }
1859 }