[u/mrichter/AliRoot.git] / ITS / AliITSandTPCrecoV2.C

////////////////////////////////////////////////////////////////////////
//
// 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;

}
Commit	Line	Data
13bc8df2	1	////////////////////////////////////////////////////////////////////////
	2	//
	3	// AliITSandTPCrecoV2.C
	4	//
	5	//
	6	////////////////////////////////////////////////////////////////////////
	7
	8
	9	/****************************************************************************
	10	* This macro is supposed to do reconstruction in the barrel ALICE trackers *
	11	* (Make sure you have TPC digits and ITS hits before using this macro !!!) *
	12	* It does the following steps *
	13	* 1) TPC cluster finding *
	14	* 2) TPC track finding *
	15	* 3) ITS cluster finding - fast recpoints - *
	16	* 4) ITS tracking *
	17	* TPC part taken from AliTPCTracking.C - Jiri Chudoba *
	18	* The TPC part will be removed *
	19	* M.Masera 30/09/2002 09:40 *
	20	* Last revision: 4/3/2003 *
	21	* *
	22	* Use: *
	23	* The main function is AliITSandTPCrecoV2 All the input arguments have a *
	24	* default value. *
	25	* nEvents: number of events to be processed. *
	26	* If nEvents<0, all the available events are used starting from *
	27	* the first one *
	28	* firstEvent: is the first event to be analyzed. *
	29	* It is set to 0 if nEvents<0 (default) *
	30	* fileNameHits: input root file with Kine and hits *
	31	* fileNameDigits: input root file with ITS and TPC digits and with ITS *
	32	* recpoints. Only fast points are used at the moment. *
	33	* An option to use slow points should be added *
	34	* fileNameClusters: output root file with TPC clusters *
	35	* fileNameTracks: output root file with TPC tracks *
	36	* fileNameITSClusters: output root file with ITS clusters *
	37	* fileNameITSTracks: output root file with ITS tracks *
	38	* *
	39	* AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0, *
	40	* TString fileNameHits="galice.root", *
	41	* TString fileNameDigits="galiceSDR.root", *
	42	* TString fileNameClusters="tpc.clusters.root", *
	43	* TString fileNameTracks="tpc.tracks.root", *
	44	* TString fileNameITSClusters="its.clusters.root", *
	45	* TString fileNameITSTracks="its.tracks.root"); *
	46	* *
	47	* Global variables:
	48	*
	49	* Int_t gDEBUG = 0; // NO verbose printouts if ==0 *
	50	* TArrayD gzvtx = 0; // Z coordinate of primary vertices
	51	* Int_t gCurrEve = 0; // current event number *
	52	* Bool_t gPPMode = kTRUE; // it must be set to kTRUE for pp interactions *
	53	* Bool_t gFastReco = kTRUE; // it must set to kTRUE for fast reconstr. *
	54	* const Double_t kgVertexNotFound = -100; // the Z coord. of the *
	55	* primary vertex is set to this value *
	56	* when the vertexing fails *
	57	* Bool_t gUseNewClustersV2 = kFALSE; // if kTRUE the AliITSclustereV2 class
	58	* is used to produce V2 clusters. *
	59	* ITS recpoints are used otherwise *
	60	****************************************************************************/
	61
	62	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	63	#include <TH1.h>
	64	#include <TH2.h>
65	#include <TProfile.h>
66	#include <TBranch.h>
67	#include <TParticle.h>
68	#include <TRandom.h>
69	#include "alles.h"
70	#include "AliMagF.h"
71	#include "AliMagFCM.h"
72	#include "AliTPCtracker.h"
73	#include "AliTPCclusterer.h"
74	#include "Riostream.h"
75	#include "TString.h"
76	#include "TArrayD.h"
77	#include "TLine.h"
78	#include "TTree.h"
79	#include "TNtuple.h"
80	#include "TParticle.h"
81	#include "TText.h"
82	#include "TSystem.h"
83	#include "TVector.h"
84	#include "TObjArray.h"
85	#include "TStyle.h"
86	#include "TCanvas.h"
87	#include "AliITS.h"
88	#include "AliITSclustererV2.h"
89	#include "AliITSgeom.h"
90	#include "AliITSRecPoint.h"
91	#include "AliITSclusterV2.h"
92	#include "AliITStrackerV2.h"
93	#include "AliITSVertexerPPZ.h"
94	#include "AliITSVertexerIons.h"
95	#include "AliRun.h"
96	#include "AliHeader.h"
97	#include "AliGenEventHeader.h"
98	#include "TArrayF.h"
99	#include "AliGenPythiaEventHeader.h"
100
101	#endif
102
103
104	Int_t gDEBUG = 0;
105	TArrayD *gzvtx = 0;
106	Int_t gCurrEve = 0;
107	Bool_t gPPMode = kTRUE;
108	Bool_t gFastReco = kTRUE;
109	Bool_t gUseNewClustersV2 = kFALSE;
110	const Double_t kgVertexNotFound = -100.;
111
112
113	Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
114	TString fileNameDigits="galiceSDR.root",
115	TString fileNameClusters="tpc.clusters.root");
116	Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
117	TFile* fileDigits, TFile* fileClusters,
118	AliTPCParam* paramTPC=0);
119	Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
120	TString fileNameClusters="tpc.clusters.root",
121	TString fileNameTracks="tpc.tracks.root");
122	Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
123	TFile* fileClusters, TFile* fileTracks,
124	AliTPCParam* paramTPC=0);
125
126	Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
127	TString fileNameHits="galice.root",
128	TString fileNameDigits="galiceSDR.root",
129	TString fileNameClusters="tpc.clusters.root",
130	TString fileNameTracks="tpc.tracks.root");
131	Int_t AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,
132	TString fileNameHits="galice.root",
133	TString fileNameDigits="galiceSDR.root",
134	TString fileNameClusters="tpc.clusters.root",
135	TString fileNameTracks="tpc.tracks.root",
136	TString fileNameITSClusters="its.clusters.root",
137	TString fileNameITSTracks="its.tracks.root");
138	Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first);
139	Int_t ITSFindTracks(TString fileNameTracks, TString fileNameITSClusters, TString fileNameITSTracks, Int_t n, Int_t first, Int_t vc, Int_t vc2);
140	void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints);
141	Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile in,TFile out);
142	////////////////////////////////////////////////////////////////////////
143	Int_t AliITSandTPCrecoV2(Int_t nEvents, Int_t firstEvent,
144	TString fileNameHits,
145	TString fileNameDigits,
146	TString fileNameClusters,
147	TString fileNameTracks,
148	TString fileNameITSClusters,
149	TString fileNameITSTracks) {
150	const Char_t *name="AliITSandTPCrecoV2";
151
152	gBenchmark->Start(name);
153
154	// Set the conversion constant for the Kalman filter
155	// and set the gAlice pointer
156	Char_t finame = (Char_t )fileNameHits.Data();
157	AliTracker::SetFieldFactor(finame,kFALSE);
158
159	if(nEvents < 0 ) {
160	nEvents = (Int_t)gAlice->TreeE()->GetEntries();
161	firstEvent = 0;
162	cout << "Processing events from " << firstEvent << " up to " << nEvents-1 <<endl;
163	}
164
165	if(fileNameDigits!=fileNameHits)gAlice->SetTreeDFileName(fileNameDigits.Data());
166	gzvtx = new TArrayD(nEvents);
167
168	// measure Z vertex coordinate for seeding
169
170	FindZV(nEvents,firstEvent,fileNameHits,fileNameDigits);
171
172	// TPC tracking - abort macro execution if tracking fails
173	if(AliTPCTracking(nEvents,firstEvent,fileNameHits,fileNameDigits,
174	fileNameClusters,fileNameTracks)) {
175	cerr << "TPC tracking failed \n";
176	return 1;
177	}
178
179	// ITS clustering
180	if(ITSFindClusters(fileNameHits,fileNameDigits,fileNameITSClusters,nEvents,firstEvent)) {
181	cerr << "ITS clustering failed \n";
182	return 2;
183	}
184
185	Int_t firstpass = 1; // 1=vert. constraint; 0=no vert. constr.; -1=skip pass
186	Int_t secondpass = 0; // as above
187	if(ITSFindTracks(fileNameTracks,fileNameITSClusters,fileNameITSTracks,nEvents,firstEvent,firstpass,secondpass)) {
188	cerr << "ITS tracking failed \n";
189	return 3;
190	}
191
192
193	gBenchmark->Stop(name);
194	gBenchmark->Show(name);
195	return 0;
196	}
197	////////////////////////////////////////////////////////////////////////
198	Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
199	TString fileNameHits,
200	TString fileNameDigits,
201	TString fileNameClusters,
202	TString fileNameTracks) {
203
204
205	// ******** Find TPC clusters ********* //
206	if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
207	cerr<<"Failed to get TPC clusters: !\n";
208	return 1;
209	}
210
211	// ******** Find TPC tracks ********* //
212	// if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
213	if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
214	cerr<<"Failed to get TPC tracks !\n";
215	return 2;
216	}
217
218	return 0;
219	}
220
221	////////////////////////////////////////////////////////////////////////
222	Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
223	TString fileNameDigits, TString fileNameClusters) {
224
225	Int_t rc;
226	const Char_t *name="TPCFindClusters";
227	if (gDEBUG>1) cout<<name<<" starts...\n";
228	if (gDEBUG>1) gBenchmark->Start(name);
229	TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
230	TFile *fileDigits = TFile::Open(fileNameDigits);
231	if (!fileDigits->IsOpen()) {
232	cerr<<"Cannnot open "<<fileNameDigits<<" !\n";
233	return 1;
234	}
235	if (!fileClusters->IsOpen()) {
236	cerr<<"Cannnot open "<<fileNameClusters<<" !\n";
237	return 1;
238	}
239
240	rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
241
242	fileDigits->Close();
243	fileClusters->Close();
244	delete fileDigits;
245	delete fileClusters;
246	if (gDEBUG>1) gBenchmark->Stop(name);
247	if (gDEBUG>1) gBenchmark->Show(name);
248
249	return rc;
250	}
251	////////////////////////////////////////////////////////////////////////
252	Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
253	TFile* fileDigits, TFile* fileClusters,
254	AliTPCParam* paramTPC) {
255
256	fileDigits->cd();
257	if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits);
258	if (!paramTPC) return 1;
259
260	for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
261	if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
262	AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
263	}
264	return 0;
265	}
266	////////////////////////////////////////////////////////////////////////
267	Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
268	TString fileNameClusters, TString fileNameTracks) {
269
270	Int_t rc = 0;
271	const Char_t *name="TPCFindTracks";
272	if (gDEBUG>1) cout<<name<<" starts"<<endl;
273	if (gDEBUG>1) gBenchmark->Start(name);
274	TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
275	TFile *fileClusters =TFile::Open(fileNameClusters);
276
277	rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
278
279	fileClusters->Close();
280	fileTracks->Close();
281	delete fileClusters;
282	delete fileTracks;
283	if (gDEBUG>1) gBenchmark->Stop(name);
284	if (gDEBUG>1) gBenchmark->Show(name);
285	return rc;
286
287	}
288	////////////////////////////////////////////////////////////////////////
289	Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
290	TFile fileClusters, TFile fileTracks,
291	AliTPCParam* paramTPC) {
292
293	Int_t rc = 0;
294	if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters);
295	if (!paramTPC) return 1;
296	Double_t vertex[3];
297	for(Int_t j=0; j<3; j++){vertex[j]=0.;}
298	AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
299	for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
300	if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
301	tracker->SetEventNumber(iEvent);
302	if((*gzvtx)[iEvent] > -100){
303	vertex[2] = (*gzvtx)[iEvent];
304	}
305	else {
306	vertex[2]=0;
307	}
308	tracker->SetVertex(vertex);
309	fileClusters->cd();
310	rc = tracker->Clusters2Tracks(0,fileTracks);
311	}
312	delete tracker;
313	return rc;
314	}
315
316
317	////////////////////////////////////////////////////////////////////////
318	Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first) {
319	Int_t rc=0;
320	const Char_t *name="ITSFindClusters";
321	cerr<<'\n'<<name<<"...\n";
322	gBenchmark->Start(name);
323	TFile *out=TFile::Open(fileNameITSClusters,"recreate");
324	TFile in = (TFile)gROOT->GetListOfFiles()->FindObject(inname.Data());
325	if(gUseNewClustersV2){
326	cout<<"Direct method\n";
327	return DirectClusterFinder(n,first,in,out);
328	}
329	else {
330	cout<<"Clusters through rec points \n";
331	}
332
333	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
334	if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
335	AliITSgeom *geom=ITS->GetITSgeom();
336	AliITSclustererV2 *clusterer = new AliITSclustererV2(geom);
337	out->cd();
338	geom->Write();
339
340	Int_t ev;
341	if(gFastReco && gDEBUG>0)cout <<"Using fast points\n";
342	else cout<<"Using slow points\n";
343	for (ev = first; ev<n; ev++){
344	cout<<"ITSFindClusters: processing event "<<ev<<endl;
345	in->cd();
346	gAlice->GetEvent(ev);
347
348	TTree *pTree=gAlice->TreeR();
349	if (!pTree){
350	cerr << "ITSFindClusters: can't get TreeR\n";
351	return 1;
352	}
353	TBranch *branch = 0;
354	if (gFastReco) {
355	branch = pTree->GetBranch("ITSRecPointsF");
356	}
357	else {
358	branch = pTree->GetBranch("ITSRecPoints");
359	}
360	if (!branch) {
361	cerr << "ITSFindClusters: can't get ITSRecPoints branch \n";
362	return 2;
363	}
364
365	out->cd();
366	TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
367	char cname[100];
368	sprintf(cname,"TreeC_ITS_%d",ev);
369
370	TTree *cTree=new TTree(cname,"ITS clusters");
371	cTree->Branch("Clusters",&clusters);
372
373	TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
374	branch->SetAddress(&points);
375
376	TClonesArray &cl=*clusters;
377	Int_t nclusters=0;
378	Int_t nentr=(Int_t)pTree->GetEntries();
379	AliITSgeom *geom=ITS->GetITSgeom();
380	Float_t lp[5];
381	Double_t rot[9];
382	Int_t lab[6];
383	for (Int_t i=0; i<nentr; i++) {
384	branch->GetEvent(i);
385	clusterer->RecPoints2Clusters(points,i,clusters);
386	cTree->Fill(); clusters->Delete();
387	points->Clear();
388	}
389	cTree->Write();
390	delete cTree; delete clusters; points->Delete(); delete points;
391	}
392	delete clusterer;
393
394	out->Close();
395
396	return rc;
397	}
398	////////////////////////////////////////////////////////////////////////
399	Int_t ITSFindTracks(TString inname, TString inname2, TString outname, Int_t n, Int_t first, Int_t vc, Int_t vc2) {
400	Int_t rc=0;
401	if(gPPMode){
402	AliITStrackV2::SetSigmaYV(0.015);
403	AliITStrackV2::SetSigmaZV(0.017);
404	}
405	TFile *out=TFile::Open(outname.Data(),"recreate");
406	if (!out->IsOpen()) {
407	cerr<<"Can't open file "<<outname.Data()<<endl;
408	return 1;
409	}
410	TFile *in =TFile::Open(inname.Data());
411	if (!in->IsOpen()) {
412	cerr<<"Can't open file "<<inname.Data()<<endl;
413	return 1;
414	}
415	TFile *in2 =TFile::Open(inname2.Data());
416	if (!in2->IsOpen()) {
417	cerr<<"Can't open file "<<inname2.Data()<<endl;
418	return 1;
419	}
420
421	AliITSgeom geom=(AliITSgeom)gFile->Get("AliITSgeom");
422	if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
423	Double_t vrtc[3];
424	AliITStrackerV2 tracker(geom);
425	for (Int_t i=first;i<n;i++){
426	tracker.SetEventNumber(i);
427	Double_t vtx=(*gzvtx)[i];
428	cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
429	cout<<"ITSFindTracks -- event "<<i<<endl;
430	cout<<"Computed vertex position: "<<vtx<<endl;
431	in2->cd();
432	vrtc[0]=0.;
433	vrtc[1]=0.;
434	vrtc[2]=vtx;
435	if(vtx != kgVertexNotFound){
436	tracker.SetVertex(vrtc); // set vertex only if it was computed
437	}
438	else {
439	vc = 0; // use vertex contraint only if vertex info is available
440	vc2 = -1;
441	}
442
443	tracker.SetupFirstPass(&vc);
444
445	tracker.SetupSecondPass(&vc2);
446
447	rc=tracker.Clusters2Tracks(in,out);
448	}
449
450	in->Close();
451	in2->Close();
452	out->Close();
453
454	return rc;
455	}
456
457	//////////////////////////////////////////////////////////////////
458	void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints){
459	TFile in = (TFile)gROOT->GetListOfFiles()->FindObject(FileNameHits.Data());
460	TFile *fo = new TFile("AliITSVertices.root","recreate");
461	if(FileNameHits!=FileWithRecPoints)gAlice->SetTreeRFileName(FileWithRecPoints);
462	AliITSVertexer *findVert;
463	if(gPPMode){
464	findVert = new AliITSVertexerPPZ(in,fo);
465	}
466	else {
467	findVert = new AliITSVertexerIons(in,fo);
468	}
469	Int_t last = first + nEvents - 1;
470	AliITSVertex *vert = 0;
471	for(Int_t i=first; i<=last; i++){
472	gAlice->GetEvent(i);
473	// The true Z coord. is fetched for comparison
474	AliHeader *header = gAlice->GetHeader();
475	AliGenEventHeader* genEventHeader = header->GenEventHeader();
476	TArrayF primaryVertex(3);
477	genEventHeader->PrimaryVertex(primaryVertex);
478	vert = findVert->FindVertexForCurrentEvent(i);
479	if(vert){
480	Double_t pos[3];
481	for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
482	vert->SetTruePos(pos);
483	Double_t found = vert->GetZv();
484	gzvtx->AddAt(found,i);
485	findVert->WriteCurrentVertex();
486	}
487	else {
488	gzvtx->AddAt(kgVertexNotFound,i);
489	}
490	}
491	delete findVert;
492	fo->Close();
493	delete fo;
494	}
495
496
497
498	Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile in,TFile out){
499
500
501	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
502	if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
503
504	AliITSgeom *geom=ITS->GetITSgeom();
505
506	geom->Write();
507
508	TStopwatch timer;
509	AliITSclustererV2 clusterer(geom);
510	for (Int_t i=first; i<eventn; i++) {
511	cerr<<"Processing event number: "<<i<<endl;
512	gAlice->GetEvent(i);
513	//ITS->MakeTreeC(); //To make the V1 cluster finders happy
514	clusterer.SetEvent(i);
515	if (gFastReco) clusterer.Digits2Clusters(in,out);
516	else clusterer.Hits2Clusters(in,out);
517	}
518	timer.Stop(); timer.Print();
519
520	out->Close();
521
522	return 0;
523
524	}
525