From 9040b78941059237d08148bf836a2660d42183e6 Mon Sep 17 00:00:00 2001 From: barbera Date: Tue, 22 Oct 2002 18:29:34 +0000 Subject: [PATCH] Tracking V1 ported to the HEAD --- ITS/AliITSComparisonV1.C | 158 ++++++++++++++ ITS/AliITSPlotTracksV1.C | 273 ++++++++++++------------ ITS/AliITSStoreFindableTracks.C | 8 + ITS/AliITSStoreFindableTracksCompiled.C | 182 ++++++++++++++++ ITS/AliITSTrackerV1.cxx | 5 +- ITS/AliITSTrackingV1.C | 79 ++++++- ITS/AliITStestV1.sh | 29 ++- 7 files changed, 583 insertions(+), 151 deletions(-) create mode 100644 ITS/AliITSComparisonV1.C create mode 100644 ITS/AliITSStoreFindableTracks.C create mode 100644 ITS/AliITSStoreFindableTracksCompiled.C diff --git a/ITS/AliITSComparisonV1.C b/ITS/AliITSComparisonV1.C new file mode 100644 index 00000000000..beac5ee9173 --- /dev/null +++ b/ITS/AliITSComparisonV1.C @@ -0,0 +1,158 @@ +#include "iostream.h" + +void AliITSComparisonV1 +(const char *foundfile = "itstracks.root", const char *truefile = "galice_tracks_0.root", Int_t evnum = 0) { + + // Make sure that ALICE objects are loaded + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + else { + delete gAlice; + gAlice=0; + } + + // Load the tree with track data + TFile *fileTrueTracks = new TFile(truefile); + TTree *treeTrueTracks = (TTree*)fileTrueTracks->Get("Tracks"); + TFile *fileFoundTracks = new TFile(foundfile); + TTree *treeFoundTracks = (TTree*)fileFoundTracks->Get(Form("TreeT%d", evnum)); + + Int_t nTrue = (Int_t)treeTrueTracks->GetEntries(); + Int_t *cnt = new Int_t[nTrue]; + for (Int_t i = 0; i < nTrue; i++) cnt[i] = 0; + + Int_t label, tpc_ok, entry, pdg_code; + Double_t px, py, pz, pt; + treeTrueTracks->SetBranchAddress("label", &label); + treeTrueTracks->SetBranchAddress("entry", &entry); + treeTrueTracks->SetBranchAddress("tpc_ok", &tpc_ok); + treeTrueTracks->SetBranchAddress("pdg_code", &pdg_code); + treeTrueTracks->SetBranchAddress("px", &px); + treeTrueTracks->SetBranchAddress("py", &py); + treeTrueTracks->SetBranchAddress("pz", &pz); + treeTrueTracks->SetBranchAddress("pt", &pt); + + AliITSIOTrack *iotrack = 0; + treeFoundTracks->SetBranchAddress("ITStracks", &iotrack); + + const Int_t npos = 36, nneg = 31; + Int_t pos[npos] = {2212, -11, -13, 211, 321, 3222, 213, 323, 10323, 3224, + 2224, 2214, -1114, -3112, -3312, 3224, -3114, -3314, 411, + 431, 413, 433, -15, 4232, 4222, 4322, 4422, 4412, 4432, + 4224, 4214, 4324, 4424, 4414, 4434, 4444}; + Int_t neg[nneg] = {2212, 11, 13, -211, -321, 3112, -213, -323, -10323, 3114, + 1114, -2224, -2214, 33112, -3222, 3114, 3314, 3334, -3224, + -411, -431, -413, -433, 15, -4422, -4432, -4214, -4324, + -4424, -4434, -444}; + + // Evaluation tree definition + Int_t labITS, labTPC, signC; + Double_t difpt, diflambda, difphi, Dz, Dr, Dtot, ptg; + TTree *treeEvaluation = new TTree("Eval", "Parameters for tracking evaluation"); + treeEvaluation->Branch("labITS" , &labITS , "labITS/I" ); + treeEvaluation->Branch("labTPC" , &labTPC , "labTPC/I" ); + treeEvaluation->Branch("difpt" , &difpt , "difpt/D" ); + treeEvaluation->Branch("diflambda", &diflambda, "diflambda/D"); + treeEvaluation->Branch("difphi" , &difphi , "difphi/D" ); + treeEvaluation->Branch("Dz" , &Dz , "Dz/D" ); + treeEvaluation->Branch("Dr" , &Dr , "Dr/D" ); + treeEvaluation->Branch("Dtot" , &Dtot , "Dtot/D" ); + treeEvaluation->Branch("ptg" , &ptg , "ptg/D" ); + treeEvaluation->Branch("signC" , &signC , "signC/I" ); + + // Make comparison + Double_t *var = 0; + Bool_t isGood; + Int_t nOK, trueEntry; + Double_t found_px, found_py, found_pz, found_pt; + Double_t found_tgl, found_lambda, found_phi, found_curv; + Double_t true_lambda, true_phi, true_px, true_py, true_pz, true_pt; + Double_t duepi = 2.*TMath::Pi(); + Bool_t ispos, isneg; + Int_t countpos = 0, countneg = 0, found_charge; + for (Int_t i = 0; i < treeFoundTracks->GetEntries(); i++) { + treeFoundTracks->GetEntry(i); + labITS = iotrack->GetLabel(); + labTPC = iotrack->GetTPCLabel(); + isGood = (labITS >= 0); + nOK = treeTrueTracks->Draw("px:py:pz", Form("label==%d && tpc_ok==1", abs(labITS)), "goff"); + if (!nOK) { + cerr << "ITS label not found among findable tracks:"; + cerr << " labITS = " << labITS; + cerr << " labTPC = " << labTPC; + cerr << endl; + isGood = kFALSE; + } + if (nOK > 1) { + cerr << "More than 1 tracks with label " << labITS << ": taking the first" << endl; + } + true_px = *treeTrueTracks->GetV1(); + true_py = *treeTrueTracks->GetV2(); + true_pz = *treeTrueTracks->GetV3(); + true_pt = TMath::Sqrt(true_px * true_px + true_py * true_py); + true_phi = TMath::ATan2(true_py, true_px); + if(true_phi < 0) true_phi += duepi; + true_phi *= 1000.0; + true_lambda = TMath::ATan(true_pz / true_pt) * 1000.0; + ptg = true_pt; + + // checks if two found good tracks have the same label in ITS + treeTrueTracks->Draw("entry", Form("label==%d && tpc_ok==1", abs(labITS)), "goff"); + trueEntry = (Int_t)*treeTrueTracks->GetV1(); + if (isGood && cnt[trueEntry] == 0) + cnt[trueEntry] = 1; + else if (isGood) { + cout << "Entry " << trueEntry << " already analyzed!" << endl; + continue; + } + + // charge matching + found_charge = (Int_t)iotrack->GetCharge(); + ispos = isneg = kFALSE; + for (Int_t j = 0; j < npos; j++) ispos = (pdg_code == pos[j]); + for (Int_t j = 0; j < nneg; j++) isneg = (pdg_code == neg[j]); + if (ispos) countpos++; + if (isneg) countneg++; + + // pt resolution (%) + found_px = iotrack->GetPx(); + found_py = iotrack->GetPy(); + found_pz = iotrack->GetPz(); + found_pt = TMath::Sqrt(found_px*found_px + found_py*found_py); + difpt = ((2. * found_pt - true_pt) / true_pt) * 100.; + + // lambda (mrad) + found_tgl = iotrack->GetStateTgl(); + found_lambda = TMath::ATan(found_tgl) * 1000.0; + diflambda = found_lambda - true_lambda; +// cout << "lambda " << found_lambda << " " << true_lambda << " " << diflambda << endl; + + // phi (mrad) + found_phi = TMath::ACos(found_px / found_pt); + if(found_phi > duepi) found_phi -= duepi; + if(found_phi < 0.) found_phi += duepi; + found_phi *= 1000.0; + difphi = found_phi - true_phi; +// cout << "phi " << found_phi << " " << true_phi << " " << difphi << endl; + + // impact parameters (microns) + Dr = iotrack->GetStateD() * 1.e4; + Dz = iotrack->GetDz() * 1.e4; + Dtot = sqrt(Dr*Dr + Dz*Dz); + + // curvature sign + found_curv = iotrack->GetStateC(); + signC = (found_curv > 0.) ? 1 : -1; + + // fill tree + treeEvaluation->Fill(); + } + + + TFile *outFile = new TFile("AliITSComparisonV1.root", "recreate"); + treeEvaluation->Write("Eval", TObject::kOverwrite); + outFile->Close(); +} + diff --git a/ITS/AliITSPlotTracksV1.C b/ITS/AliITSPlotTracksV1.C index ad583475ad2..7c7742b6d30 100644 --- a/ITS/AliITSPlotTracksV1.C +++ b/ITS/AliITSPlotTracksV1.C @@ -1,143 +1,140 @@ #include #include -void AliITSPlotTracksV1(){ +void AliITSPlotTracksV1() +{ + TFile *file = new TFile("AliITSComparisonV1.root"); + TTree *tree = (TTree*)file->Get("Eval"); + + if (!tree) { + cerr << "Evaluation tree not found!" << endl; + return; + } + + // histogram definition - I (efficiency) + TH1D *hFindables = new TH1D("hFindables", "Findable tracks", 10, 0, 2); + TH1D *hGood = new TH1D("hGood", "Good found tracks", 10, 0, 2); + TH1D *hFake = new TH1D("hFake", "Fake found tracks", 10, 0, 2); + TH1D *hRatioG = new TH1D("hRatioG", "", 10, 0, 2); + TH1D *hRatioF = new TH1D("hRatioF", "", 10, 0, 2); + hGood->Sumw2(); + hFake->Sumw2(); + hFindables->Sumw2(); + hRatioG->SetLineColor(kBlue); + hRatioG->SetLineWidth(2); + hRatioF->SetLineColor(kRed); + hRatioF->SetLineWidth(2); + + // histograms definition - II (resolution) + TH1D *hPhi = new TH1D("hPhi", "#Phi resolution", 50, -15., 15.); + TH1D *hLambda = new TH1D("hLambda", "#lambda resolution", 50, -15., 15.); + TH1D *hPt = new TH1D("hPt", "Relative Pt resolution", 40, -10., 10.); + TH1D *hDtot = new TH1D("hDtot", "Total impact parameter distribution", 100, 0., 2000.); + TH1D *hDr = new TH1D("hDr", "Transv. impact parameter distribution", 50, -1000., 1000.); + TH1D *hDz = new TH1D("hDz", "Long. impact parameter distribution", 50, -1000., 1000.); + hPhi->SetFillColor(4); + hLambda->SetFillColor(4); + hPt->SetFillColor(2); + hDtot->SetFillColor(6); + hDr->SetFillColor(kGreen); + hDz->SetFillColor(kBlue); + + // Evaluation tree settings + Int_t labITS, labTPC, signC; + Int_t i, j, tot = (Int_t)tree->GetEntries(); + Double_t difpt, diflambda, difphi, Dz, Dr, Dtot, ptg; + tree->SetBranchAddress("labITS" , &labITS ); + tree->SetBranchAddress("labTPC" , &labTPC ); + tree->SetBranchAddress("difpt" , &difpt ); + tree->SetBranchAddress("diflambda", &diflambda); + tree->SetBranchAddress("difphi" , &difphi ); + tree->SetBranchAddress("Dz" , &Dz ); + tree->SetBranchAddress("Dr" , &Dr ); + tree->SetBranchAddress("Dtot" , &Dtot ); + tree->SetBranchAddress("ptg" , &ptg ); + tree->SetBranchAddress("signC" , &signC ); + + // Filling the histogram of findable tracks (w.r. to momentum) + for(i = 0; i < tot; i++) { + tree->GetEntry(i); + hFindables->Fill(ptg); + } + + // Filling the evaluation histograms + Int_t neglabs = 0; + for(i = 0; i < tot; i++) { + tree->GetEntry(i); +// if(signC<0) continue; + if (labITS < 0) neglabs++; + if (labITS >= 0) { + hGood->Fill(ptg); + hPt->Fill(difpt); + hLambda->Fill(diflambda); + hPhi->Fill(difphi); + hDtot->Fill(Dtot); + hDr->Fill(Dr); + hDz->Fill(Dz); + } + else { + hFake->Fill(ptg); + neglabs++; + } + } - ifstream in ("AliITSTra.out"); - - TVector DataOut(10); - - ///////////////////////////////// Histograms definition /////////////////////////////////////////// - - - TH1F *hp=new TH1F("hp","PHI resolution",50,-15.,15.); hp->SetFillColor(4); - TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-15.,15.); hl->SetFillColor(4); - - // TH1F *hp=new TH1F("hp","PHI resolution",50,-5.,5.); hp->SetFillColor(4); - //TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-5.,5.); hl->SetFillColor(4); - - // TH1F *hpt=new TH1F("hpt","Relative Pt resolution",50,-5.,5.); - TH1F *hpt=new TH1F("hpt","Relative Pt resolution",40,-10.,10.); - hpt->SetFillColor(2); - TH1F *hd=new TH1F("hd","Impact parameter distribution ",100,0,2000); - hd->SetFillColor(6); - - TH1F *hdr=new TH1F("hdr","Dr ",50,-1000,1000); - hdr->SetFillColor(kGreen); - TH1F *hdz=new TH1F("hdz","Dz ",50,-1000,1000); - hdz->SetFillColor(kBlue); - - /* - TH1F *hdr=new TH1F("hdr","Dr ",50,-100,100); - hdr->SetFillColor(kGreen); - TH1F *hdz=new TH1F("hdz","Dz ",50,-1000,1000); - hdz->SetFillColor(kBlue); - */ - - TH1F *hgood=new TH1F("hgood","Good tracks",10,0,2); - hgood->Sumw2(); // aggiunto il 15-01-01 - TH1F *hfound=new TH1F("hfound","Found tracks",10,0,2); - hfound->Sumw2(); // aggiunto il 15-01-01 - TH1F *hfake=new TH1F("hfake","Fake tracks",10,0,2); - hfake->Sumw2(); - TH1F *hg=new TH1F("hg","",10,0,2); //efficiency for good tracks - hg->SetLineColor(4); hg->SetLineWidth(2); - TH1F *hf=new TH1F("hf","Efficiency for fake tracks",10,0,2); - /*hf->SetFillColor(1); hf->SetFillStyle(3013);*/ hf->SetLineColor(4); hf->SetLineWidth(2); - - /////////////////////////////////////////////////////////////////////////////////////////////////// - - ifstream in1 ("AliITSTrag.out"); - - Double_t ptg; - for(;;) { - in1 >> ptg; - if( in1.eof() ) break; - hgood->Fill(ptg); - } - in1.close(); - - Int_t neglabs=0; - for (;;){ - for (int r=0; r<10; r++) in>>DataOut(r); - if( in.eof() ) break; - - Double_t ptg=DataOut(0); Double_t labITS=DataOut(1); Double_t labTPC=DataOut(2); Double_t ptperc=DataOut(3); - Double_t deltalam=DataOut(4); Double_t deltaphi=DataOut(5); - Double_t Dtot=DataOut(6); Double_t Dr=DataOut(7); Double_t Dz=DataOut(8); - Double_t signC=DataOut(9); - - // if(signC<0) continue; - - if(labITS<0) neglabs++; - if(labITS>=0) hfound->Fill(ptg); else { hfake->Fill(ptg);} - - if(labITS>=0 ) { - hpt->Fill(ptperc); - hl->Fill(deltalam); - hp->Fill(deltaphi); - hd->Fill(Dtot); - hdr->Fill(Dr); - hdz->Fill(Dz); - } - } - - - in.close(); - Stat_t ngood=hgood->GetEntries(); cerr<<"Good tracks "<GetEntries(); cerr<<"Found tracks "<GetEntries(); cerr<<"Fake tracks "<SetOptStat(111110); - gStyle->SetOptFit(1); - TCanvas *c1=new TCanvas("c1","",0,0,700,700); - TPad *p1=new TPad("p1","",0,0.5,0.5,1); p1->Draw(); hp->SetXTitle("(mrad)"); - p1->cd(); hp->Draw(); hp->Fit("gaus"); c1->cd(); - TPad *p2=new TPad("p2","",0.5,0.5,1,1); p2->Draw(); hl->SetXTitle("(mrad)"); - p2->cd(); hl->Draw(); hl->Fit("gaus"); c1->cd(); - TPad *p3=new TPad("p3","",0,0,0.5,0.5); p3->Draw(); hpt->SetXTitle("(%)"); - p3->cd(); hpt->Draw(); hpt->Fit("gaus"); c1->cd(); - TPad *p4=new TPad("p4","",0.5,0,1,0.5); p4->Draw(); hd->SetXTitle("(micron)"); - p4->cd(); hd->Draw(); c1->cd(); - - TCanvas *c3=new TCanvas("c3","",200,200,800,500); - hfound->Print("all"); // aggiunto il 16-01-01 - hgood->Print("all"); // aggiunto il 16-01-01 - TPad *p7=new TPad("p7","",0,0,0.333,1); p7->Draw(); p7->cd(); hfound->Draw(); c3->cd(); - TPad *p8=new TPad("p8","",0.333,0,0.666,1); p8->Draw(); p8->cd(); hfake->Draw(); c3->cd(); - TPad *p9=new TPad("p9","",0.666,0,1,1); p9->Draw(); p9->cd(); hgood->Draw(); c3->cd(); - - TCanvas *c4=new TCanvas("c4","",300,300,800,500); - hg->Divide(hfound,hgood,1.,1.); //,"b"); - hf->Divide(hfake,hgood,1.,1.); //,"b"); - hg->SetMaximum(1.4); - hg->SetYTitle("Tracking efficiency"); - hg->SetXTitle("Pt (GeV/c)"); - hg->Print("all"); // aggiunto il 16-01-01 - hg->Draw(); // to not plot the erro bars hg->Draw("histo"); - TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4); - line1->Draw("same"); - TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4); - line2->Draw("histosame"); - - - hf->SetFillColor(1); - hf->SetFillStyle(3013); - hf->SetLineColor(2); - hf->SetLineWidth(2); - hf->Draw("same"); // to not plot the error bars hf->Draw("histosame"); - - TText *text = new TText(0.461176,0.248448,"Fake tracks"); - text->SetTextSize(0.05); - text->Draw(); - text = new TText(0.453919,1.11408,"Good tracks"); - text->SetTextSize(0.05); - text->Draw(); - - TCanvas *c2=new TCanvas("c2","",100,100,700,400); - TPad *p5=new TPad("p5","",0,0,0.5,1); p5->Draw(); hdr->SetXTitle("(micron)"); - p5->cd(); hdr->Draw(); hdr->Fit("gaus"); c2->cd(); - TPad *p6=new TPad("p6","",0.5,0,1,1); p6->Draw(); hdz->SetXTitle("(micron)"); - p6->cd(); hdz->Draw(); hdz->Fit("gaus"); c2->cd(); - - cout<<"neglabs = "<GetEntries() << endl; + cerr << "Good found tracks : " << hGood->GetEntries() << endl; + cerr << "Fake found tracks : " << hFake->GetEntries() << endl; + + gStyle->SetOptStat(111110); + gStyle->SetOptFit(1); + + TCanvas *c1 = new TCanvas("c1","Parameter resolutions",0,0,700,700); + c1->Divide(2, 2, 0.001, 0.001); + c1->cd(1); hPhi->SetXTitle("(mrad)"); hPhi->Draw(); hPhi->Fit("gaus", "Q"); + c1->cd(2); hLambda->SetXTitle("(mrad)"); hLambda->Draw(); hLambda->Fit("gaus", "Q"); + c1->cd(3); hPt->SetXTitle("(%)"); hPt->Draw(); hPt->Fit("gaus", "Q"); + c1->cd(4); hDtot->SetXTitle("(micron)"); hDtot->Draw(); + c1->Update(); + + TCanvas *c2 = new TCanvas("c2","Impact parameters resolution",100,100,700,400); + c2->Divide(2, 1, 0.001, 0.001); + c2->cd(1); hDr->SetXTitle("(micron)"); hDr->Draw(); hDr->Fit("gaus", "Q"); + c2->cd(2); hDz->SetXTitle("(micron)"); hDz->Draw(); hDz->Fit("gaus", "Q"); + c2->Update(); + + TCanvas *c3 = new TCanvas("c3","Momentum distributions",200,200,800,500); + c3->Divide(3, 1, 0.001, 0.001); + c3->cd(1); hGood->Draw(); + c3->cd(2); hFake->Draw(); + c3->cd(3); hFindables->Draw(); + c3->Update(); + + TCanvas *c4 = new TCanvas("c4","Tracking efficiency",300,300,800,500); + hRatioG->Divide(hGood, hFindables, 1., 1.); + hRatioF->Divide(hFake, hFindables, 1., 1.); + hRatioG->SetMaximum(1.4); + hRatioG->SetYTitle("Tracking efficiency"); + hRatioG->SetXTitle("Pt (GeV/c)"); + hRatioG->Draw(); // to not plot the erro bars hg->Draw("histo"); + hRatioF->SetFillColor(1); + hRatioF->SetFillStyle(3013); + hRatioF->SetLineColor(2); + hRatioF->SetLineWidth(2); + hRatioF->Draw("same"); // to not plot the error bars hRatioF->Draw("histosame"); + // a line to mark the best efficiency + TLine *line1 = new TLine(0,1.0,2,1.0); line1->SetLineStyle(4); + line1->Draw("same"); + TLine *line2 = new TLine(0,0.9,2,0.9); line2->SetLineStyle(4); + line2->Draw("histosame"); + // a text explanation + TText *text = new TText(0.461176,0.248448,"Fake tracks"); + text->SetTextSize(0.05); + text->Draw(); + text = new TText(0.453919,1.11408,"Good tracks"); + text->SetTextSize(0.05); + text->Draw(); + c4->Update(); + + cout << "neglabs = " << neglabs << endl; // provvisoria } diff --git a/ITS/AliITSStoreFindableTracks.C b/ITS/AliITSStoreFindableTracks.C new file mode 100644 index 00000000000..eb3bd425471 --- /dev/null +++ b/ITS/AliITSStoreFindableTracks.C @@ -0,0 +1,8 @@ +Int_t AliITSStoreFindableTracks +(Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0) +{ + gSystem->SetIncludePath("-I- -I$ALICE_ROOT/ITS -I$ALICE_ROOT/STEER"); + + gROOT->LoadMacro("$ALICE_ROOT/ITS/AliITSStoreFindableTracksCompiled.C+"); + AliITSStoreFindableTracksCompiled(nMinClusters, evname, evnum); +} diff --git a/ITS/AliITSStoreFindableTracksCompiled.C b/ITS/AliITSStoreFindableTracksCompiled.C new file mode 100644 index 00000000000..19d6800934f --- /dev/null +++ b/ITS/AliITSStoreFindableTracksCompiled.C @@ -0,0 +1,182 @@ +#include + +#include +#include +#include +#include +#include +#include + +#include "AliRun.h" +#include "AliITS.h" +#include "AliITSgeom.h" +#include "AliITSRecPoint.h" + +Int_t AliITSStoreFindableTracksCompiled +(Int_t nMinClusters = 5, const Text_t *evname = "galice", Int_t evnum = 0) +{ + // Make sure that ALICE objects are loaded + if (gAlice) { + delete gAlice; + gAlice = 0; + } + + // Define the names of all involved files + TString strEventFile(evname); + TString strOutputFile(evname); + strEventFile.Append(".root"); + strOutputFile.Append("_tracks_"); + strOutputFile += evnum; + strOutputFile.Append(".root"); + + // Connect the Root Galice file containing Geometry, Kine and Hits + TFile *fileEvent = (TFile*)gROOT->GetListOfFiles()->FindObject(strEventFile); + if (!fileEvent) fileEvent = new TFile(strEventFile,"UPDATE"); + + // Get AliRun object from file + gAlice = (AliRun*)fileEvent->Get("gAlice"); + if (gAlice) cout << "OK, found an AliRun object in file" << endl; + + // Get ITS related objects and data + AliITS* ITS =(AliITS *)gAlice->GetDetector("ITS"); + if (!ITS) { + cerr << "ITS object not found!" << endl; + return 1; + } + AliITSgeom *geometry = ITS->GetITSgeom(); + if (!geometry) { + cerr << "ITS geometry object not found!" << endl; + return 2; + } + + // Count the number of modules per layer + Int_t nLadders, nDetectors, mod_min[6], mod_max[6]; + for(Int_t i = 0; i < 6; i++) { + nLadders = geometry->GetNladders(i + 1); + nDetectors = geometry->GetNdetectors(i + 1); + mod_min[i] = geometry->GetModuleIndex(i + 1, 1, 1); + mod_max[i] = geometry->GetModuleIndex(i + 1, nLadders, nDetectors); + } + + // Load event and ITS recpoints + Int_t nParticles = gAlice->GetEvent(evnum); + cout << "Event number: " << evnum << endl; + cout << "# particles : " << nParticles <RecPoints(); + TObjArray *particles = gAlice->Particles(); + Int_t nTracks = gAlice->GetNtrack(); //FCA correction + Bool_t *hitITSLayer[6]; + for (Int_t i = 0; i < 6; i++) { + hitITSLayer[i] = new Bool_t[nTracks]; + for (Int_t j = 0; j < nTracks; j++) hitITSLayer[i][j] = kFALSE; + } + + // Load recpoints in event + TTree *TR = gAlice->TreeR(); + if (!TR) { + cerr << "TreeR object not found!" << endl; + return 4; + } + + // Scan recpoints and define findable tracks + Int_t nModules = (Int_t)TR->GetEntries(), nPoints = 0; + cout << "Found " << nModules; + cout << " entries in the TreeR (must be one per module!)" << endl; + for (Int_t layer = 1; layer <= 6; layer++) { + for (Int_t mod = mod_min[layer - 1]; mod <= mod_max[layer - 1]; mod++) { + ITS->ResetRecPoints(); + TR->GetEntry(mod); + nPoints = recPoints->GetEntries(); + if(!nPoints) { + cout << "Module " << mod << " is empty..." << endl; + continue; + } + for (Int_t point = 0; point < nPoints; point++) { + recp = (AliITSRecPoint*)recPoints->UncheckedAt(point); + for (Int_t it = 0; it < 3; it++) { + Int_t track = recp->GetLabel(it); + if(track < 0) continue; + if(track > nTracks) { + cout << "Found track index " << track; + cout << " whilw gAlice->GetNtrack() = " << nTracks << endl; + continue; + } + hitITSLayer[layer - 1][track] = kTRUE; + } // loop over recpoint labels + } //loop over points + } //loop over modules + } //loop over layers + + // Scan the file of tracks in TPC to retrieve the findable TPC tracks + TString strLabelsTPC; + Int_t label, pdg_code, nFindablesTPC = 0; + Double_t dummy; + ifstream tpc("good_tracks_tpc"); + while (tpc >> label >> pdg_code) { + for (Int_t i = 0; i < 6; i++) tpc >> dummy; + nFindablesTPC++; + strLabelsTPC.Append(Form("[%d]", label)); + } + + // Define the TTree with tracks data by means of a set of variables + Int_t nFindablesITS = 0, nFindablesITSTPC = 0; + Int_t nhits, tpc_ok, entry = 0; + Double_t vx, vy, vz; + Double_t px, py, pz, pt; + + TTree *tree = new TTree("Tracks", "Findable tracks in ITS"); + + tree->Branch("vx", &vx, "vx/D"); + tree->Branch("vy", &vy, "vy/D"); + tree->Branch("vz", &vz, "vz/D"); + tree->Branch("px", &px, "px/D"); + tree->Branch("py", &py, "py/D"); + tree->Branch("pz", &pz, "pz/D"); + tree->Branch("pt", &pt, "pt/D"); + tree->Branch("label", &label, "label/I"); + tree->Branch("entry", &entry, "entry/I"); + tree->Branch("pdg_code", &pdg_code, "pdg_code/I"); + tree->Branch("nhits", &nhits, "nhits/I"); + tree->Branch("tpc_ok", &tpc_ok, "tpc_ok/I"); + + // Fill the tree + cout << endl; + TParticle *p = 0; + for (Int_t i = 0; i < nTracks; i++) { + p = gAlice->Particle(i); + px = p->Px(); + py = p->Py(); + pz = p->Pz(); + pt = p->Pt(); + vx = p->Vx(); + vy = p->Vy(); + vz = p->Vz(); + nhits = 0; + for (Int_t j = 0; j < 6; j++) if (hitITSLayer[j][i]) nhits++; + if (nhits < nMinClusters) continue; + cout << "Track " << i << " stored\r" << flush; + tpc_ok = (strLabelsTPC.Contains(Form("[%d]", i))); + pdg_code = p->GetPdgCode(); + label = i; + nFindablesITS++; + if (tpc_ok) nFindablesITSTPC++; + tree->Fill(); + entry++; + } + + // Save into a file + TFile *fileOutput = new TFile(strOutputFile, "recreate"); + tree->Write(); + fileOutput->Close(); + + cout << "# findable tracks in TPC : " << nFindablesTPC << endl; + cout << "# findable tracks in ITS : " << nFindablesITS << endl; + cout << "# findable tracks in ITS+TPC : " << nFindablesITSTPC << endl; + + return 0; +} diff --git a/ITS/AliITSTrackerV1.cxx b/ITS/AliITSTrackerV1.cxx index dc6a7c20645..4b761734e4a 100644 --- a/ITS/AliITSTrackerV1.cxx +++ b/ITS/AliITSTrackerV1.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.22 2002/10/22 14:45:36 alibrary +Introducing Riostream.h + Revision 1.21 2002/02/05 09:12:26 hristov Small mods for gcc 3.02 @@ -550,7 +553,7 @@ void AliITSTrackerV1::DoTracking(Int_t evNumber,Int_t minTr,Int_t maxTr, kkprov->SetConvConst(100/0.299792458/0.2/fFieldFactor); TFile *cf=TFile::Open("AliTPCclusters.root"); - AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60"); + AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60_150x60"); if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();} cf->cd(); diff --git a/ITS/AliITSTrackingV1.C b/ITS/AliITSTrackingV1.C index c2b0f24bd45..c49c4784b3f 100644 --- a/ITS/AliITSTrackingV1.C +++ b/ITS/AliITSTrackingV1.C @@ -2,9 +2,76 @@ #include "iostream.h" #endif +Bool_t TPCSortTracks(Int_t event = 0) +{ + TFile *fileTracks = TFile::Open("AliTPCtracks.root"); + TFile *fileClusters = TFile::Open("AliTPCclusters.root"); + TFile *fileEvent = TFile::Open("galice.root"); + + // get TPC parameterization + AliTPCParam *param=(AliTPCParam *)fileEvent->Get("75x40_100x60_150x60"); + if (!param) { + cerr << "(TPCSortTracks) ERROR: TPC parameters have not been found!" << endl; + return kFALSE; + } + + // read and sort tracks + Int_t i; + TSortedList tracks_list; + AliTPCtrack *iotrack = 0; + TTree *tracktree = (TTree*)fileTracks->Get(Form("TreeT_TPC_%d", event)); + Int_t nentr = (Int_t)tracktree->GetEntries(); + for (i = 0; i < nentr; i++) { + iotrack = new AliTPCtrack; + tracktree->SetBranchAddress("tracks", &iotrack); + tracktree->GetEvent(i); + tracks_list.Add(iotrack); + } + delete tracktree; + + // assign to each track its GEANT label + fileClusters->cd(); + AliTPCtracker *tracker = new AliTPCtracker(param, event); + tracker->LoadInnerSectors(); + tracker->LoadOuterSectors(); + TListIter iter(&tracks_list); + for (i = 0; i < nentr; i++) { + iotrack = (AliTPCtrack*)iter.Next(); + if (!iotrack) { + cerr << "(TPCSortTracks) WARNING: Track no. " << i << " is NULL!!!" << endl; + continue; + } + tracker->CookLabel(iotrack, 0.1); + } + delete tracker; + + // create the new TTree of TPC tracks sorted w.r. to Pt + tracktree = new TTree(Form("TreeT_TPC_%d", event),"Tree with TPC tracks sorted w.r to pt"); + tracktree->Branch("tracks", "AliTPCtrack", &iotrack, 32000, 0); + iter.Reset(); + for (i = 0; i < nentr; i++) { + iotrack = (AliTPCtrack*)iter.Next(); + if (!iotrack) { + cerr << "(TPCSortTracks) WARNING: Track no. " << i << " is NULL!!!" << endl; + continue; + } + tracktree->Fill(); + } + + // save the new tree into new file + TFile *fileOutput = TFile::Open("AliTPCtracksSorted.root","recreate"); + tracktree->Write(); + fileOutput->Close(); + fileEvent->Close(); + fileClusters->Close(); + fileTracks->Close(); + + return kTRUE; +} + + void AliITSTrackingV1(Int_t evNumber1=0,Int_t evNumber2=0, Int_t min_t=-1, Int_t max_t=0,Bool_t flagvert=1, Bool_t realmass=0, const char *filename="galice.root") { - ///////////////// Dynamically link some shared libs //////////////////////////////// if (gClassTable->GetID("AliRun") < 0) { @@ -14,6 +81,16 @@ void AliITSTrackingV1(Int_t evNumber1=0,Int_t evNumber2=0, Int_t min_t=-1, Int_t delete gAlice; gAlice=0; } + + cout << "Sorting TPC tracks w.r. to transverse momentum..."; + Bool_t success_sorting = TPCSortTracks(); + if (success_sorting) { + cout << "DONE!" << endl; + } + else { + cout << "Some error occurred..." << endl; + return 1; + } // Connect the Root Galice file containing Geometry, Kine and Hits TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename); diff --git a/ITS/AliITStestV1.sh b/ITS/AliITStestV1.sh index 9eacef7ba7a..6087a0fc97a 100755 --- a/ITS/AliITStestV1.sh +++ b/ITS/AliITStestV1.sh @@ -4,13 +4,19 @@ ./AliITSDeleteOldFiles.sh # run the hit generation -aliroot -q -b "$ALICE_ROOT/macros/grun.C($1)" +aliroot -q -b "$ALICE_ROOT/macros/grun.C" # digitize TPC -aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($1)" +aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C" -# prepare TPC tracks for matching with the ITS -aliroot -q -b "$ALICE_ROOT/ITS/AliTPCTracking4ITS.C($1)" +# find clusters in TPC +aliroot -q -b "$ALICE_ROOT/TPC/AliTPCFindClusters.C" + +# find tracks in TPC +aliroot -q -b "$ALICE_ROOT/TPC/AliTPCFindTracks.C" + +# do TPC tracking comparison +aliroot -q -b "$ALICE_ROOT/TPC/AliTPCComparison.C" # create summable digits for the ITS aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2SDigits.C" @@ -19,16 +25,17 @@ aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2SDigits.C" aliroot -q -b "$ALICE_ROOT/ITS/AliITSSDigits2Digits.C" # create reconstructed points for the ITS -aliroot -q -b "$ALICE_ROOT/ITS/AliITSDigits2RecPoints.C(0,$[$1-1])" - -# prepare for tracking -aliroot -q -b "$ALICE_ROOT/ITS/AliITSTracksV1.C(0,$[$1-1])" +aliroot -q -b "$ALICE_ROOT/ITS/AliITSDigits2RecPoints.C" # do the tracking -aliroot -q -b "$ALICE_ROOT/ITS/AliITSTrackingV1.C(0,$[$1-1])" +aliroot -q -b "$ALICE_ROOT/ITS/AliITSTrackingV1.C" + +# prepare for comparison +# aliroot -q -b "$ALICE_ROOT/ITS/AliITSTracksV1.C" +aliroot -q -b "$ALICE_ROOT/ITS/AliITSStoreFindableTracks.C" -# do the comparison -aliroot -q -b "$ALICE_ROOT/ITS/AliITSComparisonV1.C(0,$[$1-1])" +# do ITS tracking comparison +aliroot -q -b "$ALICE_ROOT/ITS/AliITSComparisonV1.C" # # after all of the above you can run AliITSPlotTracksV1.C macro under -- 2.43.0