ITS/oldmacros/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
  83 //                        cerr << "ITS label not found among findable tracks:";
  84 //                      cerr << "   labITS = " << labITS;
  85 //                      cerr << "   labTPC = " << labTPC;
  86 //                      cerr << endl;
  87
  88                         isGood = kFALSE;
  89                 }
  90                 if (nOK > 1) {
  91                         cerr << "More than 1 tracks with label " << labITS << ": taking the first" << endl;
  92                 }
  93                 true_px = *treeTrueTracks->GetV1();
  94                 true_py = *treeTrueTracks->GetV2();
  95                 true_pz = *treeTrueTracks->GetV3();
  96                 true_pt = TMath::Sqrt(true_px * true_px + true_py * true_py);
  97                 true_phi = TMath::ATan2(true_py, true_px);
  98                 if(true_phi < 0) true_phi += duepi;
  99                 true_phi *= 1000.0;
 100                 true_lambda = TMath::ATan(true_pz / true_pt) * 1000.0;
 101                 ptg = true_pt;
 102
 103                 // checks if two found good tracks have the same label in ITS
 104                 treeTrueTracks->Draw("entry", Form("label==%d && tpc_ok==1", abs(labITS)), "goff");
 105                 trueEntry = (Int_t)*treeTrueTracks->GetV1();
 106                 if (isGood && cnt[trueEntry] == 0)
 107                         cnt[trueEntry] = 1;
 108                 else if (isGood) {
 109                         cout << "Entry " << trueEntry << " already analyzed!" << endl;
 110                         continue;
 111                 }
 112
 113                 // charge matching
 114                 found_charge = (Int_t)iotrack->GetCharge();
 115                 ispos = isneg = kFALSE;
 116                 for (Int_t j = 0; j < npos; j++) ispos = (pdg_code == pos[j]);
 117                 for (Int_t j = 0; j < nneg; j++) isneg = (pdg_code == neg[j]);
 118                 if (ispos) countpos++;
 119                 if (isneg) countneg++;
 120
 121                 // pt resolution (%)
 122                 found_px = iotrack->GetPx();
 123                 found_py = iotrack->GetPy();
 124                 found_pz = iotrack->GetPz();
 125                 found_pt = TMath::Sqrt(found_px*found_px + found_py*found_py);
 126                 difpt = ((found_pt - true_pt) / true_pt) * 100.;
 127                 //cout << found_pt << " " << true_pt << " " << difpt << endl;
 128
 129                 // lambda (mrad)
 130                 found_tgl = iotrack->GetStateTgl();
 131                 found_lambda = TMath::ATan(found_tgl) * 1000.0;
 132                 diflambda = found_lambda - true_lambda;
 133 //              cout << "lambda " << found_lambda << " " << true_lambda << " " << diflambda << endl;
 134
 135                 // phi (mrad)
 136                 found_phi = TMath::ACos(found_px / found_pt);
 137                 if(found_phi > duepi) found_phi -= duepi;
 138                 if(found_phi < 0.) found_phi += duepi;
 139                 found_phi *= 1000.0;
 140                 difphi = found_phi - true_phi;
 141 //              cout << "phi " << found_phi << " " << true_phi << " " << difphi << endl;
 142
 143                 // impact parameters (microns)
 144                 Dr = iotrack->GetStateD() * 1.e4;
 145                 Dz = iotrack->GetDz() * 1.e4;
 146                 Dtot = sqrt(Dr*Dr + Dz*Dz);
 147
 148                 // curvature sign
 149                 found_curv = iotrack->GetStateC();
 150                 signC = (found_curv > 0.) ? 1 : -1;
 151
 152                 // fill tree
 153                 treeEvaluation->Fill();
 154         }
 155
 156
 157         TFile *outFile = new TFile("AliITSComparisonV1.root", "recreate");
 158         treeEvaluation->Write("Eval", TObject::kOverwrite);
 159         outFile->Close();
 160 }
 161