-/****************************************************************************
- * This macro performs track and vertex reconstruction in TPC and ITS. *
- * The ITS Kalman tracker V2 is feeded "with" parameterized TPC tracks. *
- * *
- * Reconstruction is performed in the following steps: *
- * 1) TPC tracking parameterization *
- * 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 *
- * 5) Create a reference file with simulation info (p,PDG...) *
- * *
- * 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) *
- ****************************************************************************/
-
-// structure for track references
-typedef struct {
- Int_t lab;
- Int_t pdg;
- Int_t mumlab;
- Int_t mumpdg;
- Float_t Vx,Vy,Vz;
- Float_t Px,Py,Pz;
-} RECTRACK;
-
-//===== Functions definition =================================================
-
-void CopyVtx(const Char_t *inName,const Char_t *outName);
-
-void ITSFindClustersV2(Char_t SlowOrFast);
-
-void ITSFindTracksV2(Int_t *skipEvt);
-
-void ITSMakeRefFile(Int_t *skipEvt);
-
-void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
-
-void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
-
-void TPCParamTracks(Int_t coll,Double_t Bfield);
-
-Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
-
-void VtxFromHeader(const Char_t *outName,Bool_t smear);
-
-void VtxFromTracks(const Char_t *outName);
-
-void ZvtxFromSPD(const Char_t *outName);
-
-//=============================================================================
-
-// 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') {
-
- //---------------------------------------------------------------------
- // CONFIGURATION
+void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
//
- // _Magnetic_field_
- gBfieldValue = 0.4;
+ // Macro to create AliESDs.root using parametrized TPC tracking
+ // and AliITStrackerV2
//
- // _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_ (AliESDVertex 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";
- 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;
-
-
- // 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(collcode,gBfieldValue);
-
- // ITS clusters
- ITSFindClustersV2(SlowOrFast);
-
- // Vertex for ITS tracking
- PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
-
- break;
-
- case 'B':
- printf(" ---> only tracking in ITS <---\n");
-
- // Update list of events to be skipped
- if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
-
- break;
- }
-
- // Mark events that have to be skipped (if any)
- 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);
-
- // Tracking in ITS
- ITSFindTracksV2(skipEvt);
-
- // Vertex for analysis
- PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
-
- // Create ITS tracks reference file
- ITSMakeRefFile(skipEvt);
- delete [] skipEvt;
-
-
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return;
-}
-//-----------------------------------------------------------------------------
-void CopyVtx(const Char_t *inName,const Char_t *outName) {
-
- // Open input and output files
- TFile *inFile = new TFile(inName);
- TFile *outFile = new TFile(outName,"update");
-
- TDirectory *curdir;
- Char_t vname[20];
-
-
- for(Int_t ev=0; ev<gNevents; ev++) {
- sprintf(vname,"Vertex_%d",ev);
- AliESDVertex *vertex = (AliESDVertex*)inFile->Get(vname);
- if(!vertex) continue;
- curdir = gDirectory;
- outFile->cd();
- vertex->Write();
- curdir->cd();
- vertex = 0;
- }
-
- inFile->Close();
- outFile->Close();
- delete inFile;
- delete outFile;
-
- return;
-}
-//-----------------------------------------------------------------------------
-void ITSFindClustersV2(Char_t SlowOrFast) {
-
- printf("\n------------------------------------\n");
-
- const Char_t *name="ITSFindClustersV2";
- printf("\n %s\n",name);
- gBenchmark->Start(name);
-
- //--- taken from AliITStestV2.C--------------------------------------
+ // Input files:
+ // - galice.root
+ // - Kinematics.root
+ // - TrackReferences.root
+ // - ITS.RecPoints.root (use AliRecontruction class)
+ // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
//
- if (SlowOrFast=='f') {
- //cerr<<"Fast AliITSRecPoint(s) !\n";
- //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
- //AliITSHits2FastRecPoints();
- } else {
- gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C");
- AliITSHits2SDigits();
- gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C");
- AliITSSDigits2Digits();
- //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C");
- //AliITSDigits2RecPoints();
- }
- gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
- AliITSFindClustersV2(SlowOrFast,gNevents);
+ // A. Dainese - LNL
//
- //--------------------------------------------------------------------
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
+ /**** Initialization of the NewIO *******/
+
+ if (gAlice) {
+ delete gAlice->GetRunLoader();
+ delete gAlice;
+ gAlice=0;
+ }
+
+ AliRunLoader *rl = AliRunLoader::Open("galice.root");
+ if (rl == 0x0) {
+ cerr<<"Can not open session"<<endl;
+ return;
+ }
+ Int_t retval = rl->LoadgAlice();
+ if (retval) {
+ cerr<<"LoadgAlice returned error"<<endl;
+ delete rl;
+ return;
+ }
+ retval = rl->LoadHeader();
+ if (retval) {
+ cerr<<"LoadHeader returned error"<<endl;
+ delete rl;
+ return;
+ }
+ gAlice=rl->GetAliRun();
+
+
+ TDatabasePDG *DataBase = TDatabasePDG::Instance();
+
+ // Get field from galice.root
+ AliMagF *fiel = (AliMagF*)gAlice->Field();
+ Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+ // Set the conversion constant between curvature and Pt
+ AliTracker::SetFieldMap(fiel,kTRUE);
+
+ /**** The TPC corner ********************/
+
+ Int_t collcode = 1; // pp collisions
+ AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
+ tpcTrackerPar.Init();
+
+ //**** Switch on the PID class (mandatory!)
+ AliPID pid;
+
+ /**** The ITS corner ********************/
+
+ AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
+ if (itsl == 0x0) {
+ cerr<<"Cannot get the ITS loader"<<endl;
+ return;
+ }
+ itsl->LoadRecPoints("read");
+
+ AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
+ if (!dITS) {
+ cerr<<"Cannot find the ITS detector !"<<endl;
+ return;
+ }
+ AliITSgeom *geom = dITS->GetITSgeom();
+
+ //An instance of the ITS tracker
+ AliITStrackerV2 itsTracker(geom);
+ Int_t ITSclusters[6] = {1,1,1,1,1,1};
+ itsTracker.SetLayersNotToSkip(ITSclusters);
+
+ /***** The TREE is born *****/
+
+ TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
+ AliESD *event=0;
+ esdTree->Branch("ESD","AliESD",&event);
+
+ if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
+ if(lastEvent>rl->GetNumberOfEvents()) lastEvent=rl->GetNumberOfEvents()-1;
+ cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
+
+ TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
+ AliESDVertex *myvertex = new AliESDVertex();
+ Char_t zver[100];
+ Double_t vtx[3];
+ Double_t cvtx[6];
+
+
+ //<----------------------------------The Loop over events begins
+ TStopwatch timer;
+ Int_t trc;
+ for(Int_t i=firstEvent; i<=lastEvent; i++) {
+
+ cerr<<" Processing event number : "<<i<<endl;
+ AliESD *event = new AliESD();
+ event->SetRunNumber(gAlice->GetRunNumber());
+ event->SetEventNumber(i);
+ event->SetMagneticField(gAlice->Field()->SolenoidField());
+ rl->GetEvent(i);
+
+ //***** Primary vertex
+ sprintf(zver,"Event%d/Vertex",i);
+ myvertex = (AliESDVertex*)ppZ->Get(zver);
+ if(!myvertex) {
+ esdTree->Fill(); delete event;
+ continue;
+ }
+ event->SetVertex(myvertex);
+ myvertex->GetXYZ(vtx);
+ myvertex->GetCovMatrix(cvtx);
+
+ //***** TPC tracking
+ if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
+ printf("exiting TPC tracker with code %d in event %d\n",trc,i);
+ esdTree->Fill(); delete event;
+ continue;
+ }
+
+ //***** ITS tracking
+ itsTracker.SetVertex(vtx,cvtx);
+ TTree *itsTree=itsl->TreeR();
+ if (!itsTree) {
+ cerr<<"Can't get the ITS cluster tree !\n";
+ esdTree->Fill(); delete event;
+ return;
+ }
+ itsTracker.UnloadClusters();
+ itsTracker.LoadClusters(itsTree);
+ if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
+ printf("exiting ITS tracker with code %d in event %d\n",trc,i);
+ esdTree->Fill(); delete event;
+ continue;
+ }
+
+
+ esdTree->Fill();
+ delete event;
+
+ }//<-----------------------------------The Loop over events ends here
+ timer.Stop(); timer.Print();
+
+ // The AliESDs.root is born
+ TFile *ef = TFile::Open("AliESDs.root","RECREATE");
+ if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
+
+ esdTree->Write(); //Write the TREE and close everything
+ delete esdTree;
+ ef->Close();
+
+ delete rl;
return;
}
-//-----------------------------------------------------------------------------
-Int_t ITSFindTracksV2(Int_t *skipEvt) {
-
- printf("\n------------------------------------\n");
-
- const Char_t *name="ITSFindTracksV2";
- printf("\n %s\n",name);
- gBenchmark->Start(name);
-
-
- 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*)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<gNevents; ev++){
- // write to logfile of begun events
- fprintf(logfile,"%d\n",ev);
-
- if(skipEvt[ev]) continue;
- printf(" --- Processing event %d ---\n",ev);
-
- // pass event number to the tracker
- tracker.SetEventNumber(ev);
-
- // set position of primary vertex
- sprintf(vname,"Vertex_%d",ev);
- AliESDVertex *vertex = (AliESDVertex*)inVertex->Get(vname);
- if(vertex) {
- vertex->GetXYZ(vtx);
- delete vertex;
- } else {
- printf(" AliESDVertex 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
-
- // no vtx constraint if vertex not found
- if(vtx[2]<-999.) {
- flag1stPass=0;
- vtx[2]=0.;
- }
-
- tracker.SetVertex(vtx);
-
- // setup vertex constraint in the two tracking passes
- Int_t flags[2];
- flags[0]=flag1stPass;
- tracker.SetupFirstPass(flags);
- flags[0]=flag2ndPass;
- tracker.SetupSecondPass(flags);
-
- // find the tracks
- tracker.Clusters2Tracks(inTPCtrks,outFile);
-
- } // loop on events
-
- fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
- fclose(logfile);
-
- delete geom;
-
- inTPCtrks->Close();
- inClusters->Close();
- inVertex->Close();
- outFile->Close();
-
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return;
-}
-//-----------------------------------------------------------------------------
-void ITSMakeRefFile(Int_t *skipEvt) {
-
- printf("\n------------------------------------\n");
-
- const Char_t *name="ITSMakeRefFile";
- printf("\n %s\n",name);
- gBenchmark->Start(name);
-
-
- TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
- TFile *trk = TFile::Open("AliITStracksV2.root");
- TFile *kin = TFile::Open("galice.root");
-
-
- // Get gAlice object from file
- gAlice=(AliRun*)kin->Get("gAlice");
-
- Int_t label;
- TParticle *Part;
- TParticle *Mum;
- RECTRACK rectrk;
-
-
- for(Int_t ev=0; ev<gNevents; ev++){
- if(skipEvt[ev]) continue;
- printf(" --- Processing event %d ---\n",ev);
-
- gAlice->GetEvent(ev);
-
- trk->cd();
-
- // Tree with ITS tracks
- char tname[100];
- sprintf(tname,"TreeT_ITS_%d",ev);
-
- TTree *tracktree=(TTree*)trk->Get(tname);
- if(!tracktree) continue;
- AliITStrackV2 *itstrack=new AliITStrackV2;
- tracktree->SetBranchAddress("tracks",&itstrack);
- Int_t nentr=(Int_t)tracktree->GetEntries();
-
- // Tree for true track parameters
- char ttname[100];
- sprintf(ttname,"Tree_Ref_%d",ev);
- TTree *reftree = new TTree(ttname,"Tree with true track params");
- reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
-
- for(Int_t i=0; i<nentr; i++) {
- tracktree->GetEvent(i);
- label = TMath::Abs(itstrack->GetLabel());
-
- Part = (TParticle*)gAlice->Particle(label);
- rectrk.lab=label;
- rectrk.pdg=Part->GetPdgCode();
- rectrk.mumlab = Part->GetFirstMother();
- if(Part->GetFirstMother()>=0) {
- Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
- rectrk.mumpdg=Mum->GetPdgCode();
- } else {
- rectrk.mumpdg=-1;
- }
- rectrk.Vx=Part->Vx();
- rectrk.Vy=Part->Vy();
- rectrk.Vz=Part->Vz();
- rectrk.Px=Part->Px();
- rectrk.Py=Part->Py();
- rectrk.Pz=Part->Pz();
-
- reftree->Fill();
- } // loop on tracks
-
- out->cd();
- reftree->Write();
-
- delete itstrack;
- delete reftree;
- } // loop on events
-
- trk->Close();
- kin->Close();
- out->Close();
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
-
- return;
-}
-//-----------------------------------------------------------------------------
-void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
-
- printf("\n------------------------------------\n");
- printf("\nChecking for events to skip...\n");
-
- Int_t evt,ncol;
-
- 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);
-
- return;
-}
-//-----------------------------------------------------------------------------
-void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
-
- printf("\n------------------------------------\n");
-
- 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;
- }
-
- gBenchmark->Stop(name);
- gBenchmark->Show(name);
-
- return;
-}
-//-----------------------------------------------------------------------------
-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");
- }
-
- return 1;
-}
-//-----------------------------------------------------------------------------
-void VtxFromHeader(const Char_t *outName,Bool_t smear) {
-
- TDatime t;
- UInt_t seed = t.Get();
- gRandom->SetSeed(seed);
-
- TFile *galice = new TFile("galice.root");
- TFile *outFile = new TFile(outName,"update");
-
- 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];
-
- 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 AliESDVertex
- AliESDVertex *vertex = new AliESDVertex(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 AliESDVertex
- AliESDVertex *vertex = new AliESDVertex(pos,sigma,vname);
- }
- delete header;
- // write AliESDVertex 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;
- }
-
- 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;
-}
-//-----------------------------------------------------------------------------
-
-
-
-
-
-
-
-
-
* - All the code necessary to create a BataBase has been included in *
* class (see the macro AliTPCtrackingParamDB.C for the details). *
* *
- * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
+ * 2006/03/16: Adapted to ESD input/output *
* *
+ * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it *
+ * adapted to ESD output by Marcello Lunardon, Padova *
**************************************************************************/
// *
// This is a dummy comment
#include <TH2F.h>
#include <TLegend.h>
#include <TLine.h>
+#include <TMatrixD.h>
#include <TParticle.h>
#include <TStyle.h>
#include <TSystem.h>
//-----------------------------------------------------------------------------
AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
- Int_t kn, const char* evfoldname):
+ const char* evfoldname):
fEvFolderName(evfoldname) {
//-----------------------------------------------------------------------------
// This is the class conctructor
//-----------------------------------------------------------------------------
- fNevents = kn; // events to be processed
fBz = kBz; // value of the z component of L3 field (Tesla)
fColl = kcoll; // collision code (0: PbPb6000; 1: pp)
fSelAndSmear = kTRUE; // by default selection and smearing are done
- if(fBz!=0.4) {
- cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n";
- cerr<<" Available: 0.4"<<endl;
+ if(fBz!=0.4 && fBz!=0.5) {
+ Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n Available: 0.4 or 0.5");
}
- if(fColl!=0 && fColl!=1) {
- cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n";
- cerr<<" Available: 0 -> PbPb6000"<<endl;
- cerr<<" 1 -> pp"<<endl;
+ if(fColl!=0 && fColl!=1) {
+ Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n Available: 0 -> PbPb6000\n 1 -> pp");
}
fDBfileName = gSystem->Getenv("ALICE_ROOT");
if(fColl==0) fDBfileName.Append("PbPb6000");
if(fColl==1) fDBfileName.Append("pp");
if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
+ // use same DB for 0.4 and 0.5 T; for 0.5 T, correction done in CookTrack()
+ if(fBz==0.5) fDBfileName.Append("_B0.4T.root");
}
//-----------------------------------------------------------------------------
AliTPCtrackerParam::~AliTPCtrackerParam() {}
fAlpha *= TMath::Pi();
}
//-----------------------------------------------------------------------------
-Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out) {
+Int_t AliTPCtrackerParam::Init() {
//-----------------------------------------------------------------------------
-// This function creates the TPC parameterized tracks
+// This function reads the DB from file
//-----------------------------------------------------------------------------
- Error("BuildTPCtracks","in and out parameters ignored. new io");
-
-/********************************************/
- AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
- if (rl == 0x0)
- {
- Error("BuildTPCtracks","Can not get Run Loader from event folder named %s.",
- fEvFolderName.Data());
- return 2;
- }
- AliLoader* tpcloader = rl->GetLoader("TPCLoader");
- if (tpcloader == 0x0)
- {
- Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
- return 3;
- }
-
-/********************************************/
-
- TFile *fileDB=0;
- TTree *covTreePi[50];
- TTree *covTreeKa[50];
- TTree *covTreePr[50];
- TTree *covTreeEl[50];
- TTree *covTreeMu[50];
-
if(fSelAndSmear) {
- cerr<<"+++\n+++ Reading DataBase from:\n+++ "<<
- fDBfileName.Data()<<"\n+++\n";
+ printf("+++\n+++ Reading DataBase from:%s\n+++\n+++\n",fDBfileName.Data());
// Read paramters from DB file
if(!ReadAllData(fDBfileName.Data())) {
- cerr<<"AliTPCtrackerParam::BuildTPCtracks: \
- Could not read data from DB\n\n"; return 1;
+ printf("AliTPCtrackerParam::BuildTPCtracks: \
+ Could not read data from DB\n\n"); return 1;
}
+
+ } else printf("\n ! Creating ALL TRUE tracks at TPC inner radius !\n\n");
- // Read the trees with regularized cov. matrices from DB
- TString str;
- fileDB = TFile::Open(fDBfileName.Data());
- Int_t nBinsPi = fDBgridPi.GetTotBins();
- for(Int_t l=0; l<nBinsPi; l++) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- str+=l;
- covTreePi[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsKa = fDBgridKa.GetTotBins();
- for(Int_t l=0; l<nBinsKa; l++) {
- str = "/CovMatrices/Kaons/CovTreeKa_bin";
- str+=l;
- covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsPr = fDBgridPr.GetTotBins();
- for(Int_t l=0; l<nBinsPr; l++) {
- if(fColl==0) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- } else {
- str = "/CovMatrices/Protons/CovTreePr_bin";
- }
- str+=l;
- covTreePr[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsEl = fDBgridEl.GetTotBins();
- for(Int_t l=0; l<nBinsEl; l++) {
- str = "/CovMatrices/Electrons/CovTreeEl_bin";
- str+=l;
- covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsMu = fDBgridMu.GetTotBins();
- for(Int_t l=0; l<nBinsMu; l++) {
- if(fColl==0) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- } else {
- str = "/CovMatrices/Muons/CovTreeMu_bin";
- }
- str+=l;
- covTreeMu[l] = (TTree*)fileDB->Get(str.Data());
- }
- } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
+ // Check if value in the galice file is equal to selected one (fBz)
+ AliMagF *fiel = (AliMagF*)gAlice->Field();
+ Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+ printf("Magnetic field is %6.2f Tesla\n",fieval);
+ if(fBz!=fieval) {
+ printf("AliTPCtrackerParam::BuildTPCtracks: Invalid field!");
+ printf("Field selected is: %f T\n",fBz);
+ printf("Field found on file is: %f T\n",fieval);
+ return 1;
+ }
+
+ // Set the conversion constant between curvature and Pt
+ AliTracker::SetFieldMap(fiel,kTRUE);
+
+ return 0;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::BuildTPCtracks(AliESD *event) {
+//-----------------------------------------------------------------------------
+// This function creates the TPC parameterized tracks and writes them
+// as AliESDtrack objects in the ESD event
+//-----------------------------------------------------------------------------
- TFile *infile=(TFile*)inp;
- infile->cd();
- // Get gAlice object from file
+ if(!event) return -1;
+
+ AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName);
+ if (rl == 0x0) {
+ Error("BuildTPCtracks","Can not get Run Loader from event folder named %s",
+ fEvFolderName.Data());
+ return 2;
+ }
+ AliLoader* tpcloader = rl->GetLoader("TPCLoader");
+ if (tpcloader == 0x0) {
+ Error("BuildTPCtracks","Can not get TPC Loader from Run Loader.");
+ return 3;
+ }
- rl->LoadgAlice();
- rl->LoadHeader();
+ // Get gAlice object from file
+ if(!gAlice) rl->LoadgAlice();
+ //rl->LoadHeader();
+ rl->LoadKinematics();
tpcloader->LoadHits("read");
if(!(gAlice=rl->GetAliRun())) {
- cerr<<"Can not get gAlice from Run Loader !\n";
+ printf("Cannot get gAlice from Run Loader !\n");
return 1;
}
- // Check if value in the galice file is equal to selected one (fBz)
- AliMagF *fiel = (AliMagF*)gAlice->Field();
- Double_t fieval=(Double_t)fiel->SolenoidField()/10.;
- printf("Magnetic field is %6.2f Tesla\n",fieval);
- if(fBz!=fieval) {
- cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!"<<endl;
- cerr<<"Field selected is: "<<fBz<<" T\n";
- cerr<<"Field found on file is: "<<fieval<<" T\n";
- return 0;
- }
-
// Get TPC detector
AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
- Int_t ver = tpc->IsVersion();
- cerr<<"+++ TPC version "<<ver<<" has been found !\n";
rl->CdGAFile();
- AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
+ AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60");
if(digp){
delete digp;
digp = new AliTPCParamSR();
}
- else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60");
+ else digp=(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60");
if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; }
tpc->SetParam(digp);
- // Set the conversion constant between curvature and Pt
- AliTracker::SetFieldMap(fiel,kTRUE);
-
TParticle *part=0;
AliTPCseedGeant *seed=0;
AliTPCtrack *tpctrack=0;
Double_t sPt,sEta;
- Int_t cl=0,bin,label,pdg,charge;
+ Int_t bin,label,pdg,charge;
Int_t tracks;
Int_t nParticles,nSeeds,arrentr;
- Char_t tname[100];
//Int_t nSel=0,nAcc=0;
+ Int_t evt=event->GetEventNumber();
+
+ tracks=0;
- // loop over first n events in file
- for(Int_t evt=0; evt<fNevents; evt++){
- cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
- tracks=0;
-
- // tree for TPC tracks
- sprintf(tname,"TreeT_TPC_%d",evt);
- TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
- tracktree->Branch("tracks","AliTPCtrack",&tpctrack,20000,0);
-
- // array for TPC tracks
- TObjArray tArray(20000);
-
- // array for TPC seeds with geant info
- TObjArray sArray(20000);
-
- // get the particles stack
- nParticles = (Int_t)gAlice->GetEvent(evt);
+ // array for TPC tracks
+ TObjArray tArray(20000);
+
+ // array for TPC seeds with geant info
+ TObjArray sArray(20000);
+
+ // get the particles stack
+ nParticles = (Int_t)gAlice->GetEvent(evt);
- Bool_t *done = new Bool_t[nParticles];
- Int_t *pdgCodes = new Int_t[nParticles];
- Double_t *ptkine = new Double_t[nParticles];
- Double_t *pzkine = new Double_t[nParticles];
-
- // loop on particles and store pdg codes
- for(Int_t l=0; l<nParticles; l++) {
- part = (TParticle*)gAlice->GetMCApp()->Particle(l);
- pdgCodes[l] = part->GetPdgCode();
- ptkine[l] = part->Pt();
- pzkine[l] = part->Pz();
- done[l] = kFALSE;
- }
- cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nParticles<<
- "\n+++\n";
-
- cerr<<"\n ********** MAKING SEEDS *****************"<<endl<<endl;
-
-
- // Create the seeds for the TPC tracks at the inner radius of TPC
- if(fColl==0) {
- // Get TreeH with hits
- TTree *th = tpcloader->TreeH();
- MakeSeedsFromHits(tpc,th,sArray);
- } else {
- // Get TreeTR with track references
- TTree *ttr = rl->TreeTR();
- MakeSeedsFromRefs(ttr,sArray);
- }
-
-
- nSeeds = sArray.GetEntries();
- cerr<<"\n\n+++\n+++ Number of seeds: "<<nSeeds<<"\n+++\n";
-
-
- cerr<<"\n ********** BUILDING TRACKS **************"<<endl<<endl;
-
- // loop over entries in sArray
- for(Int_t l=0; l<nSeeds; l++) {
- //if(l%1000==0) cerr<<" --- Processing seed "
- // <<l<<" of "<<nSeeds<<" ---\r";
-
- seed = (AliTPCseedGeant*)sArray.At(l);
-
- // this is TEMPORARY: only for reconstruction of pp production for charm
- if(fColl==1) cl = CheckLabel(seed,nParticles,ptkine,pzkine);
- if(cl) continue;
+ Bool_t *done = new Bool_t[nParticles];
+ Int_t *pdgCodes = new Int_t[nParticles];
+
+ // loop on particles and store pdg codes
+ for(Int_t l=0; l<nParticles; l++) {
+ part = (TParticle*)gAlice->GetMCApp()->Particle(l);
+ pdgCodes[l] = part->GetPdgCode();
+ done[l] = kFALSE;
+ }
- // Get track label
- label = seed->GetLabel();
-
- // check if this track has already been processed
- if(done[label]) continue;
- // PDG code & electric charge
- pdg = pdgCodes[label];
- if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
- else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
- else continue;
- pdg = TMath::Abs(pdg);
- if(pdg>3000) pdg=211;
-
- if(fSelAndSmear) SetParticle(pdg);
+ printf("+++ Number of particles: %d\n",nParticles);
+ // Create the seeds for the TPC tracks at the inner radius of TPC
+ if(fColl==0) {
+ // Get TreeH with hits
+ TTree *th = tpcloader->TreeH();
+ MakeSeedsFromHits(tpc,th,sArray);
+ } else {
+ // Get TreeTR with track references
+ rl->LoadTrackRefs();
+ TTree *ttr = rl->TreeTR();
+ if (!ttr) return -3;
+ MakeSeedsFromRefs(ttr,sArray);
+ }
- sPt = seed->GetPt();
- sEta = seed->GetEta();
-
- // Apply selection according to TPC efficiency
- //if(TMath::Abs(pdg)==211) nAcc++;
- if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue;
- //if(TMath::Abs(pdg)==211) nSel++;
-
- // create AliTPCtrack object
- BuildTrack(seed,charge);
-
- if(fSelAndSmear) {
- bin = fDBgrid->GetBin(sPt,sEta);
- switch (pdg) {
- case 211:
- fCovTree = covTreePi[bin];
- break;
- case 321:
- fCovTree = covTreeKa[bin];
- break;
- case 2212:
- fCovTree = covTreePr[bin];
- break;
- case 11:
- fCovTree = covTreeEl[bin];
- break;
- case 13:
- fCovTree = covTreeMu[bin];
- break;
- }
- // deal with covariance matrix and smearing of parameters
- CookTrack(sPt,sEta);
-
- // assign the track a dE/dx and make a rough PID
- CookdEdx(sPt,sEta);
+ nSeeds = sArray.GetEntries();
+ printf("+++ Number of seeds: %d\n",nSeeds);
+
+ // loop over entries in sArray
+ for(Int_t l=0; l<nSeeds; l++) {
+ //if(l%1==0) printf(" --- Processing seed %d of %d ---\n",l,nSeeds);
+
+ seed = (AliTPCseedGeant*)sArray.At(l);
+
+ // Get track label
+ label = seed->GetLabel();
+
+ // check if this track has already been processed
+ if(done[label]) continue;
+
+ // PDG code & electric charge
+ pdg = pdgCodes[label];
+ if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; }
+ else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; }
+ else continue;
+ pdg = TMath::Abs(pdg);
+ if(pdg>3000) pdg=211;
+
+ if(fSelAndSmear) SetParticle(pdg);
+
+ sPt = seed->GetPt();
+ sEta = seed->GetEta();
+
+ // Apply selection according to TPC efficiency
+ //if(TMath::Abs(pdg)==211) nAcc++;
+ if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue;
+ //if(TMath::Abs(pdg)==211) nSel++;
+
+ // create AliTPCtrack object
+ BuildTrack(seed,charge);
+
+ if(fSelAndSmear) {
+ bin = fDBgrid->GetBin(sPt,sEta);
+ switch (pdg) {
+ case 211:
+ //fCovTree = &(fCovTreePi[bin]);
+ fCovTree = fCovTreePi[bin];
+ break;
+ case 321:
+ //fCovTree = &(fCovTreeKa[bin]);
+ fCovTree = fCovTreeKa[bin];
+ break;
+ case 2212:
+ //fCovTree = &(fCovTreePr[bin]);
+ fCovTree = fCovTreePr[bin];
+ break;
+ case 11:
+ //fCovTree = &(fCovTreeEl[bin]);
+ fCovTree = fCovTreeEl[bin];
+ break;
+ case 13:
+ //fCovTree = &(fCovTreeMu[bin]);
+ fCovTree = fCovTreeMu[bin];
+ break;
}
+ // deal with covariance matrix and smearing of parameters
+ CookTrack(sPt,sEta);
- // put track in array
- AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
- iotrack->SetLabel(label);
- tArray.AddLast(iotrack);
- // Mark track as "done" and register the pdg code
- done[label] = kTRUE;
- tracks++;
-
- } // loop over entries in sArray
-
-
- // sort array with TPC tracks (decreasing pT)
- tArray.Sort();
-
- arrentr = tArray.GetEntriesFast();
- for(Int_t l=0; l<arrentr; l++) {
- tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
- tracktree->Fill();
+ // assign the track a dE/dx and make a rough PID
+ CookdEdx(sPt,sEta);
}
-
-
- // write the tree with tracks in the output file
- out->cd();
- tracktree->Write();
- delete tracktree;
- delete [] done;
- delete [] pdgCodes;
- delete [] ptkine;
- delete [] pzkine;
-
- printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
- //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
-
- sArray.Delete();
- tArray.Delete();
+ // put track in array
+ AliTPCtrack *iotrack = new AliTPCtrack(fTrack);
+ iotrack->SetLabel(label);
+ tArray.AddLast(iotrack);
+ // Mark track as "done" and register the pdg code
+ done[label] = kTRUE;
+ tracks++;
+
+ } // loop over entries in sArray
+
+ // sort array with TPC tracks (decreasing pT)
+ tArray.Sort();
+
+ // convert to AliESDtrack and write to AliESD
+ arrentr = tArray.GetEntriesFast();
+ Int_t idx;
+ Double_t wgts[5];
+ for(Int_t l=0; l<arrentr; l++) {
+ tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
+ AliESDtrack ioESDtrack;
+ ioESDtrack.UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
+ // rough PID
+ wgts[0]=0.; wgts[1]=0.; wgts[2]=0.; wgts[3]=0.; wgts[4]=0.;
+ if(TMath::Abs(tpctrack->GetMass()-0.9)<0.1) {
+ idx = 4; // proton
+ } else if(TMath::Abs(tpctrack->GetMass()-0.5)<0.1) {
+ idx = 3; // kaon
+ } else {
+ idx = 2; // pion
+ }
+ wgts[idx] = 1.;
+ ioESDtrack.SetESDpid(wgts);
+ event->AddTrack(&ioESDtrack);
+ }
+
+
+ delete [] done;
+ delete [] pdgCodes;
+
+ printf("+++ Number of TPC tracks: %d\n",tracks);
+ //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
+
+ sArray.Delete();
+ tArray.Delete();
- infile->cd();
- } // loop on events
-
- if(fileDB) fileDB->Close();
-
return 0;
}
//-----------------------------------------------------------------------------
cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.;
// Magnetic field
- TVector3 bfield(0.,0.,fBz);
+ TVector3 bfield(0.,0.,-fBz);
// radius [cm] of track projection in (x,y)
- Double_t rho = s->GetPt()*100./0.299792458/bfield.Z();
+ Double_t rho = s->GetPt()*100./0.299792458/TMath::Abs(bfield.Z());
// center of track projection in local reference frame
TVector3 sMom,sPos;
xx[0] = s->GetYL();
xx[1] = s->GetZL();
xx[3] = s->GetPz()/s->GetPt();
- xx[4] = -ch/rho;
+ xx[4] = ch/rho;
xx[2] = xx[4]*x0;
// create the object AliTPCtrack
return;
}
//-----------------------------------------------------------------------------
-Int_t AliTPCtrackerParam::CheckLabel(AliTPCseedGeant *s,Int_t nPart,
- Double_t *ptkine,Double_t *pzkine) const {
-//-----------------------------------------------------------------------------
-// This function checks if the label of the seed has been correctly
-// assigned (to do only for pp charm production with AliRoot v3-08-02)
-//-----------------------------------------------------------------------------
-
- Int_t sLabel = s->GetLabel();
- Double_t sPt = s->GetPt();
- Double_t sPz = s->GetPz();
-
- // check if the label is correct (comparing momentum)
- if(sLabel<nPart &&
- TMath::Abs(sPt-ptkine[sLabel])*
- TMath::Abs(sPz-pzkine[sLabel])<0.001) return 0;
-
- if((sLabel-30)>=nPart) return 1;
-
- Double_t diff=0,mindiff=1000.;
- Int_t bestLabel=0;
-
- for(Int_t i=sLabel-30; i<sLabel; i++) {
- if(i<0 || i>=nPart) continue;
- diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]);
- if(diff<mindiff) { mindiff = diff; bestLabel = i; }
- }
-
- if(mindiff>0.001) return 1;
- s->SetLabel(bestLabel);
-
- return 0;
-}
-//-----------------------------------------------------------------------------
void AliTPCtrackerParam::CompareTPCtracks(
+ Int_t nEvents,
const Char_t* galiceName,
const Char_t* trkGeaName,
const Char_t* trkKalName,
Char_t tname[100];
// loop on events in file
- for(Int_t evt=0; evt<fNevents; evt++) {
+ for(Int_t evt=0; evt<nEvents; evt++) {
cerr<<"\n --- Reading tracks for event "<<evt<<" ---\n\n";
sprintf(tname,"TreeT_TPC_%d",evt);
//Very rough PID
Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt;
+ Double_t massPi = (Double_t)TDatabasePDG::Instance()->GetParticle(211)->Mass();
+ Double_t massKa = (Double_t)TDatabasePDG::Instance()->GetParticle(321)->Mass();
+ Double_t massPr = (Double_t)TDatabasePDG::Instance()->GetParticle(2212)->Mass();
+
if (p<0.6) {
if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
- t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
}
if (dEdx < 39.+ 12./p/p) {
- t.AssignMass(AliPID::ParticleMass(AliPID::kKaon)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massKa); new(&fTrack) AliTPCtrack(t); return;
}
- t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
}
if (p<1.2) {
if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) {
- t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
}
- t.AssignMass(AliPID::ParticleMass(AliPID::kProton)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return;
}
- t.AssignMass(AliPID::ParticleMass(AliPID::kPion)); new(&fTrack) AliTPCtrack(t); return;
+ t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return;
}
//-----------------------------------------------------------------------------
void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) {
cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
cc[13]= covmat.c43;
for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
- cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
-
- TMatrixD covMatSmear(5,5);
-
+ cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
+
+
+ TMatrixD covMatSmear(5,5);
covMatSmear = GetSmearingMatrix(cc,pt,eta);
// get track original parameters
xx[4]=fTrack.GetC();
// use smearing matrix to smear the original parameters
+ xxsm[0]=xref;
SmearTrack(xx,xxsm,covMatSmear);
AliTPCtrack track(0,xxsm,cc,xref,alpha);
Int_t label;
Int_t nTracks=(Int_t)th->GetEntries();
- cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
- "\n+++\n\n";
+ cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
+ "\n";
AliTPChit *tpcHit=0;
if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
b->SetAddress(&tkRefArray);
Int_t nTkRef = (Int_t)b->GetEntries();
- cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<<nTkRef<<
- "\n+++\n\n";
+ cerr<<"+++ Number of entries in TreeTR(TPC): "<<nTkRef<<"\n";
// loop on track references
for(Int_t l=0; l<nTkRef; l++){
- if(l%1000==0) cerr<<" --- Processing primary track "
- <<l<<" of "<<nTkRef<<" ---\r";
+ //if(l%1000==0) cerr<<" --- Processing primary track "
+ // <<l<<" of "<<nTkRef<<" ---\r";
if(!b->GetEvent(l)) continue;
nnn = tkRefArray->GetEntriesFast();
AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
// reject if not at the inner part of TPC
- if(TMath::Abs(ioseed->GetXL()-82.9701) > 0.1) {
+ if(TMath::Abs(ioseed->GetXL()-83.8) > 0.2) {
+ //cerr<<"not at TPC inner part: XL = "<<ioseed->GetXL()<<endl;
delete ioseed; continue;
}
if(!ReaddEdx(inName,11)) return 0;
if(!ReaddEdx(inName,13)) return 0;
if(!ReadDBgrid(inName)) return 0;
+ if(!ReadCovTrees(inName)) return 0;
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadCovTrees(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the covariance matrix trees from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
+ return 0; }
+
+ TString str;
+
+ TFile *inFile = TFile::Open(inName);
+
+
+ Int_t nBinsPi = fDBgridPi.GetTotBins();
+ for(Int_t l=0; l<nBinsPi; l++) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ str+=l;
+ fCovTreePi[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreePi[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsKa = fDBgridKa.GetTotBins();
+ for(Int_t l=0; l<nBinsKa; l++) {
+ str = "/CovMatrices/Kaons/CovTreeKa_bin";
+ str+=l;
+ fCovTreeKa[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeKa[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsPr = fDBgridPr.GetTotBins();
+ for(Int_t l=0; l<nBinsPr; l++) {
+ if(fColl==0) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ } else {
+ str = "/CovMatrices/Protons/CovTreePr_bin";
+ }
+ str+=l;
+ fCovTreePr[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreePr[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsEl = fDBgridEl.GetTotBins();
+ for(Int_t l=0; l<nBinsEl; l++) {
+ str = "/CovMatrices/Electrons/CovTreeEl_bin";
+ str+=l;
+ fCovTreeEl[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeEl[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsMu = fDBgridMu.GetTotBins();
+ for(Int_t l=0; l<nBinsMu; l++) {
+ if(fColl==0) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ } else {
+ str = "/CovMatrices/Muons/CovTreeMu_bin";
+ }
+ str+=l;
+ fCovTreeMu[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeMu[l] = new TTree(*fCovTree);
+ }
+
+ //inFile->Close();
return 1;
}
if(gSystem->AccessPathName(inName,kFileExists)) {
cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
return 0; }
+ // Introduced change to "reverse" the TMatrixD read from file.
+ // Needed because storage mode of TMatrixD changed from column-wise
+ // to rwo-wise in ROOT.
+ //
+ // A.Dainese 03/06/05
+
+ TMatrixD dummyMatrix(9,3);
TFile *inFile = TFile::Open(inName);
switch (pdg) {
case 211:
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
new(&fRegParPi) TMatrixD(*fRegPar);
+ //printf("before\n");
+ //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
+ dummyMatrix = fRegParPi;
+ fRegParPi(0,0) = dummyMatrix(0,0);
+ fRegParPi(0,1) = dummyMatrix(0,1);
+ fRegParPi(0,2) = dummyMatrix(0,2);
+ fRegParPi(1,0) = dummyMatrix(3,0);
+ fRegParPi(1,1) = dummyMatrix(1,1);
+ fRegParPi(1,2) = dummyMatrix(1,2);
+ fRegParPi(2,0) = dummyMatrix(6,0);
+ fRegParPi(2,1) = dummyMatrix(3,2);
+ fRegParPi(2,2) = dummyMatrix(2,2);
+ fRegParPi(3,0) = dummyMatrix(1,0);
+ fRegParPi(3,1) = dummyMatrix(4,0);
+ fRegParPi(3,2) = dummyMatrix(7,0);
+ fRegParPi(4,0) = dummyMatrix(3,1);
+ fRegParPi(4,1) = dummyMatrix(4,1);
+ fRegParPi(4,2) = dummyMatrix(7,1);
+ fRegParPi(5,0) = dummyMatrix(6,1);
+ fRegParPi(5,1) = dummyMatrix(4,2);
+ fRegParPi(5,2) = dummyMatrix(7,2);
+ fRegParPi(6,0) = dummyMatrix(2,0);
+ fRegParPi(6,1) = dummyMatrix(5,0);
+ fRegParPi(6,2) = dummyMatrix(8,0);
+ fRegParPi(7,0) = dummyMatrix(2,1);
+ fRegParPi(7,1) = dummyMatrix(5,1);
+ fRegParPi(7,2) = dummyMatrix(8,1);
+ fRegParPi(8,0) = dummyMatrix(6,2);
+ fRegParPi(8,1) = dummyMatrix(5,2);
+ fRegParPi(8,2) = dummyMatrix(8,2);
+ //printf("after\n");
+ //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2));
break;
case 321:
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
new(&fRegParKa) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParKa;
+ fRegParKa(0,0) = dummyMatrix(0,0);
+ fRegParKa(0,1) = dummyMatrix(0,1);
+ fRegParKa(0,2) = dummyMatrix(0,2);
+ fRegParKa(1,0) = dummyMatrix(3,0);
+ fRegParKa(1,1) = dummyMatrix(1,1);
+ fRegParKa(1,2) = dummyMatrix(1,2);
+ fRegParKa(2,0) = dummyMatrix(6,0);
+ fRegParKa(2,1) = dummyMatrix(3,2);
+ fRegParKa(2,2) = dummyMatrix(2,2);
+ fRegParKa(3,0) = dummyMatrix(1,0);
+ fRegParKa(3,1) = dummyMatrix(4,0);
+ fRegParKa(3,2) = dummyMatrix(7,0);
+ fRegParKa(4,0) = dummyMatrix(3,1);
+ fRegParKa(4,1) = dummyMatrix(4,1);
+ fRegParKa(4,2) = dummyMatrix(7,1);
+ fRegParKa(5,0) = dummyMatrix(6,1);
+ fRegParKa(5,1) = dummyMatrix(4,2);
+ fRegParKa(5,2) = dummyMatrix(7,2);
+ fRegParKa(6,0) = dummyMatrix(2,0);
+ fRegParKa(6,1) = dummyMatrix(5,0);
+ fRegParKa(6,2) = dummyMatrix(8,0);
+ fRegParKa(7,0) = dummyMatrix(2,1);
+ fRegParKa(7,1) = dummyMatrix(5,1);
+ fRegParKa(7,2) = dummyMatrix(8,1);
+ fRegParKa(8,0) = dummyMatrix(6,2);
+ fRegParKa(8,1) = dummyMatrix(5,2);
+ fRegParKa(8,2) = dummyMatrix(8,2);
break;
case 2212:
if(fColl==0) {
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
}
new(&fRegParPr) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParPr;
+ fRegParPr(0,0) = dummyMatrix(0,0);
+ fRegParPr(0,1) = dummyMatrix(0,1);
+ fRegParPr(0,2) = dummyMatrix(0,2);
+ fRegParPr(1,0) = dummyMatrix(3,0);
+ fRegParPr(1,1) = dummyMatrix(1,1);
+ fRegParPr(1,2) = dummyMatrix(1,2);
+ fRegParPr(2,0) = dummyMatrix(6,0);
+ fRegParPr(2,1) = dummyMatrix(3,2);
+ fRegParPr(2,2) = dummyMatrix(2,2);
+ fRegParPr(3,0) = dummyMatrix(1,0);
+ fRegParPr(3,1) = dummyMatrix(4,0);
+ fRegParPr(3,2) = dummyMatrix(7,0);
+ fRegParPr(4,0) = dummyMatrix(3,1);
+ fRegParPr(4,1) = dummyMatrix(4,1);
+ fRegParPr(4,2) = dummyMatrix(7,1);
+ fRegParPr(5,0) = dummyMatrix(6,1);
+ fRegParPr(5,1) = dummyMatrix(4,2);
+ fRegParPr(5,2) = dummyMatrix(7,2);
+ fRegParPr(6,0) = dummyMatrix(2,0);
+ fRegParPr(6,1) = dummyMatrix(5,0);
+ fRegParPr(6,2) = dummyMatrix(8,0);
+ fRegParPr(7,0) = dummyMatrix(2,1);
+ fRegParPr(7,1) = dummyMatrix(5,1);
+ fRegParPr(7,2) = dummyMatrix(8,1);
+ fRegParPr(8,0) = dummyMatrix(6,2);
+ fRegParPr(8,1) = dummyMatrix(5,2);
+ fRegParPr(8,2) = dummyMatrix(8,2);
break;
case 11:
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
new(&fRegParEl) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParEl;
+ fRegParEl(0,0) = dummyMatrix(0,0);
+ fRegParEl(0,1) = dummyMatrix(0,1);
+ fRegParEl(0,2) = dummyMatrix(0,2);
+ fRegParEl(1,0) = dummyMatrix(3,0);
+ fRegParEl(1,1) = dummyMatrix(1,1);
+ fRegParEl(1,2) = dummyMatrix(1,2);
+ fRegParEl(2,0) = dummyMatrix(6,0);
+ fRegParEl(2,1) = dummyMatrix(3,2);
+ fRegParEl(2,2) = dummyMatrix(2,2);
+ fRegParEl(3,0) = dummyMatrix(1,0);
+ fRegParEl(3,1) = dummyMatrix(4,0);
+ fRegParEl(3,2) = dummyMatrix(7,0);
+ fRegParEl(4,0) = dummyMatrix(3,1);
+ fRegParEl(4,1) = dummyMatrix(4,1);
+ fRegParEl(4,2) = dummyMatrix(7,1);
+ fRegParEl(5,0) = dummyMatrix(6,1);
+ fRegParEl(5,1) = dummyMatrix(4,2);
+ fRegParEl(5,2) = dummyMatrix(7,2);
+ fRegParEl(6,0) = dummyMatrix(2,0);
+ fRegParEl(6,1) = dummyMatrix(5,0);
+ fRegParEl(6,2) = dummyMatrix(8,0);
+ fRegParEl(7,0) = dummyMatrix(2,1);
+ fRegParEl(7,1) = dummyMatrix(5,1);
+ fRegParEl(7,2) = dummyMatrix(8,1);
+ fRegParEl(8,0) = dummyMatrix(6,2);
+ fRegParEl(8,1) = dummyMatrix(5,2);
+ fRegParEl(8,2) = dummyMatrix(8,2);
break;
case 13:
if(fColl==0) {
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
}
new(&fRegParMu) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParMu;
+ fRegParMu(0,0) = dummyMatrix(0,0);
+ fRegParMu(0,1) = dummyMatrix(0,1);
+ fRegParMu(0,2) = dummyMatrix(0,2);
+ fRegParMu(1,0) = dummyMatrix(3,0);
+ fRegParMu(1,1) = dummyMatrix(1,1);
+ fRegParMu(1,2) = dummyMatrix(1,2);
+ fRegParMu(2,0) = dummyMatrix(6,0);
+ fRegParMu(2,1) = dummyMatrix(3,2);
+ fRegParMu(2,2) = dummyMatrix(2,2);
+ fRegParMu(3,0) = dummyMatrix(1,0);
+ fRegParMu(3,1) = dummyMatrix(4,0);
+ fRegParMu(3,2) = dummyMatrix(7,0);
+ fRegParMu(4,0) = dummyMatrix(3,1);
+ fRegParMu(4,1) = dummyMatrix(4,1);
+ fRegParMu(4,2) = dummyMatrix(7,1);
+ fRegParMu(5,0) = dummyMatrix(6,1);
+ fRegParMu(5,1) = dummyMatrix(4,2);
+ fRegParMu(5,2) = dummyMatrix(7,2);
+ fRegParMu(6,0) = dummyMatrix(2,0);
+ fRegParMu(6,1) = dummyMatrix(5,0);
+ fRegParMu(6,2) = dummyMatrix(8,0);
+ fRegParMu(7,0) = dummyMatrix(2,1);
+ fRegParMu(7,1) = dummyMatrix(5,1);
+ fRegParMu(7,2) = dummyMatrix(8,1);
+ fRegParMu(8,0) = dummyMatrix(6,2);
+ fRegParMu(8,1) = dummyMatrix(5,2);
+ fRegParMu(8,2) = dummyMatrix(8,2);
break;
}
inFile->Close();
//-----------------------------------------------------------------------------
// This function smears track parameters according to streched cov. matrix
//-----------------------------------------------------------------------------
+ Double_t xref=xxsm[0]; xxsm[0]=0;
+
AliGausCorr *corgen = new AliGausCorr(cov,5);
TArrayD corr(5);
corgen->GetGaussN(corr);
+ // check on fP2(ESD) = fX*fP4-fP2(TPC)
+ // with fX=xref (not smeared), fP4=xx[4]+corr[4] e fP2=xx[2]+corr[2];
+ // if fP2(ESD)^2 > 1 -> resmear...
+ Double_t fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
+ while ( (fp2esd*fp2esd) > 1.0 ) {
+ Warning("SmearTrack()","Badly smeared track, retrying...");
+ corgen->GetGaussN(corr);
+ fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]);
+ }
+
delete corgen;
corgen = 0;
- for(Int_t l=0;l<5;l++) {
- xxsm[l] = xx[l]+corr[l];
- }
+ for(Int_t l=0;l<5;l++) xxsm[l] = xx[l]+corr[l];
return;
}