HLT/MUON/OnlineAnalysis/AliHLTMUONDecisionComponent.cxx

   1 /**************************************************************************
   2  * This file is property of and copyright by the ALICE HLT Project        *
   3  * All rights reserved.                                                   *
   4  *                                                                        *
   5  * Primary Authors:                                                       *
   6  *   Artur Szostak <artursz@iafrica.com>                                  *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id: $ */
  18
  19 ///
  20 ///  @file   AliHLTMUONDecisionComponent.cxx
  21 ///  @author Artur Szostak <artursz@iafrica.com>
  22 ///  @date   30 April 2008
  23 ///  @brief  Implementation of the decision component for dimuon HLT triggering.
  24 ///
  25 // class documentation is in the header file.
  26
  27 #include "AliHLTMUONDecisionComponent.h"
  28 #include "AliHLTMUONConstants.h"
  29 #include "AliHLTMUONUtils.h"
  30 #include "AliHLTMUONCalculations.h"
  31 #include "AliHLTMUONDataBlockReader.h"
  32 #include "AliHLTMUONDataBlockWriter.h"
  33 #include "AliCDBEntry.h"
  34 #include "AliCDBManager.h"
  35 #include "TObjString.h"
  36 #include "TString.h"
  37 #include <cstdlib>
  38 #include <cstring>
  39 #include <cerrno>
  40 #include <cmath>
  41 #include <new>
  42
  43
  44 // Helper type for memory allocation.
  45 typedef const AliHLTMUONMansoTrackStruct* AliHLTMUONMansoTrackStructP;
  46
  47
  48 ClassImp(AliHLTMUONDecisionComponent);
  49
  50
  51 AliHLTMUONDecisionComponent::AliHLTMUONDecisionComponent() :
  52         AliHLTMUONProcessor(),
  53         fMaxTracks(1),
  54         fTrackCount(0),
  55         fTracks(new AliHLTMUONMansoTrackStructP[fMaxTracks]),
  56         fLowPtCut(1.),  // 1 GeV/c cut
  57         fHighPtCut(2.),  // 2 GeV/c cut
  58         fLowMassCut(2.5),  // 2.7 GeV/c^2 cut
  59         fHighMassCut(7.),  // 8 GeV/c^2 cut
  60         fWarnForUnexpecedBlock(false)
  61 {
  62         ///
  63         /// Default constructor.
  64         ///
  65 }
  66
  67
  68 AliHLTMUONDecisionComponent::~AliHLTMUONDecisionComponent()
  69 {
  70         ///
  71         /// Default destructor deletes the fTracks array.
  72         ///
  73
  74         assert(fTracks != NULL);
  75         delete [] fTracks;
  76 }
  77
  78
  79 const char* AliHLTMUONDecisionComponent::GetComponentID()
  80 {
  81         ///
  82         /// Inherited from AliHLTComponent. Returns the component ID.
  83         ///
  84
  85         return AliHLTMUONConstants::DecisionComponentId();
  86 }
  87
  88
  89 void AliHLTMUONDecisionComponent::GetInputDataTypes(
  90                 vector<AliHLTComponentDataType>& list
  91         )
  92 {
  93         ///
  94         /// Inherited from AliHLTProcessor. Returns the list of expected input data types.
  95         ///
  96
  97         assert( list.empty() );
  98         list.push_back( AliHLTMUONConstants::MansoTracksBlockDataType() );
  99 }
 100
 101
 102 AliHLTComponentDataType AliHLTMUONDecisionComponent::GetOutputDataType()
 103 {
 104         /// Inherited from AliHLTComponent. Returns kAliHLTMultipleDataType
 105         /// refer to GetOutputDataTypes for all returned data types.
 106
 107         return kAliHLTMultipleDataType;
 108 }
 109
 110
 111 int AliHLTMUONDecisionComponent::GetOutputDataTypes(AliHLTComponentDataTypeList& list)
 112 {
 113         /// Inherited from AliHLTComponent. Returns the output data types.
 114
 115         assert( list.empty() );
 116         list.push_back( AliHLTMUONConstants::SinglesDecisionBlockDataType() );
 117         list.push_back( AliHLTMUONConstants::PairsDecisionBlockDataType() );
 118         return 1;
 119 }
 120
 121
 122 void AliHLTMUONDecisionComponent::GetOutputDataSize(
 123                 unsigned long& constBase, double& inputMultiplier
 124         )
 125 {
 126         ///
 127         /// Inherited from AliHLTComponent. Returns an estimate of the expected output data size.
 128         ///
 129
 130         constBase = sizeof(AliHLTMUONSinglesDecisionBlockStruct);
 131         constBase += sizeof(AliHLTMUONPairsDecisionBlockStruct);
 132         inputMultiplier = 100;
 133 }
 134
 135
 136 AliHLTComponent* AliHLTMUONDecisionComponent::Spawn()
 137 {
 138         ///
 139         /// Inherited from AliHLTComponent. Creates a new object instance.
 140         ///
 141
 142         return new AliHLTMUONDecisionComponent;
 143 }
 144
 145
 146 int AliHLTMUONDecisionComponent::DoInit(int argc, const char** argv)
 147 {
 148         ///
 149         /// Inherited from AliHLTComponent.
 150         /// Parses the command line parameters and initialises the component.
 151         ///
 152
 153         HLTInfo("Initialising dHLT trigger decision component.");
 154
 155         bool lowPtCutSet = false;
 156         bool highPtCutSet = false;
 157         bool lowMassCutSet = false;
 158         bool highMassCutSet = false;
 159         fWarnForUnexpecedBlock = false;
 160
 161         for (int i = 0; i < argc; i++)
 162         {
 163                 if (strcmp( argv[i], "-lowptcut" ) == 0)
 164                 {
 165                         if (lowPtCutSet)
 166                         {
 167                                 HLTWarning("Low pT cut parameter was already specified."
 168                                         " Will replace previous value given by -lowptcut."
 169                                 );
 170                         }
 171
 172                         if (argc <= i+1)
 173                         {
 174                                 HLTError("The value for the low pT cut was not specified.");
 175                                 return -EINVAL;
 176                         }
 177
 178                         char* cpErr = NULL;
 179                         double num = strtod(argv[i+1], &cpErr);
 180                         if (cpErr == NULL or *cpErr != '\0')
 181                         {
 182                                 HLTError("Cannot convert '%s' to a floating point number.", argv[i+1]);
 183                                 return -EINVAL;
 184                         }
 185                         fLowPtCut = (AliHLTFloat32_t)num;
 186                         lowPtCutSet = true;
 187
 188                         i++;
 189                         continue;
 190                 }
 191
 192                 if (strcmp( argv[i], "-highptcut" ) == 0)
 193                 {
 194                         if (lowPtCutSet)
 195                         {
 196                                 HLTWarning("High pT cut parameter was already specified."
 197                                         " Will replace previous value given by -highptcut."
 198                                 );
 199                         }
 200
 201                         if (argc <= i+1)
 202                         {
 203                                 HLTError("The value for the high pT cut was not specified.");
 204                                 return -EINVAL;
 205                         }
 206
 207                         char* cpErr = NULL;
 208                         double num = strtod(argv[i+1], &cpErr);
 209                         if (cpErr == NULL or *cpErr != '\0')
 210                         {
 211                                 HLTError("Cannot convert '%s' to a floating point number.", argv[i+1]);
 212                                 return -EINVAL;
 213                         }
 214                         fHighPtCut = (AliHLTFloat32_t)num;
 215                         highPtCutSet = true;
 216
 217                         i++;
 218                         continue;
 219                 }
 220
 221                 if (strcmp( argv[i], "-lowmasscut" ) == 0)
 222                 {
 223                         if (lowPtCutSet)
 224                         {
 225                                 HLTWarning("Low invariant mass cut parameter was already specified."
 226                                         " Will replace previous value given by -lowmasscut."
 227                                 );
 228                         }
 229
 230                         if (argc <= i+1)
 231                         {
 232                                 HLTError("The value for the low invariant mass cut was not specified.");
 233                                 return -EINVAL;
 234                         }
 235
 236                         char* cpErr = NULL;
 237                         double num = strtod(argv[i+1], &cpErr);
 238                         if (cpErr == NULL or *cpErr != '\0')
 239                         {
 240                                 HLTError("Cannot convert '%s' to a floating point number.", argv[i+1]);
 241                                 return -EINVAL;
 242                         }
 243                         fLowMassCut = (AliHLTFloat32_t)num;
 244                         lowMassCutSet = true;
 245
 246                         i++;
 247                         continue;
 248                 }
 249
 250                 if (strcmp( argv[i], "-highmasscut" ) == 0)
 251                 {
 252                         if (lowPtCutSet)
 253                         {
 254                                 HLTWarning("High invariant mass cut parameter was already specified."
 255                                         " Will replace previous value given by -highmasscut."
 256                                 );
 257                         }
 258
 259                         if (argc <= i+1)
 260                         {
 261                                 HLTError("The value for the high invariant mass cut was not specified.");
 262                                 return -EINVAL;
 263                         }
 264
 265                         char* cpErr = NULL;
 266                         double num = strtod(argv[i+1], &cpErr);
 267                         if (cpErr == NULL or *cpErr != '\0')
 268                         {
 269                                 HLTError("Cannot convert '%s' to a floating point number.", argv[i+1]);
 270                                 return -EINVAL;
 271                         }
 272                         fHighMassCut = (AliHLTFloat32_t)num;
 273                         highMassCutSet = true;
 274
 275                         i++;
 276                         continue;
 277                 }
 278
 279                 if (strcmp(argv[i], "-warn_on_unexpected_block") == 0)
 280                 {
 281                         fWarnForUnexpecedBlock = true;
 282                         continue;
 283                 }
 284
 285                 HLTError("Unknown option '%s'.", argv[i]);
 286                 return -EINVAL;
 287         }
 288
 289         // Read cut parameters from CDB if they were not specified on the command line.
 290         if (not lowPtCutSet or not highPtCutSet or not lowMassCutSet or not highMassCutSet)
 291         {
 292                 int result = ReadConfigFromCDB(
 293                                 NULL,
 294                                 not lowPtCutSet, not highPtCutSet,
 295                                 not lowMassCutSet, not highMassCutSet
 296                         );
 297                 if (result != 0) return result;
 298         }
 299
 300         HLTDebug("Using the following cut parameters:");
 301         HLTDebug("              Low pT cut = %f GeV/c", fLowPtCut);
 302         HLTDebug("             High pT cut = %f GeV/c", fHighPtCut);
 303         HLTDebug("  Low invariant mass cut = %f GeV/c^2", fLowMassCut);
 304         HLTDebug(" High invariant mass cut = %f GeV/c^2", fHighMassCut);
 305
 306         return 0;
 307 }
 308
 309
 310 int AliHLTMUONDecisionComponent::DoDeinit()
 311 {
 312         ///
 313         /// Inherited from AliHLTComponent. Performs a cleanup of the component.
 314         ///
 315
 316         HLTInfo("Deinitialising dHLT trigger decision component.");
 317         return 0;
 318 }
 319
 320
 321 int AliHLTMUONDecisionComponent::Reconfigure(const char* cdbEntry, const char* componentId)
 322 {
 323         /// Inherited from AliHLTComponent. Reconfigures the component from CDB.
 324
 325         if (strcmp(componentId, GetComponentID()) == 0)
 326         {
 327                 HLTInfo("Reading new entries for cut parameters from CDB.");
 328                 int result = ReadConfigFromCDB(cdbEntry);
 329                 HLTDebug("Using the following new cut parameters:");
 330                 HLTDebug("              Low pT cut = %f GeV/c", fLowPtCut);
 331                 HLTDebug("             High pT cut = %f GeV/c", fHighPtCut);
 332                 HLTDebug("  Low invariant mass cut = %f GeV/c^2", fLowMassCut);
 333                 HLTDebug(" High invariant mass cut = %f GeV/c^2", fHighMassCut);
 334                 return result;
 335         }
 336         else
 337                 return 0;
 338 }
 339
 340
 341 int AliHLTMUONDecisionComponent::DoEvent(
 342                 const AliHLTComponentEventData& evtData,
 343                 const AliHLTComponentBlockData* blocks,
 344                 AliHLTComponentTriggerData& /*trigData*/,
 345                 AliHLTUInt8_t* outputPtr,
 346                 AliHLTUInt32_t& size,
 347                 std::vector<AliHLTComponentBlockData>& outputBlocks
 348         )
 349 {
 350         ///
 351         /// Inherited from AliHLTProcessor. Processes the new event data.
 352         ///
 353
 354         AliHLTUInt32_t specification = 0;  // Contains the output data block spec bits.
 355
 356         // Loop over all input blocks in the event with track data and add pointers
 357         // to the tracks into the tracks array. These will be used later by the
 358         // trigger algorithm to get to the individual tracks.
 359         fTrackCount = 0; // reset number of tracks in array.
 360         for (AliHLTUInt32_t n = 0; n < evtData.fBlockCnt; n++)
 361         {
 362                 HLTDebug("Handling block: %u, with fDataType = '%s', fPtr = %p and fSize = %u bytes.",
 363                         n, DataType2Text(blocks[n].fDataType).c_str(), blocks[n].fPtr, blocks[n].fSize
 364                 );
 365
 366                 if (blocks[n].fDataType == AliHLTMUONConstants::MansoTracksBlockDataType())
 367                 {
 368                         // Build up the specification which indicates what DDLs
 369                         // contributed to the output data.
 370                         specification |= blocks[n].fSpecification;
 371
 372                         AliHLTMUONMansoTracksBlockReader inblock(
 373                                         reinterpret_cast<char*>(blocks[n].fPtr) + blocks[n].fOffset,
 374                                         blocks[n].fSize
 375                                 );
 376                         if (not BlockStructureOk(inblock)) continue;
 377
 378                         for (AliHLTUInt32_t i = 0; i < inblock.Nentries(); i++)
 379                         {
 380                                 int result = AddTrack(&inblock[i]);
 381                                 if (result != 0)
 382                                 {
 383                                         size = 0; // Important to tell framework that nothing was generated.
 384                                         return result;
 385                                 }
 386                         }
 387                 }
 388                 else
 389                 {
 390                         // Log a message indicating that we got a data block that we
 391                         // do not know how to handle.
 392                         if (fWarnForUnexpecedBlock)
 393                                 HLTWarning("Received a data block of a type we cannot handle: '%s', spec: 0x%X",
 394                                         DataType2Text(blocks[n].fDataType).c_str(), blocks[n].fSpecification
 395                                 );
 396                         else
 397                                 HLTDebug("Received a data block of a type we cannot handle: '%s', spec: 0x%X",
 398                                         DataType2Text(blocks[n].fDataType).c_str(), blocks[n].fSpecification
 399                                 );
 400                 }
 401         }
 402
 403         // Now we can create our two new output data blocks for the single tracks
 404         // and track pairs.
 405         AliHLTMUONSinglesDecisionBlockWriter singlesBlock(outputPtr, size);
 406
 407         if (not singlesBlock.InitCommonHeader())
 408         {
 409                 Logging(kHLTLogError,
 410                         "AliHLTMUONDecisionComponent::DoEvent",
 411                         "Buffer overflow",
 412                         "The buffer is only %d bytes in size. We need a minimum of"
 413                         " %d bytes for the singles output data block.",
 414                         size, sizeof(AliHLTMUONSinglesDecisionBlockWriter::HeaderType)
 415                 );
 416                 size = 0; // Important to tell framework that nothing was generated.
 417                 return -ENOBUFS;
 418         }
 419
 420         if (not singlesBlock.SetNumberOfEntries(fTrackCount))
 421         {
 422                 AliHLTUInt32_t bytesneeded = sizeof(AliHLTMUONSinglesDecisionBlockWriter::HeaderType)
 423                         + fTrackCount * sizeof(AliHLTMUONSinglesDecisionBlockWriter::ElementType);
 424                 HLTError("The buffer is only %d bytes in size. We need a minimum of"
 425                         " %d bytes for the singles output data block.",
 426                         size, bytesneeded
 427                 );
 428                 size = 0; // Important to tell framework that nothing was generated.
 429                 return -ENOBUFS;
 430         }
 431
 432         AliHLTMUONPairsDecisionBlockWriter pairsBlock(
 433                         outputPtr + singlesBlock.BytesUsed(),
 434                         size - singlesBlock.BytesUsed()
 435                 );
 436
 437         if (not pairsBlock.InitCommonHeader())
 438         {
 439                 Logging(kHLTLogError,
 440                         "AliHLTMUONDecisionComponent::DoEvent",
 441                         "Buffer overflow",
 442                         "The buffer is only %d bytes in size. We need a minimum of"
 443                         " %d bytes for the pairs output data block.",
 444                         size,
 445                         sizeof(AliHLTMUONPairsDecisionBlockWriter::HeaderType) + singlesBlock.BytesUsed()
 446                 );
 447                 size = 0; // Important to tell framework that nothing was generated.
 448                 return -ENOBUFS;
 449         }
 450
 451         AliHLTUInt32_t numOfPairs = fTrackCount * (fTrackCount-1) / 2;
 452         if (not pairsBlock.SetNumberOfEntries(numOfPairs))
 453         {
 454                 AliHLTUInt32_t bytesneeded = sizeof(AliHLTMUONPairsDecisionBlockWriter::HeaderType)
 455                         + numOfPairs * sizeof(AliHLTMUONPairsDecisionBlockWriter::ElementType)
 456                         + singlesBlock.BytesUsed();
 457                 HLTError("The buffer is only %d bytes in size. We need a minimum of"
 458                         " %d bytes for the pairs output data block.",
 459                         size, bytesneeded
 460                 );
 461                 size = 0; // Important to tell framework that nothing was generated.
 462                 return -ENOBUFS;
 463         }
 464
 465         ApplyTriggerAlgorithm(
 466                         singlesBlock.BlockHeader(),
 467                         singlesBlock.GetArray(),
 468                         pairsBlock.BlockHeader(),
 469                         pairsBlock.GetArray()
 470                 );
 471
 472         AliHLTComponentBlockData sbd;
 473         FillBlockData(sbd);
 474         sbd.fPtr = outputPtr;
 475         sbd.fOffset = 0;
 476         sbd.fSize = singlesBlock.BytesUsed();
 477         sbd.fDataType = AliHLTMUONConstants::SinglesDecisionBlockDataType();
 478         sbd.fSpecification = specification;
 479         outputBlocks.push_back(sbd);
 480         size = singlesBlock.BytesUsed();
 481
 482         AliHLTComponentBlockData pbd;
 483         FillBlockData(pbd);
 484         pbd.fPtr = outputPtr;
 485         pbd.fOffset = singlesBlock.BytesUsed();
 486         pbd.fSize = pairsBlock.BytesUsed();
 487         pbd.fDataType = AliHLTMUONConstants::PairsDecisionBlockDataType();
 488         pbd.fSpecification = specification;
 489         outputBlocks.push_back(pbd);
 490         size += pairsBlock.BytesUsed();
 491
 492         return 0;
 493 }
 494
 495
 496 int AliHLTMUONDecisionComponent::ReadConfigFromCDB(
 497                 const char* path,
 498                 bool setLowPtCut, bool setHighPtCut,
 499                 bool setLowMassCut, bool setHighMassCut
 500         )
 501 {
 502         /// Reads the cut parameters from the CDB.
 503         /// \param path Indicates the partial (or relative) path to the CDB entry.
 504
 505         assert(AliCDBManager::Instance() != NULL);
 506
 507         const char* pathToEntry = AliHLTMUONConstants::DecisionComponentCDBPath();
 508         if (path != NULL)
 509                 pathToEntry = path;
 510
 511         TMap* map = NULL;
 512         int result = FetchTMapFromCDB(pathToEntry, map);
 513         if (result != 0) return result;
 514
 515         if (setLowPtCut)
 516         {
 517                 Double_t value = 0;
 518                 result = GetFloatFromTMap(map, "lowptcut", value, pathToEntry, "low pT cut");
 519                 if (result != 0) return result;
 520                 fLowPtCut = (AliHLTFloat32_t) value;
 521         }
 522
 523         if (setHighPtCut)
 524         {
 525                 Double_t value = 0;
 526                 result = GetFloatFromTMap(map, "highptcut", value, pathToEntry, "high pT cut");
 527                 if (result != 0) return result;
 528                 fHighPtCut = (AliHLTFloat32_t) value;
 529         }
 530
 531         if (setLowMassCut)
 532         {
 533                 Double_t value = 0;
 534                 result = GetFloatFromTMap(map, "lowmasscut", value, pathToEntry, "low invariant mass cut");
 535                 if (result != 0) return result;
 536                 fLowMassCut = (AliHLTFloat32_t) value;
 537         }
 538
 539         if (setHighMassCut)
 540         {
 541                 Double_t value = 0;
 542                 result = GetFloatFromTMap(map, "highmasscut", value, pathToEntry, "high invariant mass cut");
 543                 if (result != 0) return result;
 544                 fHighMassCut = (AliHLTFloat32_t) value;
 545         }
 546
 547         return 0;
 548 }
 549
 550
 551 int AliHLTMUONDecisionComponent::AddTrack(const AliHLTMUONMansoTrackStruct* track)
 552 {
 553         /// Adds a track to the internal track list for future reference in
 554         /// ApplyTriggerAlgorithm when we actually apply the trigger algorithm.
 555
 556         assert(fTrackCount <= fMaxTracks);
 557         assert(fTracks != NULL);
 558
 559         if (fTrackCount == fMaxTracks)
 560         {
 561                 // Buffer full so we need to resize it.
 562                 const AliHLTMUONMansoTrackStruct** tmp = NULL;
 563                 try
 564                 {
 565                         tmp = new AliHLTMUONMansoTrackStructP[fMaxTracks+1];
 566                 }
 567                 catch (const std::bad_alloc&)
 568                 {
 569                         HLTError("Could not allocate more memory for the track array.");
 570                         return -ENOMEM;
 571                 }
 572
 573                 // Copy over the exisiting data and then delete the old array.
 574                 memcpy(tmp, fTracks, sizeof(AliHLTMUONMansoTrackStructP)*fTrackCount);
 575                 delete [] fTracks;
 576                 fTracks = tmp;
 577                 fMaxTracks = fMaxTracks+1;
 578         }
 579
 580         fTracks[fTrackCount] = track;
 581         fTrackCount++;
 582         return 0;
 583 }
 584
 585
 586 void AliHLTMUONDecisionComponent::ApplyTriggerAlgorithm(
 587                 AliHLTMUONSinglesDecisionBlockStruct& singlesHeader,
 588                 AliHLTMUONTrackDecisionStruct* singlesDecision,
 589                 AliHLTMUONPairsDecisionBlockStruct& pairsHeader,
 590                 AliHLTMUONPairDecisionStruct* pairsDecision
 591         )
 592 {
 593         /// This method applies the dHLT trigger decision algorithm to all the
 594         /// tracks found in the input data.
 595
 596         // Zero the trigger counters for single tracks.
 597         singlesHeader.fNlowPt = 0;
 598         singlesHeader.fNhighPt = 0;
 599
 600         // Zero the trigger counters for pairs.
 601         pairsHeader.fNunlikeAnyPt = 0;
 602         pairsHeader.fNunlikeLowPt = 0;
 603         pairsHeader.fNunlikeHighPt = 0;
 604         pairsHeader.fNlikeAnyPt = 0;
 605         pairsHeader.fNlikeLowPt = 0;
 606         pairsHeader.fNlikeHighPt = 0;
 607         pairsHeader.fNmassAny = 0;
 608         pairsHeader.fNmassLow = 0;
 609         pairsHeader.fNmassHigh = 0;
 610
 611         // For the single tracks we check if a track has pT larger than either
 612         // the low or high pT cut. If it does then we increment the appropriate
 613         // counters in the header.
 614         for (AliHLTUInt32_t n = 0; n < fTrackCount; n++)
 615         {
 616                 const AliHLTMUONMansoTrackStruct* track = fTracks[n];
 617                 AliHLTMUONTrackDecisionStruct& decision = singlesDecision[n];
 618
 619                 bool passedHighPtCut = false;
 620                 bool passedLowPtCut = false;
 621
 622                 AliHLTFloat32_t pt = sqrt(track->fPx * track->fPx + track->fPy * track->fPy);
 623
 624                 if (pt > fHighPtCut)
 625                 {
 626                         passedHighPtCut = true;
 627                         singlesHeader.fNhighPt++;
 628                 }
 629                 if (pt > fLowPtCut)
 630                 {
 631                         passedLowPtCut = true;
 632                         singlesHeader.fNlowPt++;
 633                 }
 634
 635                 decision.fTrackId = track->fId;
 636                 decision.fTriggerBits = AliHLTMUONUtils::PackTrackDecisionBits(
 637                                 passedHighPtCut, passedLowPtCut
 638                         );
 639                 decision.fPt = pt;
 640         }
 641
 642         // Now we generate all the possible pairs of tracks and fill in the
 643         // trigger information. This will consist of calculating the invariant
 644         // mass for the pair, checking if it passes the low or high mass cut
 645         // and incrementing the appropriate statistics.
 646         AliHLTUInt32_t currentPair = 0;
 647         for (AliHLTUInt32_t i = 0; i < fTrackCount; i++)
 648         for (AliHLTUInt32_t j = i+1; j < fTrackCount; j++)
 649         {
 650                 const AliHLTMUONMansoTrackStruct* tracki = fTracks[i];
 651                 const AliHLTMUONMansoTrackStruct* trackj = fTracks[j];
 652                 const AliHLTMUONTrackDecisionStruct& trackidecision = singlesDecision[i];
 653                 const AliHLTMUONTrackDecisionStruct& trackjdecision = singlesDecision[j];
 654                 AliHLTMUONPairDecisionStruct& decision = pairsDecision[currentPair];
 655
 656                 AliHLTFloat32_t muMass = 0.1056583568; // muon mass in GeV/c^2
 657
 658                 AliHLTFloat32_t mass = AliHLTMUONCalculations::ComputeMass(
 659                                 muMass, tracki->fPx, tracki->fPy, tracki->fPz,
 660                                 muMass, trackj->fPx, trackj->fPy, trackj->fPz
 661                         );
 662
 663                 AliHLTMUONParticleSign signi, signj;
 664                 bool hitset[4];
 665                 AliHLTMUONUtils::UnpackMansoTrackFlags(tracki->fFlags, signi, hitset);
 666                 AliHLTMUONUtils::UnpackMansoTrackFlags(trackj->fFlags, signj, hitset);
 667
 668                 AliHLTUInt8_t highPtCount = 0;
 669                 if (trackidecision.fPt > fHighPtCut) highPtCount++;
 670                 if (trackjdecision.fPt > fHighPtCut) highPtCount++;
 671                 AliHLTUInt8_t lowPtCount = 0;
 672                 if (trackidecision.fPt > fLowPtCut) lowPtCount++;
 673                 if (trackjdecision.fPt > fLowPtCut) lowPtCount++;
 674
 675                 bool unlikeSign = (signi == kSignMinus and signj == kSignPlus) or
 676                                   (signi == kSignPlus  and signj == kSignMinus);
 677
 678                 bool passedHighMassCut = false;
 679                 bool passedLowMassCut = false;
 680                 if (unlikeSign)
 681                 {
 682                         pairsHeader.fNunlikeAnyPt++;
 683                         if (lowPtCount == 2) pairsHeader.fNunlikeLowPt++;
 684                         if (highPtCount == 2) pairsHeader.fNunlikeHighPt++;
 685
 686                         if (mass > fHighMassCut)
 687                         {
 688                                 passedHighMassCut = true;
 689                                 if (highPtCount == 2) pairsHeader.fNmassHigh++;
 690                         }
 691                         if (mass > fLowMassCut)
 692                         {
 693                                 passedLowMassCut = true;
 694                                 pairsHeader.fNmassAny++;
 695                                 if (lowPtCount == 2) pairsHeader.fNmassLow++;
 696                         }
 697                 }
 698                 else
 699                 {
 700                         pairsHeader.fNlikeAnyPt++;
 701                         if (lowPtCount == 2) pairsHeader.fNlikeLowPt++;
 702                         if (highPtCount == 2) pairsHeader.fNlikeHighPt++;
 703                 }
 704
 705                 decision.fTrackAId = tracki->fId;
 706                 decision.fTrackBId = trackj->fId;
 707                 decision.fTriggerBits = AliHLTMUONUtils::PackPairDecisionBits(
 708                                 passedHighMassCut, passedLowMassCut, unlikeSign,
 709                                 highPtCount, lowPtCount
 710                         );
 711                 decision.fInvMass = mass;
 712
 713                 currentPair++;
 714         }
 715
 716         assert( currentPair == fTrackCount * (fTrackCount-1) / 2 );
 717 }
 718