HLT/MUON/macros/MakeTriggerTable.C

   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 // Simple macro to generate N-tuple for performance analysis.
  18 // This macro must be compiled and run like so from within the aliroot command
  19 // prompt:
  20 //   root [0] gSystem->Load("libAliHLTMUON");
  21 //   root [0] .x MakeHitsTable.C+
  22
  23 #include "TVector3.h"
  24 #include "TSystem.h"
  25 #include "TROOT.h"
  26 #include "TParticle.h"
  27 #include "TArrayF.h"
  28 #include "TStopwatch.h"
  29 #include "TNtuple.h"
  30 #include "TFile.h"
  31 #include "TError.h"
  32
  33 #include "AliCDBManager.h"
  34 #include "AliLoader.h"
  35 #include "AliRunLoader.h"
  36 #include "AliStack.h"
  37 #include "AliRun.h"
  38 #include "AliMUON.h"
  39 #include "AliMUONMCDataInterface.h"
  40 #include "AliMUONHit.h"
  41 #include "AliHLTMUONEvent.h"
  42 #include "AliHLTMUONRootifierComponent.h"
  43 #include "AliHLTMUONMansoTrack.h"
  44
  45 #include <cstdlib>
  46 #include <vector>
  47 #include <iostream>
  48 using std::cout;
  49 using std::cerr;
  50 using std::endl;
  51
  52
  53 void MakeTriggerTable(
  54                 Int_t firstEvent = 0,
  55                 Int_t lastEvent = -1,
  56                 const char* L0outputfile = "output.root",
  57                 Float_t maxSigma = 4., // 4 standard deviations
  58                 Float_t sigmaXtrg = 0.5,  // 5 mm resolution
  59                 Float_t sigmaYtrg = 0.5,  // 5 mm resolution
  60                 Float_t sigmaZtrg = 0.02  // 2 microns resolution
  61         )
  62 {
  63         // Setup the CDB default storage and run number if nothing was set.
  64         AliCDBManager* cdbManager = AliCDBManager::Instance();
  65         if (cdbManager == NULL)
  66         {
  67                 cerr << "ERROR: Global CDB manager object does not exist." << endl;
  68                 return;
  69         }
  70         if (cdbManager->GetDefaultStorage() == NULL)
  71         {
  72                 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  73         }
  74         if (cdbManager->GetRun() == -1)
  75         {
  76                 cdbManager->SetRun(0);
  77         }
  78
  79         gSystem->Load("libAliHLTMUON");
  80         // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
  81
  82         TString fieldnames = "event:isprimary:pdgcode:sign:px:py:pz";
  83         for (Int_t i = 11; i <= 14; i++)
  84         {
  85                 fieldnames += ":x";
  86                 fieldnames += i;
  87                 fieldnames += ":y";
  88                 fieldnames += i;
  89                 fieldnames += ":z";
  90                 fieldnames += i;
  91         }
  92         fieldnames += ":L0sign:L0px:L0py:L0pz:L0fit";
  93         for (Int_t i = 11; i <= 14; i++)
  94         {
  95                 fieldnames += ":L0x";
  96                 fieldnames += i;
  97                 fieldnames += ":L0y";
  98                 fieldnames += i;
  99                 fieldnames += ":L0z";
 100                 fieldnames += i;
 101         }
 102         TNtuple* table = new TNtuple("triggertable", "Table of triggers.", fieldnames);
 103
 104         TFile L0file(L0outputfile, "READ");
 105
 106         AliMUONMCDataInterface data;
 107         AliStack* stack;
 108
 109         Int_t nEvents = data.NumberOfEvents();
 110         if (lastEvent < 0) lastEvent = nEvents-1;
 111         for (Int_t event = firstEvent; event <= lastEvent; event++)
 112         {
 113                 cout << "Processing event: " << event;
 114
 115                 char buf[1024];
 116                 char* str = &buf[0];
 117                 sprintf(str, "AliHLTMUONEvent;%d", event+1);
 118                 AliHLTMUONEvent* L0event = dynamic_cast<AliHLTMUONEvent*>( L0file.Get(str) );
 119                 if (L0event == NULL)
 120                 {
 121                         cout << endl;
 122                         cerr << "ERROR: Could not find " << str << " in " << L0outputfile << endl;
 123                         continue;
 124                 }
 125
 126                 data.GetEvent(event);
 127
 128                 stack = data.Stack(event);
 129                 if (stack == NULL)
 130                 {
 131                         cout << endl;
 132                         cerr << "ERROR: Stack is NULL" << endl;
 133                         continue;
 134                 }
 135
 136                 Int_t trackcount = data.NumberOfParticles();
 137                 cout << "\ttrackcount = " << trackcount << endl;
 138
 139                 TArrayF isPrimary(trackcount);
 140                 TArrayF pdgcode(trackcount);
 141                 TArrayF sign(trackcount);
 142                 TArrayF px(trackcount);
 143                 TArrayF py(trackcount);
 144                 TArrayF pz(trackcount);
 145                 TArrayF L0sign(trackcount);
 146                 TArrayF L0px(trackcount);
 147                 TArrayF L0py(trackcount);
 148                 TArrayF L0pz(trackcount);
 149                 TArrayF L0fit(trackcount);
 150
 151                 //cout << "Fill particle data" << endl;
 152                 for (Int_t i = 0; i < trackcount; i++)
 153                 {
 154                         TParticle* p = data.Particle(i);
 155                         isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
 156                         pdgcode[i] = p->GetPdgCode();
 157                         if (p->GetPDG() != NULL)
 158                         {
 159                                 Double_t charge = p->GetPDG()->Charge();
 160                                 if (charge > 0)
 161                                         sign[i] = 1;
 162                                 else if (charge < 0)
 163                                         sign[i] = -1;
 164                                 else
 165                                         sign[i] = 0;
 166                         }
 167                         else
 168                         {
 169                                 sign[i] = 0;
 170                                 cerr << "ERROR: Could not get the PDG information." << endl;
 171                         }
 172                         px[i] = p->Px();
 173                         py[i] = p->Py();
 174                         pz[i] = p->Pz();
 175                         L0sign[i] = 0;
 176                         L0px[i] = 0;
 177                         L0py[i] = 0;
 178                         L0pz[i] = 0;
 179                         L0fit[i] = -1;
 180                 }
 181
 182                 TArrayF hitX[4];
 183                 TArrayF hitY[4];
 184                 TArrayF hitZ[4];
 185                 TArrayF L0hitX[4];
 186                 TArrayF L0hitY[4];
 187                 TArrayF L0hitZ[4];
 188                 for (Int_t i = 0; i < 4; i++)
 189                 {
 190                         hitX[i].Set(trackcount);
 191                         hitY[i].Set(trackcount);
 192                         hitZ[i].Set(trackcount);
 193                         L0hitX[i].Set(trackcount);
 194                         L0hitY[i].Set(trackcount);
 195                         L0hitZ[i].Set(trackcount);
 196                         for (Int_t j = 0; j < trackcount; j++)
 197                         {
 198                                 hitX[i][j] = 0;
 199                                 hitY[i][j] = 0;
 200                                 hitZ[i][j] = 0;
 201                                 L0hitX[i][j] = 0;
 202                                 L0hitY[i][j] = 0;
 203                                 L0hitZ[i][j] = 0;
 204                         }
 205                 }
 206
 207                 //cout << "Fill hits" << endl;
 208                 for (Int_t i = 0; i < data.NumberOfTracks(); i++)
 209                 for (Int_t j = 0; j < data.NumberOfHits(i); j++)
 210                 {
 211                         AliMUONHit* hit = data.Hit(i, j);
 212                         Int_t chamber = hit->Chamber() - 1;
 213                         if (chamber < 0 || chamber >= 14)
 214                         {
 215                                 cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
 216                                 continue;
 217                         }
 218
 219                         //cout << "hit->Track() = " << hit->Track() << endl;
 220                         if (hit->Track() < 0 || hit->Track() >= trackcount)
 221                         {
 222                                 cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
 223                                         << trackcount << ")" << endl;
 224                                 continue;
 225                         }
 226                         Int_t particleindex = hit->Track();
 227
 228                         if (chamber >= 10)
 229                         {
 230                                 hitX[chamber-10][particleindex] = hit->X();
 231                                 hitY[chamber-10][particleindex] = hit->Y();
 232                                 hitZ[chamber-10][particleindex] = hit->Z();
 233                         }
 234                 } // hits loop
 235
 236                 // Ok now go through the L0 information for this event and correlate
 237                 // the L0 trigger records to the kine data.
 238                 // This is done as follows: take a trigger record. Then for each kine track
 239                 // calculate the distance between each hit of the kine track and corresponding
 240                 // L0 trigger track fragment hit.
 241                 // Then sum up the distances and divide by the number of hits matched.
 242                 // This gives us a fit quality /Chi2 type parameter. The kine track that
 243                 // has the smallest sum of differences is then chosen as the kine track
 244                 // correlating to the current L0 trigger record being processed.
 245                 for (Int_t i = 0; i < L0event->Array().GetEntriesFast(); i++)
 246                 {
 247                         if (L0event->Array().At(i)->IsA() != AliHLTMUONTriggerRecord::Class())
 248                                 continue;
 249                         AliHLTMUONTriggerRecord* trigrec = static_cast<AliHLTMUONTriggerRecord*>(L0event->Array().At(i));
 250
 251                         // Find the best fitting kine track.
 252                         Double_t bestFitQuality = 1e30;
 253                         Int_t bestFitTrack = -1;
 254                         for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
 255                         {
 256                                 Double_t sumOfDiffs = 0;
 257                                 Double_t hitsMatched = 0;
 258                                 for (Int_t j = 11; j <= 14; j++)
 259                                 {
 260                                         const TVector3& hit = trigrec->Hit(j);
 261                                         if (hit == TVector3(0,0,0)) continue;
 262                                         TVector3 hV(hitX[j-11][kinetrack], hitY[j-11][kinetrack], hitZ[j-11][kinetrack]);
 263                                         if (hV == TVector3(0,0,0)) continue;
 264                                         TVector3 diff = hV - hit;
 265                                         Double_t diffX = diff.X() / sigmaXtrg;
 266                                         Double_t diffY = diff.Y() / sigmaYtrg;
 267                                         Double_t diffZ = diff.Z() / sigmaZtrg;
 268                                         if (diffX > maxSigma) continue;
 269                                         if (diffY > maxSigma) continue;
 270                                         if (diffZ > maxSigma) continue;
 271                                         sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
 272                                         hitsMatched++;
 273                                 }
 274
 275                                 // Now check the fit quality.
 276                                 if (hitsMatched <= 0) continue;
 277                                 Double_t fitQuality = sumOfDiffs / hitsMatched;
 278                                 if (fitQuality < bestFitQuality)
 279                                 {
 280                                         bestFitQuality = fitQuality;
 281                                         bestFitTrack = kinetrack;
 282                                 }
 283                         }
 284
 285                         if (bestFitTrack < 0)
 286                         {
 287                                 // No track was matched so we need to add a fake track record.
 288                                 Float_t args[7+4*3+5+4*3];
 289                                 args[0] = event;
 290                                 args[1] = -1;
 291                                 args[2] = 0;
 292                                 args[3] = 0;
 293                                 args[4] = 0;
 294                                 args[5] = 0;
 295                                 args[6] = 0;
 296                                 for (Int_t j = 0; j < 4; j++)
 297                                 {
 298                                         args[7+j*3+0] = 0;
 299                                         args[7+j*3+1] = 0;
 300                                         args[7+j*3+2] = 0;
 301                                 }
 302                                 args[7+4*3+0] = trigrec->Sign();
 303                                 args[7+4*3+1] = trigrec->Px();
 304                                 args[7+4*3+2] = trigrec->Py();
 305                                 args[7+4*3+3] = trigrec->Pz();
 306                                 args[7+4*3+4] = -1;
 307                                 for (Int_t j = 0; j < 4; j++)
 308                                 {
 309                                         const TVector3& hit = trigrec->Hit(j+11);
 310                                         args[7+4*3+5+j*3+0] = hit.X();
 311                                         args[7+4*3+5+j*3+1] = hit.Y();
 312                                         args[7+4*3+5+j*3+2] = hit.Z();
 313                                 }
 314                                 table->Fill(args);
 315                         }
 316                         else
 317                         {
 318                                 // Fill the details about the L0 track to the best fitting track info.
 319                                 L0sign[bestFitTrack] = trigrec->Sign();
 320                                 L0px[bestFitTrack] = trigrec->Px();
 321                                 L0py[bestFitTrack] = trigrec->Py();
 322                                 L0pz[bestFitTrack] = trigrec->Pz();
 323                                 L0fit[bestFitTrack] = bestFitQuality;
 324                                 for (Int_t j = 11; j <= 14; j++)
 325                                 {
 326                                         const TVector3& hit = trigrec->Hit(j);
 327                                         if (hit == TVector3(0,0,0)) continue;
 328                                         L0hitX[j-11][bestFitTrack] = hit.X();
 329                                         L0hitY[j-11][bestFitTrack] = hit.Y();
 330                                         L0hitZ[j-11][bestFitTrack] = hit.Z();
 331                                 }
 332                         }
 333                 }
 334
 335                 //cout << "Fill table" << endl;
 336                 for (Int_t i = 0; i < trackcount; i++)
 337                 {
 338                         Float_t args[7+4*3+5+4*3];
 339                         args[0] = event;
 340                         args[1] = isPrimary[i];
 341                         args[2] = pdgcode[i];
 342                         args[3] = sign[i];
 343                         args[4] = px[i];
 344                         args[5] = py[i];
 345                         args[6] = pz[i];
 346                         for (Int_t j = 0; j < 4; j++)
 347                         {
 348                                 args[7+j*3+0] = hitX[j][i];
 349                                 args[7+j*3+1] = hitY[j][i];
 350                                 args[7+j*3+2] = hitZ[j][i];
 351                         }
 352                         args[7+4*3+0] = L0sign[i];
 353                         args[7+4*3+1] = L0px[i];
 354                         args[7+4*3+2] = L0py[i];
 355                         args[7+4*3+3] = L0pz[i];
 356                         args[7+4*3+4] = L0fit[i];
 357                         for (Int_t j = 0; j < 4; j++)
 358                         {
 359                                 args[7+4*3+5+j*3+0] = L0hitX[j][i];
 360                                 args[7+4*3+5+j*3+1] = L0hitY[j][i];
 361                                 args[7+4*3+5+j*3+2] = L0hitZ[j][i];
 362                         }
 363                         table->Fill(args);
 364                 }
 365
 366         } // event loop
 367
 368         TFile file("triggertable.root", "RECREATE");
 369         file.cd();
 370         table->Write(table->GetName(), TObject::kOverwrite);
 371
 372         cout << "Done." << endl;
 373 }