/****************************************************************************
- * This macro is supposed to do reconstruction in the ITS via Kalman *
- * tracker V2. The ITStracker is feeded with parametrized TPC tracks *
+ * This macro performs track and vertex reconstruction in TPC and ITS. *
+ * The ITS Kalman tracker V2 is feeded "with" parameterized TPC tracks. *
* *
- * It does the following steps: *
+ * Reconstruction is performed in the following steps: *
* 1) TPC tracking parameterization *
- * 2) Fast points in ITS *
- * 3) ITS cluster finding V2 *
- * 4) Determine z position of primary vertex *
- * - read from event header in PbPb events *
- * - determined using points on SPD in pp events (M.Masera) *
- * 5) ITS track finding V2 *
- * * in pp, redetermine the z position of the primary vertex *
+ * 2) ITS clusters: slow or fast *
+ * 3) Primary vertex reconstruction *
+ * - read from event header for Pb-Pb events *
+ * - determined using points in pixels for pp/pA events *
+ * 4) ITS track finding V2 *
+ * - in pp/pA, redetermine the position of primary vertex *
* using the reconstructed tracks *
- * 6) Create a reference file with simulation info (p,PDG...) *
+ * 5) Create a reference file with simulation info (p,PDG...) *
* *
- * If TString mode="A" all 6 steps are executed *
- * If TString mode="B" only steps 3-6 are executed *
+ * If mode='A' all 5 steps are executed *
+ * If mode='B' only steps 4-5 are executed *
* *
- * (Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it *
- * from AliTPCtest.C & AliITStestV2.C by I.Belikov *
+ * Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it *
+ * (from AliTPCtest.C & AliITStestV2.C by I.Belikov) *
****************************************************************************/
-#if !defined(__CINT__) || defined(__MAKECINT__)
-//-- --- standard headers-------------
-#include <iostream.h>
-//--------Root headers ---------------
-#include <TFile.h>
-#include <TStopwatch.h>
-#include <TObject.h>
-#include <TSystem.h>
-#include <TH1.h>
-#include <TH2.h>
-#include <TProfile.h>
-#include <TTree.h>
-#include <TBranch.h>
-#include <TParticle.h>
-#include <TRandom.h>
-//----- AliRoot headers ---------------
-#include "alles.h"
-#include "AliRun.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
-#include "AliMagF.h"
-#include "AliModule.h"
-#include "AliArrayI.h"
-#include "AliDigits.h"
-#include "AliITS.h"
-#include "AliTPC.h"
-#include "AliITSgeom.h"
-#include "AliITSRecPoint.h"
-#include "AliITSclusterV2.h"
-#include "AliITSsimulationFastPoints.h"
-#include "AliITStrackerV2.h"
-#include "AliKalmanTrack.h"
-#include "AliTPCtrackerParam.h"
-//-------------------------------------
-#endif
// structure for track references
typedef struct {
//===== Functions definition =================================================
-Int_t TPCParamTracks(const Char_t *galiceName,const Char_t *outName,
- const Int_t coll,const Double_t Bfield,Int_t n);
+void CopyVtx(const Char_t *inName,const Char_t *outName);
-Int_t ITSHits2FastRecPoints(const Char_t *galiceName,Int_t n);
+void ITSFindClustersV2(Char_t SlowOrFast);
-Int_t ITSFindClustersV2(const Char_t *galiceName,const Char_t *outName,
- Int_t *skipEvt,Int_t n);
+void ITSFindTracksV2(Int_t *skipEvt);
-Int_t ZvtxFromHeader(const Char_t *galiceName,Int_t n);
+void ITSMakeRefFile(Int_t *skipEvt);
-Int_t ZvtxFromSPD(const Char_t *galiceName,Int_t n);
+void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
-Int_t ZvtxFromTracks(const Char_t *trkame,Int_t n);
+void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
-Int_t ZvtxFastpp(const Char_t *galiceName,Int_t n);
+void TPCParamTracks(Int_t coll,Double_t Bfield);
-void EvalZ(TH1F *hist,Int_t sepa,Float_t &av,Float_t &sig,Int_t ncoinc,
- TArrayF *zval,ofstream *deb);
+Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
-Bool_t VtxTrack(Double_t pt,Double_t d0rphi);
+void VtxFromHeader(const Char_t *outName,Bool_t smear);
-Double_t d0zRes(Double_t pt);
+void VtxFromTracks(const Char_t *outName);
-Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName,Int_t n);
+void ZvtxFromSPD(const Char_t *outName);
-Int_t MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
-
-Int_t ITSFindTracksV2(const Char_t *galiceName,const Char_t *inName,
- const Char_t *inName2,const Char_t *outName,
- Int_t *skipEvt,TString vtxMode,Int_t n);
-
-Int_t ITSMakeRefFile(const Char_t *galiceName,const Char_t *inName,
- const Char_t *outName,Int_t *skipEvt,Int_t n);
-
-void WriteVtx(const Char_t *name,Double_t *zvtx,Int_t n);
+//=============================================================================
-void ReadVtx(const Char_t *name,Double_t *zvtx,Int_t n);
+// number of events to be processed
+Int_t gNevents;
+// magnetic field
+Double_t gBfieldValue;
-//=============================================================================
+void AliBarrelRec_TPCparam(Int_t n=-1,Char_t mode='A') {
-void AliBarrelRec_TPCparam(TString mode="A",Int_t n=5000) {
+ //---------------------------------------------------------------------
+ // CONFIGURATION
+ //
+ // _Magnetic_field_
+ gBfieldValue = 0.4;
+ //
+ // _Type_of_collision_ (needed for TPC tracking parameterization)
+ // Available choices: !!! ONLY B = 0.4 TESLA !!!
+ // collcode = 0 -> PbPb6000 (HIJING with b<2fm)
+ // collcode = 1 -> low multiplicity: pp or pA
+ Int_t collcode = 1;
+ //
+ // _ITS_clusters_reconstruction_
+ // Available choices: (from AliITStestV2.C)
+ // SlowOrFast = 's' slow points
+ // SlowOrFast = 'f' fast points
+ Char_t SlowOrFast = 'f';
+ //
+ // _Primary_vertex_for_ITS_tracking_
+ // Available choices:
+ // Vtx4Tracking = 'H' from event Header
+ // --- for Pb-Pb ---
+ // Vtx4Tracking = 'S' from event header + Smearing
+ // (x=15,y=15,z=10) micron
+ // --- for pp/pA ---
+ // Vtx4Tracking = 'P' z from pixels, x,y in(0,0)
+ Char_t Vtx4Tracking = 'P';
+ // _Primary_vertex_for_analysis_ (AliITSVertex stored in tracks file)
+ // Available choices:
+ // Vtx4Analysis = 'C' Copy the same used for tracking
+ // --- for pp/pA ---
+ // Vtx4Analysis = 'T' x,y,z from Tracks
+ Char_t Vtx4Analysis = 'T';
+ //
+ // END CONFIGURATION
+ //---------------------------------------------------------------------
const Char_t *name=" AliBarrelRec_TPCparam";
- cerr<<'\n'<<name<<"...\n";
+ printf("\n %s\n",name);
gBenchmark->Start(name);
+ if(n==-1) { // read number of events to be processed from file
+ TFile *f = new TFile("galice.root");
+ gAlice = (AliRun*)f->Get("gAlice");
+ n = gAlice->GetEventsPerRun();
+ delete gAlice;
+ gAlice=0;
+ f->Close();
+ delete f;
+ printf(" All %d events in file will be processed\n",n);
+ }
+ gNevents = n;
- // filenames
- const Char_t *galiceName = "galice.root";
- const Char_t *TPCtrkNameS = "AliTPCtracksParam.root";
- const Char_t *ITSclsName = "AliITSclustersV2.root";
- const Char_t *ITStrkName = "AliITStracksV2.root";
- const Char_t *ITStrkName2 = "AliITStracksV2_2.root";
- const Char_t *ITSrefName = "ITStracksRefFile.root";
- // set here the code for the type of collision (needed for TPC tracking
- // parameterization). available collisions:
- //
- // coll = 0 -> PbPb6000 (HIJING with b<2fm)
- // coll = 1 -> pp
- const Int_t collcode = 1;
- const Bool_t slowVtx = kFALSE;
- const Bool_t retrack = kFALSE;
- // set here the value of the magnetic field
- const Double_t BfieldValue = 0.4;
-
- // set conversion constant for kalman tracks
- AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue);
-
- Bool_t clust = kTRUE;
-
- //################ SWITCH #############################################
- switch (*(mode.Data())) {
- //############## case A #############################################
+ // conversion constant for kalman tracks
+ AliKalmanTrack::SetConvConst(100/0.299792458/gBfieldValue);
+
+ // Selection of execution mode
+ switch(mode) {
case 'A':
- // ********** Build TPC tracks with parameterization *********** //
- TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n);
+ // Build TPC tracks with parameterization
+ TPCParamTracks(collcode,gBfieldValue);
- // ********** ITS RecPoints ************************************ //
- ITSHits2FastRecPoints(galiceName,n);
+ // ITS clusters
+ ITSFindClustersV2(SlowOrFast);
- break;
- //############## case B #############################################
- case 'B':
- cerr<<" ---> only tracking in ITS <---"<<endl;
+ // Vertex for ITS tracking
+ PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
- if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat",n)) return;
+ break;
+
+ case 'B':
+ printf(" ---> only tracking in ITS <---\n");
- if(!UpdateEvtsToSkip("itsclusters.log","evtsToSkip.dat",n)) clust=kFALSE;
+ // Update list of events to be skipped
+ if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
break;
- //############## END SWITCH #########################################
}
// Mark events that have to be skipped (if any)
- Int_t *skipEvt = new Int_t[n];
- for(Int_t i=0;i<n;i++) skipEvt[i] = 0;
- if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists)) {
+ Int_t *skipEvt = new Int_t[gNevents];
+ for(Int_t i=0; i<gNevents; i++) skipEvt[i] = 0;
+ if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists))
MarkEvtsToSkip("evtsToSkip.dat",skipEvt);
- }
-
- // ********** Find ITS clusters ******************************** //
- if(clust) ITSFindClustersV2(galiceName,ITSclsName,skipEvt,n);
-
- // ********** Tracking in ITS ********************************** //
- TString vtxMode;
- switch (collcode) {
- case 0: // ***** Pb-Pb *****
- ZvtxFromHeader(galiceName,n);
- vtxMode="H";
- ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
- skipEvt,vtxMode,n);
- break;
- case 1: // ***** pp *****
- if(slowVtx) {
- ZvtxFromSPD(galiceName,n);
- vtxMode="P";
- ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
- skipEvt,vtxMode,n);
- ZvtxFromTracks(ITStrkName,n);
- if(retrack) {
- vtxMode="T";
- ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName2,
- skipEvt,vtxMode,n);
- }
- } else {
- ZvtxFastpp(galiceName,n);
- vtxMode="F";
- ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
- skipEvt,vtxMode,n);
- }
- break;
- }
-
+
+ // Tracking in ITS
+ ITSFindTracksV2(skipEvt);
- // ********** Make ITS tracks reference file ******************* //
- ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,skipEvt,n);
+ // Vertex for analysis
+ PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
+ // Create ITS tracks reference file
+ ITSMakeRefFile(skipEvt);
delete [] skipEvt;
return;
}
//-----------------------------------------------------------------------------
-Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName,Int_t n) {
-
- if(!gSystem->AccessPathName(logName,kFileExists)) {
- FILE *ifile = fopen(logName,"r");
- Int_t lEvt=0,nCol=1;
- while(nCol>0) {
- nCol = fscanf(ifile,"%d",&lEvt);
- }
- fclose(ifile);
- if(lEvt==n) {
- cerr<<" All events already reconstructed "<<endl;
- return 0;
- } else {
- FILE *ofile = fopen("evtsToSkip.dat","a");
- fprintf(ofile,"%d\n",lEvt);
- fclose(ofile);
- }
- } else {
- cerr<<"File itstracking.log not found"<<endl;
- }
+void CopyVtx(const Char_t *inName,const Char_t *outName) {
- return 1;
-}
-//-----------------------------------------------------------------------------
-Int_t MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
+ // Open input and output files
+ TFile *inFile = new TFile(inName);
+ TFile *outFile = new TFile(outName,"update");
- cerr<<"\n*******************************************************************\n";
- cerr<<"\nChecking for events to skip...\n";
+ TDirectory *curdir;
+ Char_t vname[20];
- Int_t evt,ncol;
- FILE *f = fopen(evtsName,"r");
- while(1) {
- ncol = fscanf(f,"%d",&evt);
- if(ncol<1) break;
- skipEvt[evt] = 1;
- cerr<<" ev. "<<evt<<" will be skipped\n";
+ for(Int_t ev=0; ev<gNevents; ev++) {
+ sprintf(vname,"Vertex_%d",ev);
+ AliITSVertex *vertex = (AliITSVertex*)inFile->Get(vname);
+ if(!vertex) continue;
+ curdir = gDirectory;
+ outFile->cd();
+ vertex->Write();
+ curdir->cd();
+ vertex = 0;
}
- fclose(f);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
-Int_t TPCParamTracks(const Char_t *galiceName,const Char_t *outName,
- const Int_t coll,const Double_t Bfield,Int_t n) {
-
- cerr<<"\n*******************************************************************\n";
-
- const Char_t *name="TPCParamTracks";
- cerr<<'\n'<<name<<"...\n";
- gBenchmark->Start(name);
-
- TFile *outFile=TFile::Open(outName,"recreate");
- TFile *inFile =TFile::Open(galiceName);
-
- AliTPCtrackerParam tracker(coll,Bfield,n);
- tracker.BuildTPCtracks(inFile,outFile);
-
- delete gAlice; gAlice=0;
inFile->Close();
outFile->Close();
+ delete inFile;
+ delete outFile;
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return 0;
+ return;
}
//-----------------------------------------------------------------------------
-Int_t ITSHits2FastRecPoints(const Char_t *galiceName,Int_t n) {
-
- cerr<<"\n*******************************************************************\n";
+void ITSFindClustersV2(Char_t SlowOrFast) {
- const Char_t *name="ITSHits2FastRecPoints";
- cerr<<'\n'<<name<<"...\n";
- gBenchmark->Start(name);
-
- Int_t nsignal=25;
- Int_t size=-1;
-
- // Connect the Root Galice file containing Geometry, Kine and Hits
- TFile *file = TFile::Open(galiceName,"UPDATE");
-
- gAlice = (AliRun*)file->Get("gAlice");
-
-
- AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
- if(!ITS) return 1;
-
- // Set the simulation model
- for(Int_t i=0;i<3;i++) {
- ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
- }
-
-
- //
- // Event Loop
- //
- for(Int_t ev=0; ev<n; ev++) {
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
- Int_t nparticles = gAlice->GetEvent(ev);
- cerr<<"Number of particles: "<<nparticles<<endl;
- gAlice->SetEvent(ev);
- if(!gAlice->TreeR()) gAlice->MakeTree("R");
- if(!gAlice->TreeR()->GetBranch("ITSRecPointsF")) {
- ITS->MakeBranch("RF");
- if(nparticles <= 0) return 1;
-
- Int_t bgr_ev=Int_t(ev/nsignal);
- //printf("bgr_ev %d\n",bgr_ev);
- ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
- }
- } // event loop
-
- delete gAlice; gAlice=0;
- file->Close();
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
-Int_t ITSFindClustersV2(const Char_t *galiceName,const Char_t *outName,
- Int_t *skipEvt,Int_t n) {
- //
- // This function converts AliITSRecPoint(s) to AliITSclusterV2
- //
- cerr<<"\n*******************************************************************\n";
+ printf("\n------------------------------------\n");
const Char_t *name="ITSFindClustersV2";
- cerr<<'\n'<<name<<"...\n";
+ printf("\n %s\n",name);
gBenchmark->Start(name);
-
- TFile *inFile=TFile::Open(galiceName);
- if(!inFile->IsOpen()) { cerr<<"Can't open galice.root !\n"; return 1; }
-
- if(!(gAlice=(AliRun*)inFile->Get("gAlice"))) {
- cerr<<"Can't find gAlice !\n";
- return 1;
- }
- AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
- if(!ITS) { cerr<<"Can't find the ITS !\n"; return 1; }
- AliITSgeom *geom=ITS->GetITSgeom();
-
- TFile *outFile = TFile::Open(outName,"recreate");
- geom->Write();
-
- // open logfile for done events
- FILE *logfile = fopen("itsclusters.log","w");
-
- // loop on events
- for(Int_t ev=0; ev<n; ev++) {
- // write to logfile of begun events
- fprintf(logfile,"%d\n",ev);
-
- if(skipEvt[ev]) continue;
-
- cerr<<"--- Processing event "<<ev<<" ---"<<endl;
- inFile->cd();
- gAlice->GetEvent(ev);
-
- TClonesArray *clusters = new TClonesArray("AliITSclusterV2",10000);
- Char_t cname[100];
- sprintf(cname,"TreeC_ITS_%d",ev);
- TTree *cTree = new TTree(cname,"ITS clusters");
- cTree->Branch("Clusters",&clusters);
-
- TTree *pTree=gAlice->TreeR();
- if(!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
- TBranch *branch = pTree->GetBranch("ITSRecPointsF");
- if(!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1; }
- TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
- branch->SetAddress(&points);
-
- TClonesArray &cl=*clusters;
- Int_t nclusters=0;
- Int_t nentr=(Int_t)branch->GetEntries();
-
- cerr<<"Number of entries in TreeR_RF: "<<nentr<<endl;
-
- // Float_t *lp = new Float_t[5];
- // Int_t *lab = new Int_t[6];
- Float_t lp[5];
- Int_t lab[6];
-
- for(Int_t i=0; i<nentr; i++) {
- points->Clear();
- branch->GetEvent(i);
- Int_t ncl=points->GetEntriesFast(); if(ncl==0){cTree->Fill();continue;}
- Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
- Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
- Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
- Double_t yshift = x*rot[0] + y*rot[1];
- Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
- nclusters+=ncl;
-
- Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD
- if(lay==4 || lay==3) kmip=280.;
- if(lay==6 || lay==5) kmip= 38.;
-
- for(Int_t j=0; j<ncl; j++) {
- AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
- lp[0]=-p->GetX()-yshift; if(lay==1) lp[0]=-lp[0];
- lp[1]=p->GetZ()+zshift;
- lp[2]=p->GetSigmaX2();
- lp[3]=p->GetSigmaZ2();
- lp[4]=p->GetQ(); lp[4]/=kmip;
- lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
- lab[3]=ndet;
-
- Int_t label=lab[0];
- if(label>=0) {
- TParticle *part=(TParticle*)gAlice->Particle(label);
- label=-3;
- while (part->P() < 0.005) {
- Int_t m=part->GetFirstMother();
- if(m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
- label=m;
- part=(TParticle*)gAlice->Particle(label);
- }
- if (lab[1]<0) lab[1]=label;
- else if(lab[2]<0) lab[2]=label;
- else cerr<<"No empty labels !\n";
- }
- new(cl[j]) AliITSclusterV2(lab,lp);
- }
- cTree->Fill(); clusters->Delete();
- points->Delete();
- }
-
- outFile->cd();
- cTree->Write();
-
- cerr<<"Number of clusters: "<<nclusters<<endl;
- //delete [] lp;
- //delete [] lab;
- delete cTree; delete clusters; delete points;
-
- } // end loop on events
-
- fprintf(logfile,"%d\n",n); //this means all evts are successfully completed
- fclose(logfile);
-
- delete gAlice; gAlice=0;
-
- inFile->Close();
- outFile->Close();
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
-Int_t ZvtxFromHeader(const Char_t *galiceName,Int_t n) {
-
- cerr<<"\n*******************************************************************\n";
-
- const Char_t *name="ZvtxFromHeader";
- cerr<<'\n'<<name<<"...\n";
- gBenchmark->Start(name);
-
- TFile *galice = TFile::Open(galiceName);
-
- Double_t *zvtx = new Double_t[n];
-
- for(Int_t ev=0; ev<n; ev++){
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
-
- TArrayF o = 0;
- o.Set(3);
- AliHeader* header = 0;
- TTree* treeE = (TTree*)gDirectory->Get("TE");
- treeE->SetBranchAddress("Header",&header);
- treeE->GetEntry(ev);
- AliGenEventHeader* genHeader = header->GenEventHeader();
- if(genHeader) {
- // get primary vertex position
- genHeader->PrimaryVertex(o);
- // set position of primary vertex
- zvtx[ev] = (Double_t)o[2];
- } else {
- cerr<<" ! event header not found : setting z vertex to 0 !"<<endl;
- zvtx[ev] = 0.;
- }
- delete header;
- }
-
- galice->Close();
-
- // Write vertices to file
- WriteVtx("vtxHeader.root",zvtx,n);
-
- delete [] zvtx;
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
-Int_t ZvtxFastpp(const Char_t *galiceName,Int_t n) {
-
- cerr<<"\n*******************************************************************\n";
-
- const Char_t *name="ZvtxFastpp";
- cerr<<'\n'<<name<<"...\n";
- gBenchmark->Start(name);
-
- TFile *galice = TFile::Open(galiceName);
- AliRun *gAlice = (AliRun*)galice->Get("gAlice");
-
- Double_t *zvtx = new Double_t[n];
-
- TParticle *part;
- Int_t nPart,b;
- Double_t dNchdy,E,theta,eta,zvtxTrue,sigmaz,probVtx;
-
- for(Int_t ev=0; ev<n; ev++) {
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
-
- TArrayF o = 0;
- o.Set(3);
- AliHeader* header = 0;
- TTree* treeE = (TTree*)galice->Get("TE");
- treeE->SetBranchAddress("Header",&header);
- treeE->GetEntry(ev);
- AliGenEventHeader* genHeader = header->GenEventHeader();
- if(genHeader) {
- // get primary vertex position
- genHeader->PrimaryVertex(o);
- // set position of primary vertex
- zvtxTrue = (Double_t)o[2];
- } else {
- cerr<<" ! event header not found : setting z vertex to 0 !"<<endl;
- zvtxTrue = 0.;
- }
- delete header;
-
-
- // calculate dNch/dy pf the event (charged primaries in |eta|<0.5)
- nPart = (Int_t)gAlice->GetEvent(ev);
- dNchdy = 0.;
- for(b=0; b<nPart; b++) {
- part = (TParticle*)gAlice->Particle(b);
- if(TMath::Abs(part->GetPdgCode())<10) continue; // reject partons
- if(TMath::Abs(part->Vx()*part->Vx()+part->Vy()*part->Vy())>0.01) continue; // reject secondaries
- E = part->Energy();
- if(E>6900.) continue; // reject incoming protons
- theta = part->Theta();
- if(theta<.1 || theta>3.) continue;
- eta = -TMath::Log(TMath::Tan(theta/2.));
- if(TMath::Abs(eta)>0.5) continue; // count particles in |eta|<0.5
- dNchdy+=1.;
- }
-
- // get sigma(z) corresponding to the mult of this event
- if(dNchdy>1.) {
- sigmaz = 0.0417/TMath::Sqrt(dNchdy);
- } else {
- sigmaz = 0.0500;
- }
- // smear the original position of the primary vertex
- zvtx[ev] = gRandom->Gaus(zvtxTrue,sigmaz);
-
- // compute the probability that the vertex is found
- probVtx = 1.;
- if(dNchdy<24.) probVtx = 0.85+0.00673*dNchdy;
- if(dNchdy<14.) probVtx = 0.85;
-
- if(gRandom->Rndm()>probVtx) zvtx[ev] = -1000.;// no zvtx for this event
-
- }
- galice->Close();
-
- // Write vertices to file
- WriteVtx("vtxFastpp.root",zvtx,n);
-
- delete [] zvtx;
- //delete gAlice; gAlice=0;
-
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
-Int_t ZvtxFromSPD(const Char_t *galiceName,Int_t n) {
-
-
- Double_t *zvtx = new Double_t[n];
-
- Bool_t debug = kFALSE;
- ofstream fildeb;
- if(debug)fildeb.open("FindZV.out");
- ofstream filou;
- //filou.open("zcoor.dat");
- const Float_t kPi=TMath::Pi();
- const Int_t kFirstL1=0;
- const Int_t kLastL1=79;
- const Int_t kFirstL2=80;
- const Int_t kLastL2=239;
- // const Float_t kDiffPhiMax=0.0005;
- const Float_t kDiffPhiMax=0.01;
- const Float_t kphl=0.;
- const Float_t kphh=kPi/18.;
- TDirectory *currdir = gDirectory;
- TH2F *hz2z1 = 0;
- TH1F *zvdis = 0;
- TH2F *dpvsz = 0;
- TProfile *scoc = new TProfile("scoc","Zgen - Zmeas vs occ lay 1",5,10.,60.);
- TProfile *prof = new TProfile("prof","z2 vs z1",50,-15.,15.);
- TH1F *phdiff = new TH1F("phdiff","#Delta #phi",100,kphl,kphh);
- TH1F *phdifsame = new TH1F("phdifsame","#Delta #phi same track",100,kphl,kphh);
- if(gAlice){delete gAlice; gAlice = 0;}
- TFile *file = TFile::Open(galiceName);
- gAlice = (AliRun*)file->Get("gAlice");
-
- AliITS *ITS = (AliITS*)gAlice->GetModule("ITS");
- AliITSgeom *geom = ITS->GetITSgeom();
- TTree *TR=0;
- Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
- Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
- Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
- Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
- char name[30];
- Float_t avesca=0;
- Float_t avesca2=0;
- Int_t N=0;
- for(Int_t event=0; event<n; event++){
- sprintf(name,"event_%d",event);
- hz2z1=new TH2F(name,"z2 vs z1",50,-15.,15.,50,-15.,15.);
- hz2z1->SetDirectory(currdir);
- gAlice->GetEvent(event);
- TParticle * part = gAlice->Particle(0);
- Float_t oriz=part->Vz();
- cout<<"\n ==============================================================\n";
- cout<<" Processing event "<<event<<" "<<oriz<<endl;
- cout<<" ==============================================================\n";
- if(debug){
- fildeb<<"\n ==============================================================\n";
- fildeb<<" Processing event "<<event<<" "<<oriz<<endl;
- fildeb<<" ==============================================================\n";
- }
- file->cd();
- TR = gAlice->TreeR();
- if(!TR) cerr<<"TreeR not found\n";
- TBranch *branch=TR->GetBranch("ITSRecPointsF");
- if(!branch) cerr<<"Branch ITSRecPointsF not found\n";
- TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
- branch->SetAddress(&points);
- Int_t nentr=(Int_t)branch->GetEntries();
- Float_t zave=0;
- Float_t rmszav=0;
- Float_t zave2=0;
- Int_t firipixe=0;
- cout<<endl;
- for(Int_t module= kFirstL1; module<=kLastL1;module++){
- points->Clear();
- branch->GetEvent(module);
- Int_t nrecp1 = points->GetEntriesFast();
- //cerr<<nrecp1<<endl;
- for(Int_t i=0; i<nrecp1;i++){
- AliITSRecPoint *current = (AliITSRecPoint*)points->UncheckedAt(i);
- lc[0]=current->GetX();
- lc[2]=current->GetZ();
- geom->LtoG(module,lc,gc);
- zave+=gc[2];
- zave2+=gc[2]*gc[2];
- firipixe++;
- }
- points->Delete();
- }
- if(firipixe>1){
- rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
- zave=zave/firipixe;
- cout<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl;
- if(debug)fildeb<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl;
- }
- else {
- cout<<"No rec points on first layer for this event\n";
- if(debug)fildeb<<"No rec points on first layer for this event\n";
- }
- Float_t zlim1=zave-rmszav;
- Float_t zlim2=zave+rmszav;
- Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
- zlim2=zlim1 + sepa/10.;
- sprintf(name,"pz_ev_%d",event);
- dpvsz=new TH2F(name,"#Delta #phi vs Z",sepa,zlim1,zlim2,100,0.,0.03);
- dpvsz->SetDirectory(currdir);
- sprintf(name,"z_ev_%d",event);
- zvdis=new TH1F(name,"zv distr",sepa,zlim1,zlim2);
- zvdis->SetDirectory(currdir);
- cout<<"Z limits: "<<zlim1<<" "<<zlim2<<endl;
- if(debug)fildeb<<"Z limits: "<<zlim1<<" "<<zlim2<<endl;
- Int_t sizarr=100;
- TArrayF *zval = new TArrayF(sizarr);
- Int_t ncoinc=0;
- for(Int_t module= kFirstL1; module<=kLastL1;module++){
- //cout<<"processing module "<<module<<" \r";
- branch->GetEvent(module);
- Int_t nrecp1 = points->GetEntriesFast();
- TObjArray *poiL1 = new TObjArray(nrecp1);
- for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(points->UncheckedAt(i),i);
- //ITS->ResetRecPoints();
- points->Delete();
- for(Int_t i=0; i<nrecp1;i++){
- AliITSRecPoint *current = (AliITSRecPoint*)poiL1->UncheckedAt(i);
- lc[0]=current->GetX();
- lc[2]=current->GetZ();
- geom->LtoG(module,lc,gc);
- Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
- Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
- if(phi1<0)phi1=2*kPi+phi1;
- Int_t lab1 = current->GetLabel(0);
- Float_t phi1d=phi1*180./kPi;
- //cerr<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1d<<" "<<lab1<<" \n";
- for(Int_t modul2=kFirstL2; modul2<=kLastL2; modul2++){
- branch->GetEvent(modul2);
- Int_t nrecp2 = points->GetEntriesFast();
- for(Int_t j=0; j<nrecp2;j++){
- AliITSRecPoint *recp = (AliITSRecPoint*)points->UncheckedAt(j);
- lc2[0]=recp->GetX();
- lc2[2]=recp->GetZ();
- geom->LtoG(modul2,lc2,gc2);
- Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
- Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
- Int_t lab2 = recp->GetLabel(0);
- Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
- if(phi2<0)phi2=2.*kPi+phi2;
- Float_t diff = TMath::Abs(phi2-phi1);
- if(diff>kPi)diff=2.*kPi-diff;
- if(lab1==lab2)phdifsame->Fill(diff);
- if(zr0>zlim1 && zr0<zlim2){
- phdiff->Fill(diff);
- dpvsz->Fill(zr0,diff);
- if(diff<kDiffPhiMax ){
- zvdis->Fill(zr0);
- zval->AddAt(zr0,ncoinc);
- ncoinc++;
- if(ncoinc==(sizarr-1)){
- sizarr+=100;
- zval->Set(sizarr);
- }
- Float_t phi2d=phi2*180./kPi;
- Float_t diffd=diff*180./kPi;
- //cerr<<"module2 "<<modul2<<" "<<gc2[0]<<" "<<gc2[1]<<" "<<gc2[2]<<" "<<phi2d<<" "<<diffd<<" "<<lab2<<" \n";
- // cout<<"zr0= "<<zr0<<endl;
- hz2z1->Fill(gc[2],gc2[2]);
- prof->Fill(gc[2]-oriz,gc2[2]-oriz);
- }
- }
- }
- points->Delete();
- }
- }
- delete poiL1;
- } // loop on modules
- cout<<endl;
- Float_t zcmp=0;
- Float_t zsig=0;
- EvalZ(zvdis,sepa,zcmp,zsig,ncoinc,zval,&fildeb);
- if(zcmp!=-100 && zsig!=-100){
- N++;
- Float_t scarto = (oriz-zcmp)*10000.;
- Float_t ascarto=TMath::Abs(scarto);
- scoc->Fill(firipixe,ascarto);
- avesca+=ascarto;
- avesca2+=scarto*scarto;
- if(debug)fildeb<<"Mean value: "<<zcmp<<" +/-"<<zsig<<" Zgen- Zmeas (micron)= "<<scarto<<endl;
- cout<<"Mean value: "<<zcmp<<" RMS: "<<zsig<<" Zgen- Zmeas (micron)= "<<scarto<<endl;
- //filou<<event<<" "<<zcmp<<" "<<zsig<<" "<<oriz<<" "<<scarto<<endl;
- zvtx[event] = (Double_t)zcmp;
- }
- else {
- cout<<"Not enough points in event "<<event<<endl;
- if(debug)fildeb<<"Not enough points\n";
- //filou<<event<<" "<<-1000.<<" "<<0.<<" "<<oriz<<" "<<0.<<endl;
- zvtx[event] = -1000.;
- }
- delete zval;
- } // loop on events
- file->Close();
- //hz2z1->Draw("box");
-
- // Write vertices to file
- WriteVtx("vtxSPD.root",zvtx,n);
- delete [] zvtx;
-
- if(N>1) {
- avesca2=TMath::Sqrt((avesca2/(N-1)-avesca*avesca/N/(N-1))/N);
- avesca=avesca/N;
+ //--- taken from AliITStestV2.C--------------------------------------
+ //
+ if (SlowOrFast=='f') {
+ //cerr<<"Fast AliITSRecPoint(s) !\n";
+ //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
+ //AliITSHits2FastRecPoints();
} else {
- avesca2 = 0;
- }
- cout<<"\n \n ==========================================================\n";
- cout<<"Number of events with vertex estimate "<<N<<endl;
- cout<<"Average of the (abs) Zgen-Zmeas : "<<avesca;
- cout<<" +/- "<<avesca2<<" micron"<<endl;
- if(debug){
- fildeb<<"\n \n ==========================================================\n";
- fildeb<<"Number of events with vertex estimate "<<N<<endl;
- fildeb<<"Average of the (abs) Zgen-Zmeas : "<<avesca;
- fildeb<<" +/- "<<avesca2<<endl;
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C");
+ AliITSHits2SDigits();
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C");
+ AliITSSDigits2Digits();
+ //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C");
+ //AliITSDigits2RecPoints();
}
- return 0;
-}
-
-
-//-----------------------------------------------------------------------------
-void EvalZ(TH1F *hist,Int_t sepa,Float_t &av, Float_t &sig,Int_t ncoinc,
- TArrayF *zval,ofstream *deb) {
-
- av=0;
- sig=0;
- Int_t N=0;
- Int_t NbinNotZero=0;
- for(Int_t i=1;i<=sepa;i++){
- Float_t cont=hist->GetBinContent(i);
- av+=cont;
- sig+=cont*cont;
- N++;
- if(cont!=0)NbinNotZero++;
- }
- if(NbinNotZero==0){av=-100; sig=-100; return;}
- if(NbinNotZero==1){sig=-100; return;}
- sig=TMath::Sqrt(sig/(N-1)-av*av/N/(N-1));
- av=av/N;
- Float_t val1=hist->GetBinLowEdge(sepa);
- Float_t val2=hist->GetBinLowEdge(1);
- for(Int_t i=1;i<=sepa;i++){
- Float_t cont=hist->GetBinContent(i);
- if(cont>(av+sig*2.)){
- Float_t curz=hist->GetBinLowEdge(i);
- if(curz<val1)val1=curz;
- if(curz>val2)val2=curz;
- }
- }
- val2+=hist->GetBinWidth(1);
- cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
- if(deb)(*deb)<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
- av=0;
- sig=0;
- N=0;
- for(Int_t i=0; i<ncoinc; i++){
- Float_t z=zval->At(i);
- if(z<val1)continue;
- if(z>val2)continue;
- av+=z;
- sig+=z*z;
- N++;
- }
- if(N<=1){av=-100; sig=-100; return;}
- sig=TMath::Sqrt((sig/(N-1)-av*av/N/(N-1))/N);
- av=av/N;
- return;
-}
-//-----------------------------------------------------------------------------
-Int_t ZvtxFromTracks(const Char_t *trkName,Int_t n) {
-
- cerr<<"\n*******************************************************************\n";
-
- const Char_t *name="ZvtxFromTracks";
- cerr<<'\n'<<name<<"...\n";
- gBenchmark->Start(name);
-
- Double_t *zvtx = new Double_t[n];
- Double_t *zvtxSPD = new Double_t[n];
- ReadVtx("vtxSPD.root",zvtxSPD,n);
-
- TFile *itstrk = TFile::Open(trkName);
-
-
- Int_t nTrks,i,nUsedTrks;
- Double_t pt,d0rphi,d0z;
- Double_t zvtxTrks,ezvtxTrks,sumWeights;
-
- for(Int_t ev=0; ev<n; ev++){
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
-
- // Tree with ITS tracks
- Char_t tname[100];
- sprintf(tname,"TreeT_ITS_%d",ev);
-
- TTree *tracktree=(TTree*)itstrk->Get(tname);
- if(!tracktree) continue;
- AliITStrackV2 *itstrack=new AliITStrackV2;
- tracktree->SetBranchAddress("tracks",&itstrack);
- nTrks = (Int_t)tracktree->GetEntries();
- nUsedTrks = 0;
- zvtxTrks = 0.;
- ezvtxTrks = 0.;
- sumWeights = 0.;
-
- //cerr<<" ITS tracks: "<<nTrks<<endl;
- // loop on tracks
- for(i=0; i<nTrks; i++) {
- tracktree->GetEvent(i);
- // get pt and impact parameters
- pt = 1./TMath::Abs(itstrack->Get1Pt());
- d0rphi = TMath::Abs(itstrack->GetD());
- itstrack->PropagateToVertex();
- d0z = itstrack->GetZ();
- // check if the track is from (0,0) in (x,y)
- if(!VtxTrack(pt,d0rphi)) continue;
- nUsedTrks++;
- zvtxTrks += d0z/d0zRes(pt)/d0zRes(pt);
- sumWeights += 1./d0zRes(pt)/d0zRes(pt);
-
- } // loop on tracks
-
- if(nUsedTrks>1) {
- // estimated position in z of vertex
- zvtxTrks /= sumWeights;
- // estimated error
- ezvtxTrks = TMath::Sqrt(1./sumWeights);
- } else {
- zvtxTrks = zvtxSPD[ev];
- }
-
- zvtx[ev] = zvtxTrks;
-
- } // loop on events
-
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
+ AliITSFindClustersV2(SlowOrFast,gNevents);
+ //
+ //--------------------------------------------------------------------
- // Write vertices to file
- WriteVtx("vtxTracks.root",zvtx,n);
- delete [] zvtx;
- delete [] zvtxSPD;
gBenchmark->Stop(name);
gBenchmark->Show(name);
- return 0;
-}
-//----------------------------------------------------------------------------
-Bool_t VtxTrack(Double_t pt,Double_t d0rphi) {
-
- Double_t d0rphiRes = TMath::Sqrt(11.59*11.59+65.76*65.76/TMath::Power(pt,1.878))/10000.;
- Double_t nSigma = 3.;
- if(d0rphi<nSigma*d0rphiRes) { return kTRUE; } else { return kFALSE; }
-}
-//----------------------------------------------------------------------------
-Double_t d0zRes(Double_t pt) {
- Double_t res = TMath::Sqrt(34.05*34.05+170.1*170.1/TMath::Power(pt,1.226));
- return res/10000.;
+ return;
}
//-----------------------------------------------------------------------------
-Int_t ITSFindTracksV2(const Char_t *galiceName,const Char_t *inName,
- const Char_t *inName2,const Char_t *outName,
- Int_t *skipEvt,TString vtxMode,Int_t n) {
+Int_t ITSFindTracksV2(Int_t *skipEvt) {
+ printf("\n------------------------------------\n");
- cerr<<"\n*******************************************************************\n";
-
const Char_t *name="ITSFindTracksV2";
- cerr<<'\n'<<name<<"...\n";
+ printf("\n %s\n",name);
gBenchmark->Start(name);
- // Read vertices from file
- Double_t *zvtx = new Double_t[n];
- Char_t *vtxfile="vtxHeader.root";
-
- cerr<<" Using vtxMode: ";
- switch (*(vtxMode.Data())) {
- case 'H':
- cerr<<" Header"<<endl;
- vtxfile = "vtxHeader.root";
- break;
- case 'F':
- cerr<<" fast pp"<<endl;
- vtxfile = "vtxFastpp.root";
- break;
- case 'P':
- cerr<<" SPD"<<endl;
- vtxfile = "vtxSPD.root";
- break;
- case 'T':
- cerr<<" Tracks"<<endl;
- vtxfile = "vtxTracks.root";
- break;
- }
-
- ReadVtx(vtxfile,zvtx,n);
-
-
- TFile *outFile = TFile::Open(outName,"recreate");
- TFile *inFile = TFile::Open(inName);
- TFile *inFile2 = TFile::Open(inName2);
+
+ TFile *outFile = new TFile("AliITStracksV2.root","recreate");
+ TFile *inTPCtrks = new TFile("AliTPCtracksParam.root");
+ TFile *inVertex = new TFile("Vtx4Tracking.root");
+ TFile *inClusters = new TFile("AliITSclustersV2.root");
- AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
- if(!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+ AliITSgeom *geom=(AliITSgeom*)inClusters->Get("AliITSgeom");
+ if(!geom) { printf("can't get ITS geometry !\n"); return;}
+ Double_t vtx[3];
Int_t flag1stPass,flag2ndPass;
+ Char_t vname[20];
// open logfile for done events
FILE *logfile = fopen("itstracking.log","w");
+ // Instantiate AliITStrackerV2
AliITStrackerV2 tracker(geom);
+
// loop on events
- for(Int_t ev=0; ev<n; ev++){
+ for(Int_t ev=0; ev<gNevents; ev++){
// write to logfile of begun events
fprintf(logfile,"%d\n",ev);
if(skipEvt[ev]) continue;
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
+ printf(" --- Processing event %d ---\n",ev);
- // pass event number to tracker
+ // pass event number to the tracker
tracker.SetEventNumber(ev);
// set position of primary vertex
- Double_t vtx[3];
- vtx[0]=0.; vtx[1]=0.; vtx[2]=zvtx[ev];
+ sprintf(vname,"Vertex_%d",ev);
+ AliITSVertex *vertex = (AliITSVertex*)inVertex->Get(vname);
+ if(vertex) {
+ vertex->GetXYZ(vtx);
+ delete vertex;
+ } else {
+ printf(" AliITSVertex not found for event %d\n",ev);
+ printf(" Using (0,0,0) for ITS tracking\n");
+ vtx[0] = vtx[1] = vtx[2] = 0.;
+ }
flag1stPass=1; // vtx constraint
flag2ndPass=0; // no vtx constraint
vtx[2]=0.;
}
- cerr<<"+++\n+++ Setting primary vertex z = "<<vtx[2]<<
- " cm for ITS tracking\n+++\n";
tracker.SetVertex(vtx);
// setup vertex constraint in the two tracking passes
tracker.SetupSecondPass(flags);
// find the tracks
- tracker.Clusters2Tracks(inFile,outFile);
+ tracker.Clusters2Tracks(inTPCtrks,outFile);
} // loop on events
- fprintf(logfile,"%d\n",n); //this means all evts are successfully completed
+ fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
fclose(logfile);
- inFile->Close();
- inFile2->Close();
+ delete geom;
+
+ inTPCtrks->Close();
+ inClusters->Close();
+ inVertex->Close();
outFile->Close();
- delete [] zvtx;
gBenchmark->Stop(name);
gBenchmark->Show(name);
- return 0;
+ return;
}
//-----------------------------------------------------------------------------
-Int_t ITSMakeRefFile(const Char_t *galice,const Char_t *inname,
- const Char_t *outname,Int_t *skipEvt,Int_t n) {
+void ITSMakeRefFile(Int_t *skipEvt) {
-
- cerr<<"\n*******************************************************************\n";
+ printf("\n------------------------------------\n");
- Int_t rc=0;
const Char_t *name="ITSMakeRefFile";
- cerr<<'\n'<<name<<"...\n";
+ printf("\n %s\n",name);
gBenchmark->Start(name);
- TFile *out = TFile::Open(outname,"recreate");
- TFile *trk = TFile::Open(inname);
- TFile *kin = TFile::Open(galice);
+ TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
+ TFile *trk = TFile::Open("AliITStracksV2.root");
+ TFile *kin = TFile::Open("galice.root");
// Get gAlice object from file
- if(!(gAlice=(AliRun*)kin->Get("gAlice"))) {
- cerr<<"gAlice has not been found on galice.root !\n";
- return 1;
- }
+ gAlice=(AliRun*)kin->Get("gAlice");
Int_t label;
TParticle *Part;
RECTRACK rectrk;
- for(Int_t ev=0; ev<n; ev++){
+ for(Int_t ev=0; ev<gNevents; ev++){
if(skipEvt[ev]) continue;
- cerr<<" --- Processing event "<<ev<<" ---"<<endl;
+ printf(" --- Processing event %d ---\n",ev);
gAlice->GetEvent(ev);
gBenchmark->Show(name);
- return rc;
+ return;
}
//-----------------------------------------------------------------------------
-void WriteVtx(const Char_t *name,Double_t *zvtx,Int_t n) {
+void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
- // Set Random Number seed
- TDatime dt;
- UInt_t curtime=dt.Get();
- UInt_t procid=gSystem->GetPid();
- UInt_t seed=curtime-procid;
- gRandom->SetSeed(seed);
+ printf("\n------------------------------------\n");
+ printf("\nChecking for events to skip...\n");
- Double_t xvtx,yvtx;
- Double_t sigmaVtxTransverse = 15.e-4;
- TVector3 *ioVtx = new TVector3;
- TTree *vtxtree = new TTree("TreeVtx","Tree with positions of primary vertex");
- vtxtree->Branch("vertices","TVector3",&ioVtx);
+ Int_t evt,ncol;
- for(Int_t ev=0; ev<n; ev++) {
- // x and y coordinates of the vertex are smeared according to a gaussian
- xvtx = gRandom->Gaus(0.,sigmaVtxTransverse);
- yvtx = gRandom->Gaus(0.,sigmaVtxTransverse);
+ FILE *f = fopen(evtsName,"r");
+ while(1) {
+ ncol = fscanf(f,"%d",&evt);
+ if(ncol<1) break;
+ skipEvt[evt] = 1;
+ printf(" event %d will be skipped\n",evt);
+ }
+ fclose(f);
- ioVtx->SetX(xvtx);
- ioVtx->SetY(yvtx);
- ioVtx->SetZ(zvtx[ev]);
+ return;
+}
+//-----------------------------------------------------------------------------
+void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
- //cerr<<"W ev "<<ev<<" ("<<ioVtx->X()<<","<<ioVtx->Y()<<","<<ioVtx->Z()<<")"<<endl;
+ printf("\n------------------------------------\n");
-
- vtxtree->Fill();
+ const Char_t *name="PrimaryVertex";
+ printf("\n %s\n",name);
+ gBenchmark->Start(name);
+
+ switch(vtxMode) {
+ case 'H':
+ printf(" ... from event header\n");
+ VtxFromHeader(outName,kFALSE);
+ break;
+ case 'S':
+ printf(" ... from event header + smearing\n");
+ VtxFromHeader(outName,kTRUE);
+ break;
+ case 'P':
+ printf(" ... z from pixels for pp/pA\n");
+ ZvtxFromSPD(outName);
+ break;
+ case 'T':
+ printf(" ... from tracks for pp/pA\n");
+ VtxFromTracks(outName);
+ break;
+ case 'C':
+ printf(" ... copied from Vtx4Tracking.root to AliITStracksV2.root\n");
+ CopyVtx("Vtx4Tracking.root",outName);
+ break;
}
- // write the tree to a file
- TFile *f = new TFile(name,"recreate");
- vtxtree->Write();
- f->Close();
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
return;
}
//-----------------------------------------------------------------------------
-void ReadVtx(const Char_t *name,Double_t *zvtx,Int_t n) {
+void TPCParamTracks(Int_t coll,Double_t Bfield) {
+
+ printf("\n------------------------------------\n");
+
+ const Char_t *name="TPCParamTracks";
+ printf("\n %s\n",name);
+ gBenchmark->Start(name);
+
+ TFile *outFile=TFile::Open("AliTPCtracksParam.root","recreate");
+ TFile *inFile =TFile::Open("galice.root");
+
+ AliTPCtrackerParam tracker(coll,Bfield,gNevents);
+ tracker.BuildTPCtracks(inFile,outFile);
+
+ delete gAlice; gAlice=0;
+
+ inFile->Close();
+ outFile->Close();
+
+ gBenchmark->Stop(name);
+ gBenchmark->Show(name);
+
+ return;
+}
+//-----------------------------------------------------------------------------
+Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName) {
+
+ if(!gSystem->AccessPathName(logName,kFileExists)) {
+ FILE *ifile = fopen(logName,"r");
+ Int_t lEvt=0,nCol=1;
+ while(nCol>0) {
+ nCol = fscanf(ifile,"%d",&lEvt);
+ }
+ fclose(ifile);
+ if(lEvt==gNevents) {
+ printf(" All events already reconstructed\n");
+ return 0;
+ } else {
+ FILE *ofile = fopen("evtsToSkip.dat","a");
+ fprintf(ofile,"%d\n",lEvt);
+ fclose(ofile);
+ }
+ } else {
+ printf("File itstracking.log not found\n");
+ }
- TFile *f = TFile::Open(name);
+ return 1;
+}
+//-----------------------------------------------------------------------------
+void VtxFromHeader(const Char_t *outName,Bool_t smear) {
- TVector3 *ioVtx = 0;
- TTree *vtxtree = (TTree*)f->Get("TreeVtx");
- vtxtree->SetBranchAddress("vertices",&ioVtx);
- Int_t entries = (Int_t)vtxtree->GetEntries();
+ TDatime t;
+ UInt_t seed = t.Get();
+ gRandom->SetSeed(seed);
- for(Int_t ev=0; ev<n; ev++) {
- vtxtree->GetEvent(ev);
+ TFile *galice = new TFile("galice.root");
+ TFile *outFile = new TFile(outName,"update");
- zvtx[ev] = ioVtx->Z();
+ TDirectory *curdir;
+ Double_t pos[3],sigma[3];
+ if(smear) {
+ sigma[0]=15.e-4;
+ sigma[1]=15.e-4;
+ sigma[2]=10.e-4;
+ } else {
+ sigma[0]=0.;
+ sigma[1]=0.;
+ sigma[2]=0.;
+ }
+ Char_t vname[20];
- //cerr<<"R ev "<<ev<<" ("<<ioVtx->X()<<","<<ioVtx->Y()<<","<<zvtx[ev]<<")"<<endl;
+ galice->cd();
+ for(Int_t ev=0; ev<gNevents; ev++){
+ printf(" event %d\n",ev);
+ sprintf(vname,"Vertex_%d",ev);
+ TArrayF o = 0;
+ o.Set(3);
+ AliHeader* header = 0;
+ TTree* treeE = (TTree*)gDirectory->Get("TE");
+ treeE->SetBranchAddress("Header",&header);
+ treeE->GetEntry(ev);
+ AliGenEventHeader* genHeader = header->GenEventHeader();
+ if(genHeader) {
+ // get primary vertex position
+ genHeader->PrimaryVertex(o);
+ pos[0] = (Double_t)o[0];
+ pos[1] = (Double_t)o[1];
+ pos[2] = (Double_t)o[2];
+ if(smear) {
+ pos[0] = gRandom->Gaus(pos[0],sigma[0]);
+ pos[1] = gRandom->Gaus(pos[1],sigma[1]);
+ pos[2] = gRandom->Gaus(pos[2],sigma[2]);
+ }
+ // create AliITSVertex
+ AliITSVertex *vertex = new AliITSVertex(pos,sigma,vname);
+ } else {
+ printf(" ! event header not found : setting vertex to (0,0,0) !");
+ pos[0] = 0.;
+ pos[1] = 0.;
+ pos[2] = 0.;
+ // create AliITSVertex
+ AliITSVertex *vertex = new AliITSVertex(pos,sigma,vname);
+ }
+ delete header;
+ // write AliITSVertex to file
+ curdir = gDirectory;
+ outFile->cd();
+ if(smear) {
+ vertex->SetTitle("vertex from header, smeared");
+ } else {
+ vertex->SetTitle("vertex from header");
+ }
+ vertex->Write();
+ curdir->cd();
+ vertex = 0;
}
- f->Close();
+ outFile->Close();
+ galice->Close();
+
+ delete outFile;
+ delete galice;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void VtxFromTracks(const Char_t *outName) {
+
+ // Open input and output files
+ TFile *inFile = new TFile("AliITStracksV2.root");
+ TFile *outFile = new TFile(outName,"update");
+
+ // set AliRun object to 0
+ if(gAlice) gAlice = 0;
+
+ // Create vertexer
+ AliITSVertexerTracks *vertexer =
+ new AliITSVertexerTracks(inFile,outFile,gBfieldValue);
+ vertexer->SetFirstEvent(0);
+ vertexer->SetLastEvent(gNevents-1);
+ vertexer->SetDebug(0);
+ vertexer->PrintStatus();
+ // Find vertices
+ vertexer->FindVertices();
+
+ delete vertexer;
+
+ inFile->Close();
+ outFile->Close();
+ delete inFile;
+ delete outFile;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void ZvtxFromSPD(const Char_t *outName) {
+
+ // create fast RecPoints, which are used for vertex finding
+ cerr<<"Fast AliITSRecPoint(s) !\n";
+ gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
+ AliITSHits2FastRecPoints(0,gNevents-1);
+
+ // delphi ---> azimuthal range to accept tracklets
+ // window ---> window in Z around the peak of tracklets proj. in mm
+ Float_t delphi=0.05;
+ Float_t window=3.;
+ Float_t initx=0.;
+ Float_t inity=0.;
+
+ TFile *infile = new TFile("galice.root");
+ TFile *outfile = new TFile(outName,"update");
+
+ AliITSVertexerPPZ *vertexer = new AliITSVertexerPPZ(infile,outfile,initx,inity);
+ vertexer->SetFirstEvent(0);
+ vertexer->SetLastEvent(gNevents-1);
+ vertexer->SetDebug(0);
+ vertexer->SetDiffPhiMax(delphi);
+ vertexer->SetWindow(window);
+ vertexer->PrintStatus();
+ vertexer->FindVertices();
+ delete vertexer;
+ vertexer=0;
+
+ outfile->Close();
+ infile->Close();
+ delete infile;
+ delete outfile;
+
return;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff