macro for cluster merger efficiency
authorsgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Dec 2011 06:07:21 +0000 (06:07 +0000)
committersgorbuno <sgorbuno@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Dec 2011 06:07:21 +0000 (06:07 +0000)
HLT/TPCLib/HWCFemulator/macro/ClusterMergerQA.C [new file with mode: 0644]

diff --git a/HLT/TPCLib/HWCFemulator/macro/ClusterMergerQA.C b/HLT/TPCLib/HWCFemulator/macro/ClusterMergerQA.C
new file mode 100644 (file)
index 0000000..596fdda
--- /dev/null
@@ -0,0 +1,319 @@
+/****************************************************************************
+ *           Very important, delicate and rather obscure macro.             *
+ *                                                                          *
+ *               Creates list of "trackable" tracks,                        *
+ *             calculates efficiency, resolutions etc.                      *
+ *     There is a possibility to run this macro over several events.         *
+ *                                                                          *
+ *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
+ * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de        *
+ ****************************************************************************/
+
+#if defined(__CINT__) && !defined(__MAKECINT__)
+{
+  
+  gSystem->AddIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT -I$ALICE_ROOT/TPC -I$ALICE/geant3/TGeant3 -I$ALICE_ROOT/HLT/TPCLib/HWCFemulator -I$ALICE_ROOT/HLT/TPCLib -I$ALICE_ROOT/HLT/BASE ");    
+  AliHLTPluginBase::GetInstance()->LoadComponentLibraries("libAliHLTTPC.so");    
+  gROOT->LoadMacro("ClusterMergerQA.C++");
+
+  runTest("NoMerge/TPC.RecPoints.root");
+
+}
+#else
+
+  #include <TMath.h>
+  #include <TError.h>
+  #include <Riostream.h>
+  #include <TH1F.h>
+  #include <TH2F.h>
+  #include <TTree.h>
+  #include <TParticle.h>
+  #include <TCanvas.h>
+  #include <TLine.h>
+  #include <TText.h>
+  #include <TStyle.h>
+  #include <TFile.h>
+  #include <TROOT.h>
+  #include "TString.h"
+  #include "AliTPCclusterMI.h"
+
+  #include "AliStack.h"
+  #include "AliHeader.h"
+  #include "AliTrackReference.h"
+  #include "AliRunLoader.h"
+  #include "AliRun.h"
+  #include "AliTPCtrack.h"
+  #include "AliTracker.h"
+  
+  #include "AliTPC.h"
+  #include "AliTPCClustersArray.h"
+  #include "AliTPCClustersRow.h"
+  #include "AliTPCcluster.h"
+  #include "AliTPCLoader.h"
+  #include "TParticlePDG.h"
+  #include "TDatabasePDG.h"
+  #include "AliGeomManager.h"
+  #include <iostream>
+  #include <string>
+  #include "AliCDBManager.h"
+  #include "AliHLTTPCHWClusterMerger.h"
+  #include "AliHLTTPCTransform.h"
+
+extern AliRun *gAlice;
+extern TROOT *gROOT;
+
+
+void myDraw(TH1F *h)
+{
+  Int_t minc=33; 
+  if (h->GetEntries()<minc) h->Draw(); else h->Fit("gaus","Q");
+}
+
+Int_t runTest(const char *recPointsName ) 
+{
+
+  // access clusters
+  
+  //const Char_t* recPointsName = "NoMerge/TPC.RecPoints.root";  
+
+  AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
+  if (!rl) return 1;
+  rl->LoadgAlice();
+  rl->LoadHeader();
+
+  AliCDBManager *man=AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+  //man->SetDefaultStorage("local:///lustre/alice/alien/alice/simulation/2008/v4-15-Release/Residual/");
+  man->SetSpecificStorage("GRP/GRP/Data",
+                         Form("local://%s",gSystem->pwd()));
+  man->SetRun(0);
+  //man->SetRun(137366);
+  
+  if(AliGeomManager::GetGeometry()==NULL){
+    AliGeomManager::LoadGeometry();
+  }
+
+  AliGeomManager::ApplyAlignObjsFromCDB("TPC");
+  
+  //TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG, AliMagF::kBeamTypeAA, 1300.));
+  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
+  // Check field
+  if (!TGeoGlobalMagField::Instance()) {
+    ::Error("","magnetic field not initialized, please set up TGeoGlobalMagField and AliMagF");
+    return 1;
+  }
+  
+
+  Int_t nev=rl->GetNumberOfEvents();
+  // ******* Loop over events *********
+
+
+  AliMagF *field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+
+  if(! field) {
+    ::Error("","No mag. field found!!!");
+    return 1;
+  }
+
+  AliHLTTPCHWClusterMerger merger;
+  merger.Init();
+
+  int statNGood=0, statNMerged=0;
+
+  //nev = 1;
+  // ********  Loop over generated events 
+  for (Int_t nEv=0; nEv<nev; nEv++) {
+    cout<<"Event " <<nEv<<endl;
+    
+    merger.Clear();
+    static int nClusters =0;
+    static int nAccepted = 0;
+
+    vector<AliTPCclusterMI> vClusters;
+
+    // Read clusters
+    { 
+      rl->GetEvent(nEv);
+    
+      TFile* file = new TFile(recPointsName, "READ");
+      {
+       TString str = "Event";
+       str+= nEv; 
+       file->cd(str.Data());
+      }
+      
+      TTree* tRecPoints = dynamic_cast<TTree*> (gDirectory->Get("TreeR"));
+      
+      AliTPCClustersRow *row = new AliTPCClustersRow();
+      
+      tRecPoints->GetBranch("Segment")->SetAddress(&row);
+
+      for (Int_t iRow = 0; iRow < tRecPoints->GetEntriesFast(); iRow++) {
+       tRecPoints->GetEntry(iRow);             
+       for (Int_t nCl = 0; nCl < row->GetArray()->GetEntriesFast(); nCl++){
+                 
+         AliTPCclusterMI * cl = static_cast<AliTPCclusterMI*> (row->GetArray()->UncheckedAt(nCl));
+         
+         int slice=-1, padRow=-1;
+                 
+         if( !AliHLTTPCTransform::Sector2Slice(slice, padRow, cl->GetDetector(), cl->GetRow() ) ) continue;      
+         int patch = AliHLTTPCTransform::GetPatch(padRow);
+         padRow -=AliHLTTPCTransform::GetFirstRow(patch);
+
+         double yW = AliHLTTPCTransform::GetPadPitchWidth(patch);
+         double zW = AliHLTTPCTransform::GetZWidth();
+
+         AliHLTTPCRawCluster rawCluster;
+         rawCluster.SetPadRow(padRow);
+         rawCluster.SetPad(cl->GetPad());
+         rawCluster.SetTime(cl->GetTimeBin());
+
+         rawCluster.SetSigmaY2(cl->GetSigmaY2()/yW/yW);
+         rawCluster.SetSigmaZ2(cl->GetSigmaZ2()/zW/zW);
+         rawCluster.SetCharge((UShort_t)cl->GetQ()); 
+         rawCluster.SetQMax((UShort_t)cl->GetMax());
+         
+         AliHLTTPCClusterMCLabel mc;
+         int l0 = ( cl->GetLabel(0)>=0 );
+         int l1 = ( cl->GetLabel(1)>=0 );
+         int l2 = ( cl->GetLabel(2)>=0 );
+
+         int nmc = l0+l1+l2;
+         float w = (nmc>0) ? cl->GetQ()/nmc :0;
+
+         mc.fClusterID[0].fMCID = cl->GetLabel(0);  mc.fClusterID[0].fWeight=l0*w;
+         mc.fClusterID[1].fMCID = cl->GetLabel(1);  mc.fClusterID[0].fWeight=l1*w;
+         mc.fClusterID[2].fMCID = cl->GetLabel(2);  mc.fClusterID[0].fWeight=l2*w;
+
+         nClusters++;
+
+         if( !merger.CheckCandidate(slice, patch, rawCluster) ) continue;
+         
+         
+         int id = merger.AddCandidate(slice, patch, ~AliHLTUInt32_t(0), rawCluster, mc);
+         if( id>=0 ){
+           nAccepted++;
+           vClusters.push_back(*cl);
+         }
+
+         //if( (id>=0) &&(merger.GetRecords()[id].GetBorder()>=0) ) nAccepted++;
+         
+       }
+      } // loop over clusters                    
+  
+      delete row;
+      file->Close();  
+    }
+
+
+    cout<<"Merge..."<<endl;
+
+    static int nMerged = 0;
+    nMerged+=merger.Merge();
+
+    cout<<"Read "<<nClusters<<" clusters, accepted "<<nAccepted<<" clusters "<<endl;
+    cout<<"Merged "<<2*nMerged<<" clusters "
+       <<" ("<<100.*2*nMerged/nAccepted<<" % of "<<nAccepted<< " accepted border clusters )"<<endl;
+
+    for( int iB=0; iB<merger.GetNSlices()*merger.GetNBorders(); iB+=2 ){
+       const AliHLTTPCHWClusterMerger::AliBorderRecord 
+        *b1 = merger.GetBorderClusters() + merger.GetBorderFirstCluster(iB),
+        *b2 = merger.GetBorderClusters() + merger.GetBorderFirstCluster(iB+1);
+       int 
+        n1 = merger.GetBorderNClusters(iB),
+        n2 = merger.GetBorderNClusters(iB+1);
+       if( n1<=0 || n2<=0 ) continue;
+
+       //cout<<iB<<" "<<n1<<" "<<n2<<endl;
+       //cout<<iB<<" "<<merger.GetBorderFirstCluster(iB)<<" "<<n1<<endl;
+       //cout<<iB+1<<" "<<merger.GetBorderFirstCluster(iB+1)<<" "<<n2<<endl;
+       /*
+       cout<<"Border "<<iB<<":"<<endl;
+       cout<<"  left: "<<n1<<" : ";
+       for( int i=0; i<n1; i++ ){
+        const AliTPCclusterMI &c = vClusters[b1[i].fClusterRecordID];
+        cout<<" ("<<c.GetLabel(0)<<" "<<(int)c.GetTimeBin()<<") ";
+       }
+       cout<<endl;
+       cout<<"  right: "<<n2<<" : ";
+       for( int i=0; i<n2; i++ ){
+        const AliTPCclusterMI &c = vClusters[b2[i].fClusterRecordID];
+        cout<<" ("<<c.GetLabel(0)<<" "<<(int)c.GetTimeBin()<<") ";
+       }
+       cout<<endl;
+       */
+
+       for( int i1=0; i1<n1; i1++ ){
+        const AliTPCclusterMI &c1 = vClusters[b1[i1].fClusterRecordID];
+        int lab = c1.GetLabel(0);
+        if( lab<0 ) continue;
+
+        // check if the label has been already checked
+        bool b=0;
+        for( int i=0; i<i1; i++ ){
+          const AliTPCclusterMI &c = vClusters[b1[i].fClusterRecordID];
+          if( c.GetLabel(0) == lab ){
+            b = 1;
+            break;
+          }
+        }
+        if( b ) continue;
+
+        // check if there is a partner with the same label
+        bool isPartner = 0;
+        bool isMerged = 0;
+        for( int i2=0; i2<n2; i2++ ){
+          const AliTPCclusterMI &c2 = vClusters[b2[i2].fClusterRecordID];
+          if( c2.GetLabel(0) != lab ) continue;
+          if( fabs(c2.GetTimeBin()-c1.GetTimeBin() ) > 5 ) continue;
+          isPartner = 1 ;
+          // check if it has been a merged with a cluster with the same label       
+          int id = merger.GetRecords()[b2[i2].fClusterRecordID].IsMergedTo();
+          if( id<0 ) continue;
+          const AliTPCclusterMI &c21 = vClusters[id];
+          if( c21.GetLabel(0) == lab ) isMerged = 1;    
+        }
+        if( isPartner ) statNGood++;
+        if( isMerged ) statNMerged++;
+       }
+       /*
+       cout<<" Merging efficiency: merged "<<statNMerged<<" of "<<statNGood<<" tracks : "
+          <<100.*statNMerged/( statNGood>0 ?statNGood :1)<<" %"<<endl;
+       */ 
+    }
+
+       cout<<" Merging efficiency: merged "<<statNMerged<<" of "<<statNGood<<" tracks : "
+          <<100.*statNMerged/( statNGood>0 ?statNGood :1)<<" %"<<endl;
+  }// ***** End of the loop over events
+  delete rl;
+  
+
+  /*
+  gStyle->SetOptStat(111110);
+  gStyle->SetOptFit(1);   
+  //hphidist->Draw();
+  Int_t minc=33; 
+  {
+    TCanvas *c1=new TCanvas("c0","Offline vs HLT clusters",0,0,1500,800);  
+    c1->Divide(1,3);
+    c1->cd(1);
+    hQ->Draw();
+    c1->cd(2);
+    hQ10->Draw();
+    c1->cd(3);
+    hQ30->Draw();
+    c1->cd(0);
+    c1->Update();
+  }
+  */  
+  return 0;
+}
+
+#endif