Tracking V1 ported to the HEAD
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Oct 2002 18:29:34 +0000 (18:29 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Oct 2002 18:29:34 +0000 (18:29 +0000)
ITS/AliITSComparisonV1.C [new file with mode: 0644]
ITS/AliITSPlotTracksV1.C
ITS/AliITSStoreFindableTracks.C [new file with mode: 0644]
ITS/AliITSStoreFindableTracksCompiled.C [new file with mode: 0644]
ITS/AliITSTrackerV1.cxx
ITS/AliITSTrackingV1.C
ITS/AliITStestV1.sh

diff --git a/ITS/AliITSComparisonV1.C b/ITS/AliITSComparisonV1.C
new file mode 100644 (file)
index 0000000..beac5ee
--- /dev/null
@@ -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();
+}
+
index ad583475ad2af02c1aaa8c94562000b5e240d5e9..7c7742b6d3007969b400c52437fa8817beb0ddfb 100644 (file)
 #include <iostream.h>
 #include <fstream.h>
 
 #include <iostream.h>
 #include <fstream.h>
 
-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 "<<ngood<<endl;  
-  Stat_t nfound=hfound->GetEntries(); cerr<<"Found tracks "<<nfound<<endl;
-  Stat_t nfake=hfake->GetEntries(); cerr<<"Fake tracks "<<nfake<<endl;       
-  gStyle->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 = "<<neglabs<<"\n";  // provvisoria
+       // Drawing
+       cerr << "Findable tracks   : " << hFindables->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 (file)
index 0000000..eb3bd42
--- /dev/null
@@ -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 (file)
index 0000000..19d6800
--- /dev/null
@@ -0,0 +1,182 @@
+#include <fstream.h>
+
+#include <TClassTable.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <TClonesArray.h>
+#include <TParticle.h>
+
+#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 <<endl;
+       if (nParticles <= 0) {
+               cerr << "Can't have <= 0 particles!" << endl;
+               return 3; 
+       }
+       AliITSRecPoint *recp = 0;
+       TClonesArray  *recPoints = ITS->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;
+}
index dc6a7c20645987d8d654c9bc4c2ad700d1befa9b..4b761734e4a743c88cbfbc7bfe34be8abe7d9da3 100644 (file)
@@ -15,6 +15,9 @@
  
 /*
 $Log$
  
 /*
 $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
 
 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");  
     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();
     if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
 
     cf->cd();
index c2b0f24bd45b76332e0647dd3b9146db379d78b4..c49c4784b3f9936d0dc8016f36c0384ced09558c 100644 (file)
@@ -2,9 +2,76 @@
 #include "iostream.h"
 #endif
 
 #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") {
 
 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) {
   ///////////////// 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;
   }
     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);
 
 // Connect the Root Galice file containing Geometry, Kine and Hits
    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
index 9eacef7ba7a41742ab7543fcf5b778078bcd2739..6087a0fc97a2de81bf391305e093b896ab9719f8 100755 (executable)
@@ -4,13 +4,19 @@
 ./AliITSDeleteOldFiles.sh
 
 # run the hit generation
 ./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
 
 # 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"
 
 # 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/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
 
 # 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 
 
 #
 # after all of the above you can run AliITSPlotTracksV1.C macro under