macros/CompareTimeHypothesis.C

   1 //
   2 // This macro draws time resolution for true and false particle hypothesis
   3 // on a given reference plane in a given momentum range.
   4 // True hypothesis is painted in BLUE while false in RED
   5 //
   6 // Input parameters:
   7 // nRefPlane - indicate the reference plane.
   8 //   1: in TPC
   9 //   2: out TPC
  10 //   3: out TRD
  11 //
  12 // minP, maxP - total momentum cuts. momentum from gAlice are used for cuts
  13 //
  14 // The macro is using SORTED reference points in a file "trackRef.root".
  15 // Track references can be sorted using CopyReferences.C
  16 //
  17 // Before using check filenames and names of trees invlolved
  18 //
  19 // ITS tracks (nRefPlane == 1) are compared at inner TPC ref Plane
  20 // make shure they are properly propagated to this plane.
  21 //
  22 // Sylwester Radomski, e-mail: S.Radomski@gsi.de
  23 // Feb 14, 2002
  24 //
  25
  26 CompareTimeHypothesis (Int_t nRefPlane, Double_t minP, Double_t maxP) {
  27
  28   // Check input data Consistancy
  29
  30   if (nRefPlane < 1 || nRefPlane > 3) {
  31     cout << "Wrong Reference Plane Id ( " << nRefPalne << ") [1-3] " << endl;
  32     return;
  33   }
  34
  35   if (maxP < minP) {
  36     cout << "MaxP lesser than MinP" << endl;
  37     return;
  38   }
  39
  40   // Styling
  41   gROOT->SetStyle("Plain");
  42   gStyle->SetOptStat(1111);
  43   gStyle->SetStatH(0.2);
  44   gStyle->SetStatW(0.3);
  45
  46   // magnetic field 0.4 T
  47   Double_t cc = 0.299792458;
  48   AliKalmanTrack::SetConvConst(100/cc/0.4);
  49
  50   const char* gFileName = "galice.root";
  51   const char* refFileName = "trackRef.root";
  52
  53   const char* trackFileName[] = {"AliTPCtracks.root", "AliTPCBackTracks.root", "AliTRDtracks.root"};
  54   const char *treeName[] = {"TreeT_ITSb_0", "seedsTPCtoTRD_0", "TRDb_0"};
  55   const char *refTreeName[] = {"TPC", "TPC", "TRD"};
  56
  57   // Histograms
  58
  59   const Int_t nTypes = 5;
  60   Int_t codes[] = {11, 13, 211, 321, 2212};
  61
  62   //Double_t cutoff[] = {20, 20, 20, 100, 100};
  63   Double_t cutoff[] = {30, 30, 30, 100, 200};
  64
  65   const char *names[3*nTypes] = {
  66     "Electrons", "Muons", "Pions", "Kaons", "Protons",
  67     "11s", "13s", "211s", "321s", "2212s",
  68     "11o", "13o", "211o", "321o", "2212o"
  69   };
  70
  71   THStack *stack[nTypes];
  72   TH1D *hDiffSame[nTypes], *hDiffOdd[nTypes];
  73
  74   for(Int_t i=0; i<nTypes; i++) {
  75
  76     stack[i] = new THStack(names[i], names[i]);
  77     hDiffSame[i] = new TH1D(names[nTypes+i], names[nTypes+i], 50, -cutoff[i], cutoff[i]);
  78     hDiffOdd[i] = new TH1D(names[2*nTypes+i], names[2*nTypes+i], 50, -cutoff[i], cutoff[i]);
  79
  80     hDiffSame[i]->SetFillColor(kBlue);
  81     hDiffOdd[i]->SetFillColor(kRed);
  82
  83     //hDiffSame[i]->SetFillStyle(3003);
  84     //hDiffOdd[i]->SetFillStyle(3008);
  85
  86     stack[i]->Add(hDiffOdd[i]);
  87     stack[i]->Add(hDiffSame[i]);
  88   }
  89
  90   TH1D *hLength = new TH1D("length", "Length Resolution", 100, -0.5, 0.5);
  91
  92   // Particles
  93
  94   TFile *refFile = new TFile(gFileName, "READ");
  95   gAlice = (AliRun*)refFile->Get("gAlice");
  96   gAlice->GetEvent(0);
  97
  98
  99   // Reference tracks
 100
 101   refFile = new TFile(refFileName, "READ");
 102   TTree *treeRef = (TTree*)refFile->Get("TreeTR0_Sorted");
 103   TBranch *trkRef = treeRef->GetBranch(refTreeName[nRefPlane-1]);
 104   TClonesArray *trkRefs = new TClonesArray("AliTrackReference", 100);
 105   trkRef->SetAddress(&trkRefs);
 106
 107   Int_t ntracks = trkRef->GetEntries();
 108   cout << ntracks << endl;
 109
 110
 111   // Tracks
 112   TFile *trackFile = new TFile(trackFileName[nRefPlane-1]);
 113   TTree *tracks = (TTree*)trackFile->Get(treeName[nRefPlane-1]);
 114   TBranch *trkTracks = tracks->GetBranch("tracks");
 115   AliTPCtrack *track = 0;
 116   trkTracks->SetAddress(&track);
 117
 118   for(Int_t t=0; t<tracks->GetEntries(); t++) {
 119
 120     cout << t << "\r";
 121
 122     trkTracks->GetEvent(t);
 123
 124     Int_t lab = track->GetLabel();
 125     if (lab < 0 || lab > 100000) continue;
 126
 127     Int_t pdg = gAlice->Particle(lab)->GetPdgCode();
 128     treeRef->GetEvent(lab);
 129
 130     Int_t nEntrTRK = trkRefs->GetEntries();
 131
 132     // Primaries within momentum range
 133     if (gAlice->Particle(lab)->P() < minP ) continue;
 134     if (gAlice->Particle(lab)->P() > maxP ) continue;
 135     if (gAlice->Particle(lab)->GetFirstMother() != -1) continue;
 136     if (!(nEntrTRK)) continue;
 137
 138     Double_t timeTrue, timeTrack;
 139     Double_t trueLength;
 140
 141     Int_t index =  nEntrTRK-1;
 142     if (nRefPlane == 1) index = 0;
 143
 144     timeTrue = ((AliTrackReference*)(*trkRefs)[index])->GetTime() * 1e12;
 145
 146     //timeTrue = gRandom->Gaus(timeTrue, 10);
 147
 148     trueLength = ((AliTrackReference*)(*trkRefs)[index])->GetLength();
 149     hLength->Fill(trueLength - track->GetIntegratedLength());
 150
 151     for(Int_t i=0; i<nTypes; i++) {
 152
 153       timeTrack = track->GetIntegratedTime(codes[i]);
 154       if (abs(pdg) == codes[i]) hDiffSame[i]->Fill(timeTrue - timeTrack);
 155       else hDiffOdd[i]->Fill(timeTrue - timeTrack);
 156     }
 157   }
 158
 159   TCanvas *c = new TCanvas();
 160   c->SetWindowSize(640, 800);
 161   c->Divide(2,3);
 162
 163   c->cd(1);
 164   hLength->Draw();
 165   hLength->SetXTitle("True Length - Measured Length [cm]");
 166
 167   for(Int_t i=0; i<nTypes; i++) {
 168     c->cd(2+i);
 169     gPad->SetGridx();
 170     gPad->SetGridy();
 171     stack[i]->Draw();
 172     stack[i]->GetXaxis()->SetTitle("time true - measured [ps]");
 173   }
 174
 175   if (treeRef) delete treeRef;
 176   if (tracks) delete tracks;
 177 }