macro to reconstruct N events with ITS V2 tracking. TPC is also reconstructed (provis...
authormasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Mar 2003 18:32:56 +0000 (18:32 +0000)
committermasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 Mar 2003 18:32:56 +0000 (18:32 +0000)
ITS/AliITSandTPCrecoV2.C [new file with mode: 0644]

diff --git a/ITS/AliITSandTPCrecoV2.C b/ITS/AliITSandTPCrecoV2.C
new file mode 100644 (file)
index 0000000..40d1f34
--- /dev/null
@@ -0,0 +1,525 @@
+////////////////////////////////////////////////////////////////////////
+//
+// AliITSandTPCrecoV2.C 
+//
+//
+////////////////////////////////////////////////////////////////////////
+
+
+/****************************************************************************
+ * This macro is supposed to do reconstruction in the barrel ALICE trackers *
+ * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
+ * It does the following steps                                              *
+ *                   1) TPC cluster finding                                 *
+ *                   2) TPC track finding                                   *
+ *                   3) ITS cluster finding - fast recpoints -              *
+ *                   4) ITS tracking                                        *
+ * TPC part taken from AliTPCTracking.C - Jiri Chudoba                      *
+ * The TPC part will be removed                                             *
+ * M.Masera 30/09/2002 09:40                                                *
+ * Last revision: 4/3/2003                                                  *
+ *                                                                          *
+ * Use:                                                                     *
+ * The main function is AliITSandTPCrecoV2 All the input arguments have a   *
+ * default value.                                                           *
+ * nEvents: number of events to be processed.                               *
+ *          If nEvents<0, all the available events are used starting from   *
+ *          the first one                                                   *
+ * firstEvent: is the first event to be analyzed.                           *
+ *             It is set to 0 if nEvents<0 (default)                        *
+ * fileNameHits: input root file with Kine and hits                         *
+ * fileNameDigits: input root file with ITS and TPC digits and with ITS     *
+ *                 recpoints. Only fast points are used at the moment.      *
+ *                 An option to use slow points should be added             *
+ * fileNameClusters: output root file with TPC clusters                     *
+ * fileNameTracks: output root file with TPC tracks                         *
+ * fileNameITSClusters: output root file with ITS clusters                  *
+ * fileNameITSTracks: output root file with ITS tracks                      *
+ *                                                                          *
+ *   AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,               * 
+ *                  TString fileNameHits="galice.root",                    *
+ *                  TString fileNameDigits="galiceSDR.root",               *
+ *                  TString fileNameClusters="tpc.clusters.root",          *
+ *                  TString fileNameTracks="tpc.tracks.root",              *
+ *                   TString fileNameITSClusters="its.clusters.root",       *
+ *                   TString fileNameITSTracks="its.tracks.root");          *
+ *                                                                          *
+ *  Global variables:
+ *  
+ * Int_t gDEBUG = 0;   // NO verbose printouts if ==0                       *
+ * TArrayD *gzvtx = 0; // Z coordinate of primary vertices                  *
+ * Int_t gCurrEve = 0; // current event number                              *
+ * Bool_t gPPMode = kTRUE;  // it must be set to kTRUE for pp interactions  *
+ * Bool_t gFastReco = kTRUE;  // it must set to kTRUE for fast reconstr.    *
+ * const Double_t kgVertexNotFound = -100;  // the Z coord. of the          *
+ *                                 primary vertex is set to this value      *
+ *                                 when the vertexing fails                 *
+ * Bool_t gUseNewClustersV2 = kFALSE;  // if kTRUE the AliITSclustereV2 class
+ *                                 is used to produce V2 clusters.          *
+ *                                 ITS recpoints are used otherwise         *
+ ****************************************************************************/
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TBranch.h>
+#include <TParticle.h>
+#include <TRandom.h>
+#include "alles.h"
+#include "AliMagF.h"
+#include "AliMagFCM.h"
+#include "AliTPCtracker.h"
+#include "AliTPCclusterer.h"
+#include "Riostream.h"
+#include "TString.h"
+#include "TArrayD.h"
+#include "TLine.h"
+#include "TTree.h"
+#include "TNtuple.h"
+#include "TParticle.h"
+#include "TText.h"
+#include "TSystem.h"
+#include "TVector.h"
+#include "TObjArray.h"
+#include "TStyle.h"
+#include "TCanvas.h"
+#include "AliITS.h"
+#include "AliITSclustererV2.h"
+#include "AliITSgeom.h"
+#include "AliITSRecPoint.h"
+#include "AliITSclusterV2.h"
+#include "AliITStrackerV2.h"
+#include "AliITSVertexerPPZ.h"
+#include "AliITSVertexerIons.h"
+#include "AliRun.h"
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "TArrayF.h"
+#include "AliGenPythiaEventHeader.h"
+
+#endif
+
+
+Int_t gDEBUG = 0;
+TArrayD *gzvtx = 0;
+Int_t gCurrEve = 0;
+Bool_t gPPMode = kTRUE;
+Bool_t gFastReco = kTRUE;
+Bool_t gUseNewClustersV2 = kFALSE; 
+const Double_t kgVertexNotFound = -100.;
+
+
+Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
+                     TString fileNameDigits="galiceSDR.root", 
+                     TString fileNameClusters="tpc.clusters.root");
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TFile* fileDigits, TFile* fileClusters, 
+                     AliTPCParam* paramTPC=0);
+Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
+                   TString fileNameClusters="tpc.clusters.root",
+                   TString fileNameTracks="tpc.tracks.root");
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TFile* fileClusters, TFile* fileTracks,
+                   AliTPCParam* paramTPC=0);
+
+Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
+                    TString fileNameHits="galice.root",
+                    TString fileNameDigits="galiceSDR.root",
+                    TString fileNameClusters="tpc.clusters.root",
+                    TString fileNameTracks="tpc.tracks.root");
+Int_t AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,
+                       TString fileNameHits="galice.root",
+                       TString fileNameDigits="galiceSDR.root",
+                       TString fileNameClusters="tpc.clusters.root",
+                       TString fileNameTracks="tpc.tracks.root",
+                       TString fileNameITSClusters="its.clusters.root",
+                       TString fileNameITSTracks="its.tracks.root");
+Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first);
+Int_t ITSFindTracks(TString fileNameTracks, TString fileNameITSClusters, TString fileNameITSTracks, Int_t n, Int_t first, Int_t vc, Int_t vc2);
+void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints);
+Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile *out);
+////////////////////////////////////////////////////////////////////////
+Int_t AliITSandTPCrecoV2(Int_t nEvents, Int_t firstEvent,
+                       TString fileNameHits,
+                       TString fileNameDigits,
+                       TString fileNameClusters,
+                       TString fileNameTracks,
+                       TString fileNameITSClusters,
+                       TString fileNameITSTracks) {
+  const Char_t *name="AliITSandTPCrecoV2";
+
+  gBenchmark->Start(name);
+
+  // Set the conversion constant for the Kalman filter
+  // and set the gAlice pointer 
+  Char_t *finame = (Char_t *)fileNameHits.Data();
+  AliTracker::SetFieldFactor(finame,kFALSE);
+
+  if(nEvents < 0 ) {
+    nEvents = (Int_t)gAlice->TreeE()->GetEntries();
+    firstEvent = 0;
+    cout << "Processing events from " << firstEvent << " up to " << nEvents-1 <<endl;
+  }
+
+  if(fileNameDigits!=fileNameHits)gAlice->SetTreeDFileName(fileNameDigits.Data());
+  gzvtx = new TArrayD(nEvents);
+
+  // measure Z vertex coordinate for seeding
+  FindZV(nEvents,firstEvent,fileNameHits,fileNameDigits); 
+  
+  // TPC tracking - abort macro execution if tracking fails
+  if(AliTPCTracking(nEvents,firstEvent,fileNameHits,fileNameDigits,
+                   fileNameClusters,fileNameTracks)) {
+    cerr << "TPC tracking failed \n";
+    return 1;
+  }
+  // ITS clustering
+  if(ITSFindClusters(fileNameHits,fileNameDigits,fileNameITSClusters,nEvents,firstEvent)) {
+    cerr << "ITS clustering failed \n";
+    return 2;
+  }
+
+  Int_t firstpass = 1;  // 1=vert. constraint; 0=no vert. constr.; -1=skip pass
+  Int_t secondpass = 0; // as above
+  if(ITSFindTracks(fileNameTracks,fileNameITSClusters,fileNameITSTracks,nEvents,firstEvent,firstpass,secondpass)) {
+    cerr << "ITS tracking failed \n";
+    return 3;
+  }
+
+
+  gBenchmark->Stop(name);
+  gBenchmark->Show(name);
+  return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
+                     TString fileNameHits,
+                     TString fileNameDigits,
+                     TString fileNameClusters,
+                     TString fileNameTracks) {
+
+
+  // ********** Find TPC clusters *********** //
+  if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
+    cerr<<"Failed to get TPC clusters: !\n";
+    return 1;
+  }      
+
+  // ********** Find TPC tracks *********** //
+  //   if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
+  if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
+    cerr<<"Failed to get TPC tracks !\n";
+    return 2;
+  }
+
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TString fileNameDigits, TString fileNameClusters) {
+  
+  Int_t rc;
+  const Char_t *name="TPCFindClusters";
+  if (gDEBUG>1) cout<<name<<" starts...\n";
+  if (gDEBUG>1) gBenchmark->Start(name);
+  TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
+  TFile *fileDigits = TFile::Open(fileNameDigits);
+  if (!fileDigits->IsOpen()) {
+    cerr<<"Cannnot open "<<fileNameDigits<<" !\n"; 
+    return 1;
+  }
+  if (!fileClusters->IsOpen()) {
+    cerr<<"Cannnot open "<<fileNameClusters<<" !\n"; 
+    return 1;
+  }
+
+  rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
+
+  fileDigits->Close();
+  fileClusters->Close();
+  delete fileDigits;
+  delete fileClusters;
+  if (gDEBUG>1) gBenchmark->Stop(name);
+  if (gDEBUG>1) gBenchmark->Show(name);
+
+  return rc;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
+                     TFile* fileDigits, TFile* fileClusters, 
+                     AliTPCParam* paramTPC) {
+
+  fileDigits->cd();
+  if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits);
+  if (!paramTPC) return 1;
+
+  for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
+    if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
+    AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
+  }
+  return 0;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TString fileNameClusters, TString fileNameTracks) {
+
+  Int_t rc = 0;
+  const Char_t *name="TPCFindTracks";
+  if (gDEBUG>1) cout<<name<<" starts"<<endl;
+  if (gDEBUG>1) gBenchmark->Start(name);
+  TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
+  TFile *fileClusters =TFile::Open(fileNameClusters);
+
+  rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
+
+  fileClusters->Close();
+  fileTracks->Close();
+  delete fileClusters;
+  delete fileTracks;
+  if (gDEBUG>1) gBenchmark->Stop(name);
+  if (gDEBUG>1) gBenchmark->Show(name);
+  return rc;
+
+}
+////////////////////////////////////////////////////////////////////////
+Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
+                   TFile *fileClusters, TFile * fileTracks,
+                   AliTPCParam* paramTPC) {
+
+  Int_t rc = 0;
+  if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters);
+  if (!paramTPC) return 1;
+  Double_t vertex[3];
+  for(Int_t j=0; j<3; j++){vertex[j]=0.;}
+  AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
+  for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
+    if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
+    tracker->SetEventNumber(iEvent);
+    if((*gzvtx)[iEvent] > -100){
+      vertex[2] = (*gzvtx)[iEvent];
+    }
+    else {
+      vertex[2]=0;
+    }
+    tracker->SetVertex(vertex);
+    fileClusters->cd();
+    rc = tracker->Clusters2Tracks(0,fileTracks);
+  }
+  delete tracker;
+  return rc;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first) {
+  Int_t rc=0;
+  const Char_t *name="ITSFindClusters";
+  cerr<<'\n'<<name<<"...\n";
+  gBenchmark->Start(name);
+  TFile *out=TFile::Open(fileNameITSClusters,"recreate");
+  TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(inname.Data());
+  if(gUseNewClustersV2){
+    cout<<"Direct method\n";
+    return DirectClusterFinder(n,first,in,out);
+  }
+  else {
+    cout<<"Clusters through rec points \n";
+  }
+  AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
+  if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
+  AliITSgeom *geom=ITS->GetITSgeom();
+  AliITSclustererV2 *clusterer = new AliITSclustererV2(geom);
+  out->cd();   
+  geom->Write();
+     
+  Int_t ev;
+  if(gFastReco && gDEBUG>0)cout <<"Using fast points\n";
+  else cout<<"Using slow points\n";
+  for (ev = first; ev<n; ev++){
+    cout<<"ITSFindClusters: processing event "<<ev<<endl;
+    in->cd();  
+    gAlice->GetEvent(ev);
+     
+    TTree *pTree=gAlice->TreeR();
+    if (!pTree){
+      cerr << "ITSFindClusters: can't get TreeR\n";
+      return 1;
+    }
+    TBranch *branch = 0;
+    if (gFastReco) {
+      branch = pTree->GetBranch("ITSRecPointsF");
+    }
+    else {
+      branch = pTree->GetBranch("ITSRecPoints");
+    }
+    if (!branch) {
+      cerr << "ITSFindClusters: can't get ITSRecPoints branch \n";
+      return 2;
+    }
+
+    out->cd();
+    TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+    char   cname[100];
+    sprintf(cname,"TreeC_ITS_%d",ev);
+  
+    TTree *cTree=new TTree(cname,"ITS clusters");
+    cTree->Branch("Clusters",&clusters);
+    TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
+    branch->SetAddress(&points);
+     
+    TClonesArray &cl=*clusters;
+    Int_t nclusters=0;
+    Int_t nentr=(Int_t)pTree->GetEntries();
+    AliITSgeom *geom=ITS->GetITSgeom();
+    Float_t lp[5];
+    Double_t rot[9];
+    Int_t lab[6];
+    for (Int_t i=0; i<nentr; i++) {
+      branch->GetEvent(i);
+      clusterer->RecPoints2Clusters(points,i,clusters);
+      cTree->Fill(); clusters->Delete();
+      points->Clear();
+    }
+    cTree->Write();
+    delete cTree; delete clusters; points->Delete(); delete points;
+  }
+  delete clusterer;
+
+  out->Close();
+
+  return rc;
+}
+////////////////////////////////////////////////////////////////////////
+Int_t ITSFindTracks(TString inname, TString inname2, TString outname, Int_t n, Int_t first, Int_t vc, Int_t vc2) {
+  Int_t rc=0;
+  if(gPPMode){
+    AliITStrackV2::SetSigmaYV(0.015);
+    AliITStrackV2::SetSigmaZV(0.017);
+  }
+  TFile *out=TFile::Open(outname.Data(),"recreate");
+  if (!out->IsOpen()) {
+    cerr<<"Can't open file "<<outname.Data()<<endl; 
+    return 1;
+  }
+  TFile *in =TFile::Open(inname.Data());
+  if (!in->IsOpen()) {
+    cerr<<"Can't open file "<<inname.Data()<<endl; 
+    return 1;
+  }
+  TFile *in2 =TFile::Open(inname2.Data());
+  if (!in2->IsOpen()) {
+    cerr<<"Can't open file "<<inname2.Data()<<endl; 
+    return 1;
+  }
+
+  AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+  if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+  Double_t vrtc[3];
+  AliITStrackerV2 tracker(geom);
+  for (Int_t i=first;i<n;i++){
+    tracker.SetEventNumber(i);
+    Double_t vtx=(*gzvtx)[i];
+    cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
+    cout<<"ITSFindTracks -- event "<<i<<endl;
+    cout<<"Computed vertex position: "<<vtx<<endl;
+    in2->cd();
+    vrtc[0]=0.;
+    vrtc[1]=0.;
+    vrtc[2]=vtx;
+    if(vtx != kgVertexNotFound){
+      tracker.SetVertex(vrtc);   // set vertex only if it was computed
+    }
+    else {
+      vc = 0;  // use vertex contraint only if vertex info is available
+      vc2 = -1;
+    }
+
+    tracker.SetupFirstPass(&vc);
+
+    tracker.SetupSecondPass(&vc2);
+   
+    rc=tracker.Clusters2Tracks(in,out);
+  }
+
+  in->Close();
+  in2->Close();
+  out->Close();
+
+  return rc;
+}
+
+//////////////////////////////////////////////////////////////////
+void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints){
+  TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(FileNameHits.Data());
+  TFile *fo = new TFile("AliITSVertices.root","recreate"); 
+  if(FileNameHits!=FileWithRecPoints)gAlice->SetTreeRFileName(FileWithRecPoints);
+  AliITSVertexer *findVert;
+  if(gPPMode){
+    findVert = new AliITSVertexerPPZ(in,fo);
+  }
+  else {
+    findVert = new AliITSVertexerIons(in,fo);
+  }
+  Int_t last = first + nEvents - 1;
+  AliITSVertex *vert = 0;
+  for(Int_t i=first; i<=last; i++){
+    gAlice->GetEvent(i);
+    // The true Z coord. is fetched for comparison
+    AliHeader *header = gAlice->GetHeader();
+    AliGenEventHeader* genEventHeader = header->GenEventHeader();
+    TArrayF primaryVertex(3);
+    genEventHeader->PrimaryVertex(primaryVertex);
+    vert = findVert->FindVertexForCurrentEvent(i);
+    if(vert){
+      Double_t pos[3];
+      for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
+      vert->SetTruePos(pos);
+      Double_t found = vert->GetZv();
+      gzvtx->AddAt(found,i);
+      findVert->WriteCurrentVertex();
+    }
+    else {
+      gzvtx->AddAt(kgVertexNotFound,i);
+    }
+  }
+  delete findVert;
+  fo->Close();
+  delete fo;
+}
+
+
+
+Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile* out){
+
+
+   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
+   if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
+
+   AliITSgeom *geom=ITS->GetITSgeom();
+
+   geom->Write();
+
+   TStopwatch timer;
+   AliITSclustererV2 clusterer(geom);
+   for (Int_t i=first; i<eventn; i++) {
+       cerr<<"Processing event number: "<<i<<endl;
+       gAlice->GetEvent(i);
+       //ITS->MakeTreeC(); //To make the V1 cluster finders happy
+       clusterer.SetEvent(i);
+       if (gFastReco) clusterer.Digits2Clusters(in,out);
+       else                 clusterer.Hits2Clusters(in,out);
+   }
+   timer.Stop(); timer.Print();
+
+   out->Close();
+
+   return 0;
+
+}
+