HLT/MUON/macros/MakeTrackTable.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 #include "TMath.h"
  33
  34 #include "AliCDBManager.h"
  35 #include "AliLoader.h"
  36 #include "AliRunLoader.h"
  37 #include "AliStack.h"
  38 #include "AliRun.h"
  39 #include "AliMUON.h"
  40 #include "AliMUONMCDataInterface.h"
  41 #include "AliMUONHit.h"
  42 #include "AliHLTMUONEvent.h"
  43 #include "AliHLTMUONRootifierComponent.h"
  44 #include "AliHLTMUONMansoTrack.h"
  45
  46 #include <cstdlib>
  47 #include <vector>
  48 #include <iostream>
  49 using std::cout;
  50 using std::cerr;
  51 using std::endl;
  52
  53
  54 void MakeTrackTable(
  55                 Int_t firstEvent = 0,
  56                 Int_t lastEvent = -1,
  57                 const char* dHLToutputfile = "output.root",
  58                 Float_t maxSigma = 4., // 4 standard deviations
  59                 Float_t sigmaX = 0.1,  // 1 mm resolution
  60                 Float_t sigmaY = 0.01, // 100 micron resolution
  61                 Float_t sigmaZ = 0.02,  // 200 microns resolution
  62                 Float_t sigmaXtrg = 0.5,  // 5 mm resolution
  63                 Float_t sigmaYtrg = 0.5,  // 5 mm resolution
  64                 Float_t sigmaZtrg = 0.02  // 2 microns resolution
  65         )
  66 {
  67         // Setup the CDB default storage and run number if nothing was set.
  68         AliCDBManager* cdbManager = AliCDBManager::Instance();
  69         if (cdbManager == NULL)
  70         {
  71                 cerr << "ERROR: Global CDB manager object does not exist." << endl;
  72                 return;
  73         }
  74         if (cdbManager->GetDefaultStorage() == NULL)
  75         {
  76                 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
  77         }
  78         if (cdbManager->GetRun() == -1)
  79         {
  80                 cdbManager->SetRun(0);
  81         }
  82
  83         gSystem->Load("libAliHLTMUON.so");
  84         // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
  85
  86         TString fieldnames = "event:isprimary:hastrack:cantrigger:pdgcode:sign:px:py:pz";
  87         for (Int_t i = 0; i <= 14; i++)
  88         {
  89                 fieldnames += ":x";
  90                 fieldnames += i;
  91                 fieldnames += ":y";
  92                 fieldnames += i;
  93                 fieldnames += ":z";
  94                 fieldnames += i;
  95         }
  96         fieldnames += ":dHLTsign:dHLTpx:dHLTpy:dHLTpz:dHLTfit";
  97         for (Int_t i = 1; i <= 14; i++)
  98         {
  99                 fieldnames += ":dHLTx";
 100                 fieldnames += i;
 101                 fieldnames += ":dHLTy";
 102                 fieldnames += i;
 103                 fieldnames += ":dHLTz";
 104                 fieldnames += i;
 105         }
 106         TNtuple* table = new TNtuple("tracktable", "Table of tracks.", fieldnames);
 107
 108         TFile dHLTfile(dHLToutputfile, "READ");
 109
 110         AliMUONMCDataInterface data;
 111         AliStack* stack;
 112
 113         Int_t nEvents = data.NumberOfEvents();
 114         if (lastEvent < 0) lastEvent = nEvents-1;
 115         for (Int_t event = firstEvent; event <= lastEvent; event++)
 116         {
 117                 cout << "Processing event: " << event;
 118
 119                 char buf[1024];
 120                 char* str = &buf[0];
 121                 sprintf(str, "AliHLTMUONEvent;%d", event+1);
 122                 AliHLTMUONEvent* dHLTevent = dynamic_cast<AliHLTMUONEvent*>( dHLTfile.Get(str) );
 123                 if (dHLTevent == NULL)
 124                 {
 125                         cout << endl;
 126                         cerr << "ERROR: Could not find " << str << " in " << dHLToutputfile << endl;
 127                         continue;
 128                 }
 129
 130                 data.GetEvent(event);
 131
 132                 stack = data.Stack(event);
 133                 if (stack == NULL)
 134                 {
 135                         cout << endl;
 136                         cerr << "ERROR: Stack is NULL" << endl;
 137                         continue;
 138                 }
 139
 140                 Int_t trackcount = data.NumberOfParticles();
 141                 cout << "\ttrackcount = " << trackcount << endl;
 142
 143                 TArrayF isPrimary(trackcount);
 144                 TArrayF hasTrack(trackcount);
 145                 TArrayI hitsInTrigger(trackcount);
 146                 TArrayF pdgcode(trackcount);
 147                 TArrayF sign(trackcount);
 148                 TArrayF px(trackcount);
 149                 TArrayF py(trackcount);
 150                 TArrayF pz(trackcount);
 151                 TArrayF vx(trackcount);
 152                 TArrayF vy(trackcount);
 153                 TArrayF vz(trackcount);
 154                 TArrayF dHLTsign(trackcount);
 155                 TArrayF dHLTpx(trackcount);
 156                 TArrayF dHLTpy(trackcount);
 157                 TArrayF dHLTpz(trackcount);
 158                 TArrayF dHLTfit(trackcount);
 159
 160                 //cout << "Fill particle data" << endl;
 161                 for (Int_t i = 0; i < trackcount; i++)
 162                 {
 163                         TParticle* p = data.Particle(i);
 164                         //isPrimary[i] = (i < stack->GetNprimary()) ? 1 : 0;
 165                         isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
 166                         hasTrack[i] = 0;
 167                         hitsInTrigger[i] = 0;
 168                         pdgcode[i] = p->GetPdgCode();
 169                         if (p->GetPDG() != NULL)
 170                         {
 171                                 Double_t charge = p->GetPDG()->Charge();
 172                                 if (charge > 0)
 173                                         sign[i] = 1;
 174                                 else if (charge < 0)
 175                                         sign[i] = -1;
 176                                 else
 177                                         sign[i] = 0;
 178                         }
 179                         else
 180                         {
 181                                 sign[i] = 0;
 182                                 cerr << "ERROR: Could not get the PDG information." << endl;
 183                         }
 184                         px[i] = p->Px();
 185                         py[i] = p->Py();
 186                         pz[i] = p->Pz();
 187                         vx[i] = p->Vx();
 188                         vy[i] = p->Vy();
 189                         vz[i] = p->Vz();
 190                         dHLTsign[i] = 0;
 191                         dHLTpx[i] = 0;
 192                         dHLTpy[i] = 0;
 193                         dHLTpz[i] = 0;
 194                         dHLTfit[i] = -1;
 195                 }
 196
 197                 TArrayF hitX[14];
 198                 TArrayF hitY[14];
 199                 TArrayF hitZ[14];
 200                 TArrayF dHLThitX[14];
 201                 TArrayF dHLThitY[14];
 202                 TArrayF dHLThitZ[14];
 203                 for (Int_t i = 0; i < 14; i++)
 204                 {
 205                         hitX[i].Set(trackcount);
 206                         hitY[i].Set(trackcount);
 207                         hitZ[i].Set(trackcount);
 208                         dHLThitX[i].Set(trackcount);
 209                         dHLThitY[i].Set(trackcount);
 210                         dHLThitZ[i].Set(trackcount);
 211                         for (Int_t j = 0; j < trackcount; j++)
 212                         {
 213                                 hitX[i][j] = 0;
 214                                 hitY[i][j] = 0;
 215                                 hitZ[i][j] = 0;
 216                                 dHLThitX[i][j] = 0;
 217                                 dHLThitY[i][j] = 0;
 218                                 dHLThitZ[i][j] = 0;
 219                         }
 220                 }
 221
 222                 //cout << "Fill hits" << endl;
 223                 for (Int_t i = 0; i < data.NumberOfTracks(); i++)
 224                 for (Int_t j = 0; j < data.NumberOfHits(i); j++)
 225                 {
 226                         AliMUONHit* hit = data.Hit(i, j);
 227                         Int_t chamber = hit->Chamber() - 1;
 228                         if (chamber < 0 || chamber >= 14)
 229                         {
 230                                 cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
 231                                 continue;
 232                         }
 233
 234                         //cout << "hit->Track() = " << hit->Track() << endl;
 235                         if (hit->Track() < 0 || hit->Track() >= trackcount)
 236                         {
 237                                 cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
 238                                         << trackcount << ")" << endl;
 239                                 continue;
 240                         }
 241                         Int_t particleindex = hit->Track();
 242
 243                         hitX[chamber][particleindex] = hit->X();
 244                         hitY[chamber][particleindex] = hit->Y();
 245                         hitZ[chamber][particleindex] = hit->Z();
 246
 247                         if (chamber == 0)
 248                         {
 249                                 hasTrack[particleindex] = 1;
 250                         }
 251                         if (chamber >= 10)
 252                         {
 253                                 hitsInTrigger[particleindex]++;
 254                         }
 255                 } // hits loop
 256
 257                 // Ok now go through the dHLT information for this event and correlate
 258                 // the dHLT tracks to the kine data.
 259                 // This is done as follows: take a dHLT track. Then for each kine track calculate
 260                 // the distance between each hit of the kine track and corresponding dHLT track.
 261                 // Then sum up the distances and divide by the number of hits matched.
 262                 // This gives us a fit quality /Chi2 type parameter. The kine track that
 263                 // has the smallest sum of differences is then chosen as the kine track
 264                 // correlating to the current dHLT track being processed.
 265                 for (Int_t i = 0; i < dHLTevent->Array().GetEntriesFast(); i++)
 266                 {
 267                         if (dHLTevent->Array().At(i)->IsA() != AliHLTMUONMansoTrack::Class())
 268                                 continue;
 269                         AliHLTMUONMansoTrack* mtrack = static_cast<AliHLTMUONMansoTrack*>(dHLTevent->Array().At(i));
 270
 271                         // Find the best fitting kine track.
 272                         Double_t bestFitQuality = 1e30;
 273                         Int_t bestFitTrack = -1;
 274                         for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
 275                         {
 276                                 Double_t sumOfDiffs = 0;
 277                                 Double_t hitsMatched = 0;
 278                                 for (Int_t j = 7; j <= 10; j++)
 279                                 {
 280                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j);
 281                                         if (hit == NULL) continue;
 282                                         if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
 283                                                 continue;
 284
 285                                         TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
 286                                         TVector3 diff = hV - hit->Coordinate();
 287                                         Double_t diffX = diff.X() / sigmaX;
 288                                         Double_t diffY = diff.Y() / sigmaY;
 289                                         Double_t diffZ = diff.Z() / sigmaZ;
 290                                         if (diffX > maxSigma) continue;
 291                                         if (diffY > maxSigma) continue;
 292                                         if (diffZ > maxSigma) continue;
 293                                         sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
 294                                         hitsMatched++;
 295                                 }
 296                                 if (mtrack->TriggerRecord() != NULL)
 297                                 {
 298                                         for (Int_t j = 11; j <= 14; j++)
 299                                         {
 300                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
 301                                                 if (hit == TVector3(0,0,0)) continue;
 302                                                 if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
 303                                                         continue;
 304
 305                                                 TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
 306                                                 TVector3 diff = hV - hit;
 307                                                 Double_t diffX = diff.X() / sigmaXtrg;
 308                                                 Double_t diffY = diff.Y() / sigmaYtrg;
 309                                                 Double_t diffZ = diff.Z() / sigmaZtrg;
 310                                                 if (diffX > maxSigma) continue;
 311                                                 if (diffY > maxSigma) continue;
 312                                                 if (diffZ > maxSigma) continue;
 313                                                 sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
 314                                                 hitsMatched++;
 315                                         }
 316                                 }
 317
 318                                 // Now check the fit quality.
 319                                 if (hitsMatched <= 0) continue;
 320                                 Double_t fitQuality = TMath::Sqrt(sumOfDiffs) / hitsMatched;
 321                                 if (fitQuality < bestFitQuality)
 322                                 {
 323                                         bestFitQuality = fitQuality;
 324                                         bestFitTrack = kinetrack;
 325                                 }
 326                         }
 327
 328                         if (bestFitTrack < 0)
 329                         {
 330                                 // No track was matched so we need to add a fake track record.
 331                                 Float_t args[12+14*3+5+14*3];
 332                                 args[0] = event;
 333                                 args[1] = -1;
 334                                 args[2] = 0;
 335                                 args[3] = 0;
 336                                 args[4] = 0;
 337                                 args[5] = 0;
 338                                 args[6] = 0;
 339                                 args[7] = 0;
 340                                 args[8] = 0;
 341                                 args[9] = 0;
 342                                 args[10] = 0;
 343                                 args[11] = 0;
 344                                 for (Int_t j = 0; j < 14; j++)
 345                                 {
 346                                         args[12+j*3+0] = 0;
 347                                         args[12+j*3+1] = 0;
 348                                         args[12+j*3+2] = 0;
 349                                 }
 350                                 args[12+14*3+0] = mtrack->Sign();
 351                                 args[12+14*3+1] = mtrack->Px();
 352                                 args[12+14*3+2] = mtrack->Py();
 353                                 args[12+14*3+3] = mtrack->Pz();
 354                                 args[12+14*3+4] = -1;
 355                                 for (Int_t j = 0; j < 14; j++)
 356                                 {
 357                                         args[12+14*3+5+j*3+0] = 0;
 358                                         args[12+14*3+5+j*3+1] = 0;
 359                                         args[12+14*3+5+j*3+2] = 0;
 360                                 }
 361                                 for (Int_t j = 6; j < 10; j++)
 362                                 {
 363                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j+1);
 364                                         if (hit == NULL) continue;
 365                                         args[12+14*3+5+j*3+0] = hit->X();
 366                                         args[12+14*3+5+j*3+1] = hit->Y();
 367                                         args[12+14*3+5+j*3+2] = hit->Z();
 368                                 }
 369                                 if (mtrack->TriggerRecord() != NULL)
 370                                 {
 371                                         for (Int_t j = 10; j < 14; j++)
 372                                         {
 373                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j+1);
 374                                                 args[12+14*3+5+j*3+0] = hit.X();
 375                                                 args[12+14*3+5+j*3+1] = hit.Y();
 376                                                 args[12+14*3+5+j*3+2] = hit.Z();
 377                                         }
 378                                 }
 379                                 table->Fill(args);
 380                         }
 381                         else
 382                         {
 383                                 // Fill the details about the dHLT track to the best fitting track info.
 384                                 dHLTsign[bestFitTrack] = mtrack->Sign();
 385                                 dHLTpx[bestFitTrack] = mtrack->Px();
 386                                 dHLTpy[bestFitTrack] = mtrack->Py();
 387                                 dHLTpz[bestFitTrack] = mtrack->Pz();
 388                                 dHLTfit[bestFitTrack] = bestFitQuality;
 389                                 for (Int_t j = 7; j <= 10; j++)
 390                                 {
 391                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j);
 392                                         if (hit == NULL) continue;
 393                                         dHLThitX[j-1][bestFitTrack] = hit->X();
 394                                         dHLThitY[j-1][bestFitTrack] = hit->Y();
 395                                         dHLThitZ[j-1][bestFitTrack] = hit->Z();
 396                                 }
 397                                 if (mtrack->TriggerRecord() != NULL)
 398                                 {
 399                                         for (Int_t j = 11; j <= 14; j++)
 400                                         {
 401                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
 402                                                 if (hit == TVector3(0,0,0)) continue;
 403                                                 dHLThitX[j-1][bestFitTrack] = hit.X();
 404                                                 dHLThitY[j-1][bestFitTrack] = hit.Y();
 405                                                 dHLThitZ[j-1][bestFitTrack] = hit.Z();
 406                                         }
 407                                 }
 408                         }
 409                 }
 410
 411                 //cout << "Fill table" << endl;
 412                 for (Int_t i = 0; i < trackcount; i++)
 413                 {
 414                         Float_t args[12+14*3+5+14*3];
 415                         args[0] = event;
 416                         args[1] = isPrimary[i];
 417                         args[2] = hasTrack[i];
 418                         args[3] = (hitsInTrigger[i] > 2) ? 1 : 0;
 419                         args[4] = pdgcode[i];
 420                         args[5] = sign[i];
 421                         args[6] = px[i];
 422                         args[7] = py[i];
 423                         args[8] = pz[i];
 424                         args[9] = vx[i];
 425                         args[10] = vy[i];
 426                         args[11] = vz[i];
 427                         for (Int_t j = 0; j < 14; j++)
 428                         {
 429                                 args[12+j*3+0] = hitX[j][i];
 430                                 args[12+j*3+1] = hitY[j][i];
 431                                 args[12+j*3+2] = hitZ[j][i];
 432                         }
 433                         args[12+14*3+0] = dHLTsign[i];
 434                         args[12+14*3+1] = dHLTpx[i];
 435                         args[12+14*3+2] = dHLTpy[i];
 436                         args[12+14*3+3] = dHLTpz[i];
 437                         args[12+14*3+4] = dHLTfit[i];
 438                         for (Int_t j = 0; j < 14; j++)
 439                         {
 440                                 args[12+14*3+5+j*3+0] = dHLThitX[j][i];
 441                                 args[12+14*3+5+j*3+1] = dHLThitY[j][i];
 442                                 args[12+14*3+5+j*3+2] = dHLThitZ[j][i];
 443                         }
 444                         table->Fill(args);
 445                 }
 446
 447         } // event loop
 448
 449         TFile file("tracktable.root", "RECREATE");
 450         file.cd();
 451         table->Write(table->GetName(), TObject::kOverwrite);
 452
 453         cout << "Done." << endl;
 454 }