Tracking V1 ported to the HEAD
[u/mrichter/AliRoot.git] / ITS / AliITSComparisonV1.C
1 #include "iostream.h"
2
3 void AliITSComparisonV1
4 (const char *foundfile = "itstracks.root", const char *truefile = "galice_tracks_0.root", Int_t evnum = 0) {
5
6         // Make sure that ALICE objects are loaded
7         if (gClassTable->GetID("AliRun") < 0) {
8                 gROOT->LoadMacro("loadlibs.C");
9                 loadlibs();
10         }  
11         else {
12                 delete gAlice;
13                 gAlice=0;
14         }
15         
16         // Load the tree with track data
17         TFile *fileTrueTracks = new TFile(truefile);
18         TTree *treeTrueTracks = (TTree*)fileTrueTracks->Get("Tracks");
19         TFile *fileFoundTracks = new TFile(foundfile);
20         TTree *treeFoundTracks = (TTree*)fileFoundTracks->Get(Form("TreeT%d", evnum));
21         
22         Int_t nTrue = (Int_t)treeTrueTracks->GetEntries();
23         Int_t *cnt = new Int_t[nTrue];
24         for (Int_t i = 0; i < nTrue; i++) cnt[i] = 0;
25                 
26         Int_t label, tpc_ok, entry, pdg_code;
27         Double_t px, py, pz, pt;
28         treeTrueTracks->SetBranchAddress("label", &label);
29         treeTrueTracks->SetBranchAddress("entry", &entry);
30         treeTrueTracks->SetBranchAddress("tpc_ok", &tpc_ok);
31         treeTrueTracks->SetBranchAddress("pdg_code", &pdg_code);
32         treeTrueTracks->SetBranchAddress("px", &px);
33         treeTrueTracks->SetBranchAddress("py", &py);
34         treeTrueTracks->SetBranchAddress("pz", &pz);
35         treeTrueTracks->SetBranchAddress("pt", &pt);
36         
37         AliITSIOTrack *iotrack = 0;
38         treeFoundTracks->SetBranchAddress("ITStracks", &iotrack);
39         
40         const Int_t npos = 36, nneg = 31;
41         Int_t pos[npos] = {2212, -11, -13, 211, 321, 3222, 213, 323, 10323, 3224, 
42                            2224, 2214, -1114, -3112, -3312, 3224, -3114, -3314, 411, 
43                            431, 413, 433, -15, 4232, 4222, 4322, 4422, 4412, 4432, 
44                            4224, 4214, 4324, 4424, 4414, 4434, 4444};
45         Int_t neg[nneg] = {2212, 11, 13, -211, -321, 3112, -213, -323, -10323, 3114, 
46                            1114, -2224, -2214, 33112, -3222, 3114, 3314, 3334, -3224,
47                            -411, -431, -413, -433, 15, -4422, -4432, -4214, -4324, 
48                            -4424, -4434, -444};
49                            
50         // Evaluation tree definition 
51         Int_t labITS, labTPC, signC;
52         Double_t difpt, diflambda, difphi, Dz, Dr, Dtot, ptg;
53         TTree *treeEvaluation = new TTree("Eval", "Parameters for tracking evaluation");
54         treeEvaluation->Branch("labITS"   , &labITS   , "labITS/I"   );
55         treeEvaluation->Branch("labTPC"   , &labTPC   , "labTPC/I"   );
56         treeEvaluation->Branch("difpt"    , &difpt    , "difpt/D"    ); 
57         treeEvaluation->Branch("diflambda", &diflambda, "diflambda/D");
58         treeEvaluation->Branch("difphi"   , &difphi   , "difphi/D"   );
59         treeEvaluation->Branch("Dz"       , &Dz       , "Dz/D"       );
60         treeEvaluation->Branch("Dr"       , &Dr       , "Dr/D"       );
61         treeEvaluation->Branch("Dtot"     , &Dtot     , "Dtot/D"     );
62         treeEvaluation->Branch("ptg"      , &ptg      , "ptg/D"      );
63         treeEvaluation->Branch("signC"    , &signC    , "signC/I"    );
64         
65         // Make comparison
66         Double_t *var = 0;
67         Bool_t isGood;
68         Int_t nOK, trueEntry;
69         Double_t found_px, found_py, found_pz, found_pt;
70         Double_t found_tgl, found_lambda, found_phi, found_curv;
71         Double_t true_lambda, true_phi, true_px, true_py, true_pz, true_pt;
72         Double_t duepi = 2.*TMath::Pi();         
73         Bool_t ispos, isneg;
74         Int_t countpos = 0, countneg = 0, found_charge;
75         for (Int_t i = 0; i < treeFoundTracks->GetEntries(); i++) {
76                 treeFoundTracks->GetEntry(i);
77                 labITS = iotrack->GetLabel();
78                 labTPC = iotrack->GetTPCLabel();
79                 isGood = (labITS >= 0);
80                 nOK = treeTrueTracks->Draw("px:py:pz", Form("label==%d && tpc_ok==1", abs(labITS)), "goff");
81                 if (!nOK) {
82                         cerr << "ITS label not found among findable tracks:";
83                         cerr << "   labITS = " << labITS;
84                         cerr << "   labTPC = " << labTPC;
85                         cerr << endl;
86                         isGood = kFALSE;
87                 }
88                 if (nOK > 1) {
89                         cerr << "More than 1 tracks with label " << labITS << ": taking the first" << endl;
90                 }
91                 true_px = *treeTrueTracks->GetV1();
92                 true_py = *treeTrueTracks->GetV2();
93                 true_pz = *treeTrueTracks->GetV3();
94                 true_pt = TMath::Sqrt(true_px * true_px + true_py * true_py);
95                 true_phi = TMath::ATan2(true_py, true_px);
96                 if(true_phi < 0) true_phi += duepi;
97                 true_phi *= 1000.0;
98                 true_lambda = TMath::ATan(true_pz / true_pt) * 1000.0;
99                 ptg = true_pt;
100                 
101                 // checks if two found good tracks have the same label in ITS
102                 treeTrueTracks->Draw("entry", Form("label==%d && tpc_ok==1", abs(labITS)), "goff");
103                 trueEntry = (Int_t)*treeTrueTracks->GetV1();
104                 if (isGood && cnt[trueEntry] == 0) 
105                         cnt[trueEntry] = 1;
106                 else if (isGood) {
107                         cout << "Entry " << trueEntry << " already analyzed!" << endl;
108                         continue;       
109                 }
110                                 
111                 // charge matching
112                 found_charge = (Int_t)iotrack->GetCharge();
113                 ispos = isneg = kFALSE;
114                 for (Int_t j = 0; j < npos; j++) ispos = (pdg_code == pos[j]);
115                 for (Int_t j = 0; j < nneg; j++) isneg = (pdg_code == neg[j]);
116                 if (ispos) countpos++;
117                 if (isneg) countneg++;
118                 
119                 // pt resolution (%)
120                 found_px = iotrack->GetPx();
121                 found_py = iotrack->GetPy();
122                 found_pz = iotrack->GetPz();
123                 found_pt = TMath::Sqrt(found_px*found_px + found_py*found_py);
124                 difpt = ((2. * found_pt - true_pt) / true_pt) * 100.;
125                         
126                 // lambda (mrad)
127                 found_tgl = iotrack->GetStateTgl();
128                 found_lambda = TMath::ATan(found_tgl) * 1000.0;
129                 diflambda = found_lambda - true_lambda;
130 //              cout << "lambda " << found_lambda << " " << true_lambda << " " << diflambda << endl;
131                 
132                 // phi (mrad)
133                 found_phi = TMath::ACos(found_px / found_pt);
134                 if(found_phi > duepi) found_phi -= duepi;
135                 if(found_phi < 0.) found_phi += duepi;
136                 found_phi *= 1000.0;      
137                 difphi = found_phi - true_phi;
138 //              cout << "phi " << found_phi << " " << true_phi << " " << difphi << endl;
139                 
140                 // impact parameters (microns)
141                 Dr = iotrack->GetStateD() * 1.e4;
142                 Dz = iotrack->GetDz() * 1.e4;
143                 Dtot = sqrt(Dr*Dr + Dz*Dz);
144                 
145                 // curvature sign
146                 found_curv = iotrack->GetStateC();
147                 signC = (found_curv > 0.) ? 1 : -1;
148
149                 // fill tree
150                 treeEvaluation->Fill();
151         }
152
153         
154         TFile *outFile = new TFile("AliITSComparisonV1.root", "recreate");
155         treeEvaluation->Write("Eval", TObject::kOverwrite);
156         outFile->Close();
157 }
158