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