HLT/ITS/trackingSAP/Process.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TH1F.h>
   3 #include <TH1I.h>
   4 #include <TH2D.h>
   5 #include <TFile.h>
   6 #include <Riostream.h>
   7 #include <TParticle.h>
   8 #include <TCanvas.h>
   9 #include <TLegend.h>
  10 #include <TSystem.h>
  11 #include <TROOT.h>
  12 #include <TClonesArray.h>
  13 #include <TTree.h>
  14 #include <TArrayF.h>
  15 #include <TStyle.h>
  16 #include <TString.h>
  17 #include <TStopwatch.h>
  18 #include <TParticlePDG.h>
  19 #include "AliRun.h"
  20 #include "AliRunLoader.h"
  21 #include "AliMultiplicity.h"
  22 #include "AliESDtrackCuts.h"
  23 #include "AliStack.h"
  24 #include "AliESDEvent.h"
  25 #include "AliESDtrack.h"
  26 #include "AliITSRecPointContainer.h"
  27 #include "AliITSgeomTGeo.h"
  28 #include "AliITSRecPoint.h"
  29 #include "AliCDBManager.h"
  30 #include "AliHeader.h"
  31 #include "AliGenEventHeader.h"
  32 //
  33
  34 #include "../HLT/ITS/trackingSAP/AliITSSAPAux.h"
  35 #include "../HLT/ITS/trackingSAP/AliITSSAPLayer.h"
  36 #include "../HLT/ITS/trackingSAP/AliITSSAPTracker.h"
  37 #include "TProfile.h"
  38 #include "TStopwatch.h"
  39 #endif
  40
  41
  42 AliITSSAPTracker* tracker=0;
  43
  44 TTree* treeInp = 0;
  45 AliRunLoader* runLoader = 0;
  46 AliESDEvent* esd=0;
  47 typedef struct {
  48   Bool_t  validSPD;
  49   Bool_t  validTrc;
  50   Bool_t  isPrim;
  51   Char_t  ntlC;
  52   Char_t  ntlF;
  53   Char_t  mtlC;
  54   Char_t  mtlF;
  55   Char_t  ntrC;
  56   Char_t  ntrF;
  57   Int_t   pdg;
  58   Int_t   evID;
  59   Int_t   mult;
  60   Float_t pt;
  61   Float_t eta;
  62   Float_t phi;
  63   Float_t zv;
  64   //
  65   Float_t spdDPhi;
  66   Float_t spdDTht;
  67   Float_t spdChi2;
  68 } tres_t;
  69
  70 tres_t tres;
  71 TFile* flOut = 0;
  72 TTree* trOut = 0;
  73
  74
  75 void ProcessEvent(int iev);
  76 void Process(const char* path);
  77 void ProcChunk(const char* path);
  78 void TestTracker(TTree* tRP, const AliESDVertex* vtx);
  79 void LoadClusters(TTree* tRP);
  80 Double_t* DefLogAx(double xMn,double xMx, int nbin);
  81 void CheckRecStatus();
  82 void InitOutTree(const char* fname="rcInfo.root");
  83 void SaveOutTree();
  84 void InitTracker(int runNumber);
  85
  86
  87
  88 #ifdef _TIMING_
  89 TProfile* hTiming[AliITSSAPTracker::kNSW];
  90 #endif
  91
  92 void Process(const char* inpData)
  93 {
  94   //
  95   InitOutTree();
  96   TString inpDtStr = inpData;
  97   if (inpDtStr.EndsWith(".root") || inpDtStr.EndsWith(".zip#")) {
  98     ProcChunk(inpDtStr.Data());
  99   }
 100   else {
 101     ifstream inpf(inpData);
 102     if (!inpf.good()) {
 103       printf("Failed on input filename %s\n",inpData);
 104       return;
 105     }
 106     //
 107     inpDtStr.ReadLine(inpf);
 108     while ( !inpDtStr.IsNull() ) {
 109       inpDtStr = inpDtStr.Strip(TString::kBoth,' ');
 110       if (inpDtStr.BeginsWith("//") || inpDtStr.BeginsWith("#")) {inpDtStr.ReadLine(inpf); continue;}
 111       inpDtStr = inpDtStr.Strip(TString::kBoth,',');
 112       inpDtStr = inpDtStr.Strip(TString::kBoth,'"');
 113       ProcChunk(inpDtStr.Data());
 114       inpDtStr.ReadLine(inpf);
 115     }
 116   }
 117 #ifdef _CONTROLH_
 118   tracker->SaveHistos();
 119 #endif
 120   SaveOutTree();
 121 }
 122
 123 void ProcChunk(const char* path)
 124 {
 125   //
 126   TString dir=path;
 127   //
 128   printf("Processing %s\n",dir.Data());
 129   //
 130   if (dir.EndsWith("AliESDs.root")) {
 131     dir.ReplaceAll("AliESDs.root","");
 132   }
 133   //
 134   esd = new AliESDEvent();
 135   //
 136   runLoader = AliRunLoader::Open(Form("%sgalice.root",dir.Data()));
 137   if (!runLoader) {
 138     printf("galice not found\n");
 139     return;
 140   }
 141   runLoader->LoadgAlice();
 142   gAlice = runLoader->GetAliRun();
 143   runLoader->LoadHeader();
 144   runLoader->LoadKinematics();
 145   runLoader->LoadRecPoints("ITS");
 146   // ESD
 147   TFile* esdFile = TFile::Open(Form("%sAliESDs.root",dir.Data()));
 148   if (!esdFile || !esdFile->IsOpen()) {
 149     printf("Error in opening ESD file\n");
 150     runLoader->UnloadKinematics();
 151     runLoader->UnloadHeader();
 152     runLoader->UnloadgAlice();
 153     delete runLoader;
 154     return;
 155   }
 156   treeInp = (TTree*) esdFile->Get("esdTree");
 157   if (!treeInp) {
 158     printf("Error: no ESD tree found\n");
 159     runLoader->UnloadKinematics();
 160     runLoader->UnloadHeader();
 161     runLoader->UnloadgAlice();
 162     delete runLoader;
 163     return;
 164   }
 165   esd->ReadFromTree(treeInp);
 166   //
 167   for (int iEv=0;iEv<runLoader->GetNumberOfEvents(); iEv++) {
 168     printf("ev %d\n",iEv);
 169     ProcessEvent(iEv);
 170   }
 171   runLoader->UnloadRecPoints("all");
 172   runLoader->UnloadKinematics();
 173   runLoader->UnloadHeader();
 174   runLoader->UnloadgAlice();
 175   delete runLoader;
 176   runLoader = 0;
 177   esdFile->Close();
 178   delete esdFile;
 179 }
 180
 181 //_______________________________________________
 182 void ProcessEvent(int iEv)
 183 {
 184   runLoader->GetEvent(iEv);
 185   treeInp->GetEvent(iEv);
 186   //
 187   if (!tracker) InitTracker(esd->GetRunNumber());
 188   //
 189   AliStack *stack = runLoader->Stack();
 190   Int_t nParticles = stack ? stack->GetNtrack() : 0;
 191   Int_t nTracks= esd->GetNumberOfTracks();
 192   printf("Ntr: %d NPartMC: %d\n",nTracks,nParticles);
 193   const AliESDVertex* vtx = esd->GetPrimaryVertexSPD();
 194   if (!vtx || !vtx->GetStatus()) return;
 195   //
 196   TTree* treeITSRP = runLoader->GetTreeR("ITS",kFALSE);
 197   if (!treeITSRP) {
 198     printf("Failed to fetch ITS recpoints\n");
 199     exit(1);
 200   }
 201   //
 202 #ifdef _TIMING_
 203   static Double_t cpuTime[AliITSSAPTracker::kNSW] = {0};
 204 #endif
 205   printf("\n\n\nEvent: %d\n",iEv);
 206   TestTracker(treeITSRP, vtx);
 207   //
 208 #ifdef _TIMING_
 209   if (vtx->GetStatus()==1) {
 210     double mlt = esd->GetMultiplicity()->GetNumberOfTracklets();
 211     for (int i=0;i<AliITSSAPTracker::kNSW;i++) {
 212       TStopwatch& sw = (TStopwatch&)tracker->GetStopwatch(i);
 213       double tm = sw.CpuTime();
 214       hTiming[i]->Fill(mlt, tm - cpuTime[i]);
 215       cpuTime[i] = tm;
 216     }
 217   }
 218 #endif
 219   //
 220   CheckRecStatus();
 221   //esd->Reset();
 222   //
 223 }
 224
 225 //_________________________________________________
 226 void TestTracker(TTree* tRP, const AliESDVertex* vtx)
 227 {
 228   tracker->Clear();
 229   LoadClusters(tRP);
 230   tracker->SetSPDVertex(vtx);
 231   tracker->SetBz(esd->GetMagneticField());
 232   tracker->ProcessEvent();
 233   //
 234   tracker->PrintTracklets();
 235   tracker->PrintTracks();
 236   //
 237   //
 238   //  vtx->Print();
 239   //  esd->GetMultiplicity()->Print("t");
 240   //
 241   /*
 242   AliHeader* hd = runLoader->GetHeader();
 243   AliGenEventHeader* hdmc;
 244   if (tracker->GetTrackVertex().GetStatus()==1 && hd && (hdmc=hd->GenEventHeader()) ) {
 245     TArrayF vtxMC;
 246     hdmc->PrimaryVertex(vtxMC);
 247     printf("MCvtx: %f %f %f\n",vtxMC[0],vtxMC[1],vtxMC[2]);
 248   }
 249   */
 250 }
 251
 252 //_________________________________________________
 253 void LoadClusters(TTree* tRP)
 254 {
 255   //  AliITSRecPointContainer* rpcont = AliITSRecPointContainer::Instance();
 256   //  TClonesArray* itsClusters = rpcont->FetchClusters(0,tRP);
 257   //  if(!rpcont->IsSPDActive()){
 258   //    printf("No SPD rec points found, multiplicity not calculated\n");
 259   //    tRP->Print();
 260   //    return;
 261   //  }
 262   //  else printf("NSP0: %d\n",itsClusters->GetEntriesFast());
 263   static TClonesArray* itsModAdd[2198] = {0};
 264   static Bool_t first = kTRUE;
 265   if (first) {
 266     first = 0;
 267     for (int i=0;i<2198;i++) itsModAdd[i] = new TClonesArray("AliITSRecPoint");
 268   }
 269   int nMod = AliITSgeomTGeo::GetNModules();
 270   for (int imd=0;imd<nMod;imd++) {
 271     TClonesArray* itsClusters = itsModAdd[imd];
 272     tRP->SetBranchAddress("ITSRecPoints",&itsClusters);
 273     tRP->GetEntry(imd);
 274     //    itsClusters = rpcont->UncheckedGetClusters(imd);
 275
 276     int nClusters = itsClusters->GetEntriesFast();
 277     if (!nClusters) continue;
 278     while(nClusters--) {
 279       AliITSRecPoint* cluster = (AliITSRecPoint*)itsClusters->UncheckedAt(nClusters);
 280       if (!cluster) continue;
 281       tracker->AddCluster(cluster);
 282     }
 283   }
 284   //
 285 }
 286
 287 //______________________________________________
 288 Double_t* DefLogAx(double xMn,double xMx, int nbin)
 289 {
 290   // get array for log axis
 291   if (xMn<=0 || xMx<=xMn || nbin<2) {
 292     printf("Wrong axis request: %f %f %d\n",xMn,xMx,nbin);
 293     return 0;
 294   }
 295   double dx = log(xMx/xMn)/nbin;
 296   double *xax = new Double_t[nbin+1];
 297   for (int i=0;i<=nbin;i++) xax[i]= xMn*exp(dx*i);
 298   return xax;
 299 }
 300 //______________________________________________
 301 void CheckRecStatus()
 302 {
 303   //
 304   static int nev = -1;
 305   static TBits patternMC;
 306   static vector<char> ntrCorr;
 307   static vector<char> ntrFake;
 308   static vector<char> mtlCorr;
 309   static vector<char> mtlFake;
 310   static vector<char> ntlCorr;
 311   static vector<char> ntlFake;
 312   static vector<float> spdDPhi;
 313   static vector<float> spdDTht;
 314   static vector<float> spdChi2;
 315   //
 316   AliStack* stack = 0;
 317   AliRunLoader* rl = AliRunLoader::Instance();
 318   if (!rl || !(stack=rl->Stack())) return;
 319   nev++;
 320   //
 321   enum {kIsPrim=AliITSSAPTracker::kNLrActive,kValidTracklet,kValidTrack,kRecDone,kBitPerTrack};
 322   int nTrkMC = stack->GetNtrack();
 323   patternMC.SetBitNumber(nTrkMC*kBitPerTrack,0);
 324   patternMC.ResetAllBits();
 325   //
 326   ntrCorr.clear();
 327   ntrFake.clear();
 328   ntlCorr.clear();
 329   ntlFake.clear();
 330   mtlCorr.clear();
 331   mtlFake.clear();
 332   //
 333   spdDPhi.clear();
 334   spdDTht.clear();
 335   spdChi2.clear();
 336   //
 337   ntrCorr.resize(nTrkMC,0);
 338   ntrFake.resize(nTrkMC,0);
 339   ntlCorr.resize(nTrkMC,0);
 340   ntlFake.resize(nTrkMC,0);
 341   mtlCorr.resize(nTrkMC,0);
 342   mtlFake.resize(nTrkMC,0);
 343   spdDPhi.resize(nTrkMC,0);
 344   spdDTht.resize(nTrkMC,0);
 345   spdChi2.resize(nTrkMC,0);
 346
 347   //
 348   // fill MC track patterns
 349   for (int ilr=6;ilr--;) {
 350     AliITSSAPLayer *lr = tracker->GetLayer(ilr);
 351     int ncl = lr->GetNClusters();
 352     for (int icl=ncl;icl--;) {
 353       AliITSRecPoint* cl = lr->GetClusterUnSorted(icl);
 354       for (int j=0;j<3;j++) {
 355         int lb = cl->GetLabel(j);
 356         if (lb<0 || lb>=nTrkMC) break;
 357         patternMC.SetBitNumber(lb*kBitPerTrack+ilr,kTRUE);
 358       }
 359     }
 360   }
 361   //
 362   int nTrk = tracker->GetNTracklets();
 363   if (!nTrk) return;
 364   for (int itr=0;itr<nTrk;itr++) {
 365     const AliITSSAPTracker::SPDtracklet_t& trlet = tracker->GetTracklet(itr);
 366     //    PrintTracklet(itr);
 367     //
 368     int lbl = trlet.label;
 369     if (lbl==-3141593) continue;
 370     int lblA = TMath::Abs(lbl);
 371     if (lblA==lbl) ntlCorr[lblA]++;
 372     else           ntlFake[lblA]++;
 373     //
 374     if (spdChi2[lblA]==0 || spdChi2[lblA]<trlet.chi2) {
 375       spdChi2[lblA] = trlet.chi2;
 376       spdDPhi[lblA] = trlet.dphi;
 377       spdDTht[lblA] = trlet.dtht;
 378     }
 379   }
 380   //
 381   AliITSSAPLayer* lr0 = tracker->GetLayer(0);
 382   AliITSSAPLayer* lr1 = tracker->GetLayer(1);
 383   for (int itrm=0;itrm<nTrkMC;itrm++) {
 384     if (ntlCorr[itrm]+ntlFake[itrm]<2) continue;
 385     printf("\nExtra for tr %d nC:%d nF:%d\n",itrm,ntlCorr[itrm],ntlFake[itrm]);
 386     //
 387     int cnt = 0;
 388     for (int itr=0;itr<nTrk;itr++) {
 389       const AliITSSAPTracker::SPDtracklet_t& trlet = tracker->GetTracklet(itr);
 390       if (TMath::Abs(trlet.label)!=itrm) continue;
 391       AliITSSAPLayer::ClsInfo_t* clinf0 = lr0->GetClusterInfo(trlet.id1);
 392       AliITSSAPLayer::ClsInfo_t* clinf1 = lr1->GetClusterInfo(trlet.id2);
 393       printf("#%2d%s%4d chi:%.2f [%4d/%3d] [%4d/%3d]\n",cnt++,trlet.label<0 ? "-":"+",itr,trlet.chi2,
 394              trlet.id1,clinf0->detid,
 395              trlet.id2,clinf1->detid);
 396     }
 397   }
 398   //
 399   const AliMultiplicity* mltESD = esd->GetMultiplicity();
 400   nTrk = mltESD->GetNumberOfTracklets();
 401   for (int itr=0;itr<nTrk;itr++) {
 402     int lb0 = mltESD->GetLabel(itr,0);
 403     int lb1 = mltESD->GetLabel(itr,1);
 404     if (lb0==lb1) mtlCorr[lb1]++;
 405     else          mtlFake[lb1]++;
 406   }
 407   //
 408   nTrk = tracker->GetNTracks();
 409   for (int itr=0;itr<nTrk;itr++) {
 410     const AliITSSAPTracker::ITStrack_t &track = tracker->GetTrack(itr);
 411     //
 412     int lbl = track.label;
 413     if (lbl==-3141593) continue;
 414     int lblA = TMath::Abs(lbl);
 415     if (lblA==lbl) ntrCorr[lblA]++;
 416     else           ntrFake[lblA]++;
 417   }
 418   //
 419   // set reconstructability
 420   for (int itr=nTrkMC;itr--;) {
 421     int bitoffs = itr*kBitPerTrack;
 422     //
 423     tres.validSPD = patternMC.TestBitNumber(bitoffs+AliITSSAPTracker::kALrSPD1) &&
 424       patternMC.TestBitNumber(bitoffs+AliITSSAPTracker::kALrSPD2);
 425     int nmiss = 0;
 426     for (int il=AliITSSAPTracker::kALrSDD1;il<=AliITSSAPTracker::kALrSSD2;il++)
 427       if (tracker->IsObligatoryLayer(il) && !patternMC.TestBitNumber(bitoffs+il)) nmiss++;
 428     tres.validTrc = tres.validSPD && (nmiss<=tracker->GetMaxMissedLayers());
 429     //
 430     if ( !tres.validSPD && !tres.validTrc &&
 431          (ntlCorr[itr]+ntlFake[itr])==0 &&
 432          (ntrCorr[itr]+ntrFake[itr])==0 &&
 433          (mtlCorr[itr]+mtlFake[itr])==0 ) continue;
 434     //
 435     TParticle* mctr = stack->Particle(itr);
 436     tres.isPrim = stack->IsPhysicalPrimary(itr);
 437     tres.pdg = mctr->GetPdgCode();
 438     tres.pt =  mctr->Pt();
 439     tres.eta = mctr->Eta();
 440     tres.phi = mctr->Phi();
 441     //
 442     tres.ntlC = ntlCorr[itr];
 443     tres.ntlF = ntlFake[itr];
 444     tres.ntrC = ntrCorr[itr];
 445     tres.ntrF = ntrFake[itr];
 446     //
 447     tres.mtlC = mtlCorr[itr];
 448     tres.mtlF = mtlFake[itr];
 449     //
 450     tres.evID = nev;
 451     tres.mult = tracker->GetNTracklets();
 452     tres.zv = tracker->GetSPDVertex()->GetZ();
 453     //
 454     tres.spdDPhi = spdDPhi[itr];
 455     tres.spdDTht = spdDTht[itr];
 456     tres.spdChi2 = spdChi2[itr];
 457     //
 458     trOut->Fill();
 459   }
 460   //
 461 }
 462
 463 //___________________________________________________
 464 void InitOutTree(const char* fname)
 465 {
 466   // output tree
 467   flOut = TFile::Open(fname,"recreate");
 468   trOut = new TTree("rcinfo","rcinfo");
 469   trOut->Branch("res",&tres);
 470 }
 471
 472 //___________________________________________________
 473 void SaveOutTree()
 474 {
 475   if (!trOut) return;
 476   flOut->cd();
 477   trOut->Write();
 478   delete trOut;
 479   flOut->Close();
 480   delete flOut;
 481   trOut = 0;
 482   flOut = 0;
 483 }
 484
 485
 486 //______________________________________________
 487 void InitTracker(int runNumber)
 488 {
 489   AliCDBManager* man = AliCDBManager::Instance();
 490   man->SetRun(runNumber);
 491   tracker = new AliITSSAPTracker();
 492   tracker->Init();
 493   //
 494   //
 495 #ifdef _TIMING_
 496   int nbMlt = 100;
 497   double *mltAx = DefLogAx(1,6000,nbMlt);
 498   for (int i=0;i<AliITSSAPTracker::kNSW;i++) {
 499     hTiming[i] = new TProfile(tracker->GetStopwatchName(i),tracker->GetStopwatchName(i),nbMlt,mltAx);
 500   }
 501 #endif
 502   //
 503 }