AliTOFProb and AliTOFtestProb.C added
[u/mrichter/AliRoot.git] / TOF / AliTOFProb.cxx
CommitLineData
72b127a8 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//_________________________________________________________________________
17// TTask class for TOF Probabilities (ALICE-Offline week June 2002).
18// Use case: start root and execute the following macro
19// ntuple.root is assumed to be the ntuple with the reconstructed
20// data
21/*
22{
23 // Dynamically link some shared libs
24 if (gClassTable->GetID("AliRun") < 0) {
25 gROOT->LoadMacro("loadlibs.C");
26 loadlibs();
27 }
28 // create an instance of the AliTOFProb class
29 // You have to pass the ntuple
30 AliTOFProb* tofprob=new AliTOFProb("ntuple.root");
31 tofprob->Exec("");
32}
33*/
34//
35//-- Author: F. Pierella | pierella@bo.infn.it
36//////////////////////////////////////////////////////////////////////////////
37
38#include "TStyle.h"
39#include "TTask.h"
40#include "TTree.h"
41#include "TSystem.h"
42#include "TFile.h"
43#include "TH3.h"
44//#include "TH3.h"
45#include "TCanvas.h"
46#include "TFile.h"
47#include <TF1.h>
48#include <TF2.h>
49#include "TROOT.h"
50#include "TFolder.h"
51#include "TNtuple.h"
52#include "TLeaf.h"
53
54#include "AliConst.h"
55#include "AliTOFConstants.h"
56#include "AliTOFProb.h"
57
58#include <stdlib.h>
59#include <iostream.h>
60#include <fstream.h>
61
62ClassImp(AliTOFProb)
63
64//____________________________________________________________________________
65 AliTOFProb::AliTOFProb():TTask("AliTOFProb","")
66{
67 // default ctor - set the pointer member vars to zero
68
69 fhfile = 0;
70 fNtuple = 0;
71 fgen = 0;
72 foutfileName = 0;
73}
74
75//____________________________________________________________________________
76 AliTOFProb::AliTOFProb(char* headerFile):TTask("AliTOFProb","")
77{
78
79 fhfile = TFile::Open(headerFile); // connect file with ntuple
80 foutfileName=headerFile;
81 fNtuple = 0;
82 fgen = 0;
83
84 Init();
85 // add Task to //root/Tasks folder
86 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
87 roottasks->Add(this) ;
88}
89//____________________________________________________________________________
90void AliTOFProb::Init()
91{
92
93 fNtuple= (TNtuple*)fhfile->Get("Ntuple"); // get ntuple from file
94 Int_t nvar = fNtuple->GetNvar(); cout <<"N of var.="<< nvar << endl;
95 fNtuple->GetEvent(0);
96
97}
98
99//____________________________________________________________________________
100 AliTOFProb::~AliTOFProb()
101{
102 //
103 // dtor (free used memory)
104 //
105
106 if (fhfile)
107 {
108 delete fhfile;
109 fhfile = 0;
110 }
111
112
113 if (fNtuple)
114 {
115 delete fNtuple;
116 fNtuple = 0;
117 }
118
119 if (fgen)
120 {
121 delete fgen;
122 fgen = 0;
123 }
124
125 if (foutfileName)
126 {
127 delete foutfileName;
128 foutfileName = 0;
129 }
130}
131
132
133//____________________________________________________________________________
134void AliTOFProb::Exec(const Option_t *dummyOpt)
135{
136 //
137 // Performs Prob for TOF detector
138 //
139 fTask=1;
140 fhfile->cd();
141 // create the histo
142 TH3* h3D=new TH3F("h3D", "", 100, 0., 1., 100, 0., 1., 100, 0., 1.);
143 TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);
144
145 TH1* hpion =new TH1F("hpion", "", 100, 0., 1.);
146 TH1* hkaon =new TH1F("hkaon", "", 100, 0., 1.);
147 TH1* hproton=new TH1F("hproton", "", 100, 0., 1.);
148
149 TH1* hdtpion =new TH1F("hdtpion", "", 600, -120., 120.);
150 TH1* hdtkaon =new TH1F("hdtkaon", "", 600, -120., 120.);
151 TH1* hdtproton=new TH1F("hdtproton", "", 600, -120., 120.);
152
153 const Float_t timeResolution=0.12;
154
155 Int_t nparticles = (Int_t)fNtuple->GetEntries();
156 cout << " Number of nparticles =" << nparticles << endl;
157 if (nparticles <= 0) return;
158
159 for (Int_t i=0; i < nparticles; i++) {
160 fNtuple->GetEvent(i);
161 //Int_t event=(Int_t)(fNtuple->GetLeaf("event")->GetValue());
162 Int_t pdgcode=(Int_t)(fNtuple->GetLeaf("ipart")->GetValue());
163 Float_t mass=fNtuple->GetLeaf("mext")->GetValue(0);
164 Int_t matc=(Int_t)(fNtuple->GetLeaf("matc")->GetValue(0));
165 Int_t imam=(Int_t)(fNtuple->GetLeaf("imam")->GetValue(0));
166 Float_t time=fNtuple->GetLeaf("text")->GetValue(0); // [ns]
167 // generate a gaussian around time with sigma timeresolution [ns]
168 Float_t lowTime=time-1000.*timeResolution;
169 Float_t upTime =time+1000.*timeResolution;
170 TF1* gauss=new TF1("tg","gaus",lowTime,upTime);
171
172 Double_t p1=time;
173 Double_t p2=timeResolution;
174 Float_t twopi=2*TMath::Pi();
175 Float_t sqrtwopi=TMath::Sqrt(twopi);
176 //Double_t p0=1./(p2*sqrtwopi);
177 Double_t p0=1.;
178 gauss->SetParameters(p0,p1,p2);
179
180 Float_t px=fNtuple->GetLeaf("pxvtx")->GetValue(0);
181 Float_t py=fNtuple->GetLeaf("pyvtx")->GetValue(0);
182 Float_t pz=fNtuple->GetLeaf("pzvtx")->GetValue(0);
183 // track length [cm]
184 Float_t toflen=fNtuple->GetLeaf("leng")->GetValue(0);
185
186 Float_t pvtx=TMath::Sqrt(px*px+py*py+pz*pz);
187 Float_t ptvtx=TMath::Sqrt(px*px+py*py);
188
189 Bool_t isSelected=(imam == 0 && pz !=0 && TMath::ATan(TMath::Abs(ptvtx/pz))>TMath::Pi()*45./180.);
190 Bool_t trackWithTime=(matc==2 || matc==3 || matc==4);
191
192 Int_t abspdgcode=TMath::Abs(pdgcode);
193 switch(abspdgcode){
194 case 321:
195 break;
196 case 2212:
197 break;
198 case 11:
199 break;
200 default:
201 break;
202 }
203
204 if (isSelected && trackWithTime && toflen>370.)
205 {//only primary +/-45
206 Float_t beta[3]={0.,0.,0.};
207 Float_t timeth[3]={0.,0.,0.};
208 Float_t massarray[3]={0.13957,0.493677,0.9382723};
209 for (Int_t j=0; j < 3; j++) {
210 Float_t dummy=(pvtx/massarray[j])*(pvtx/massarray[j]);
211 beta[j]=TMath::Sqrt(dummy/(1.+dummy));
212 timeth[j]=toflen/(29.9792458*beta[j]);
213 }
214 Float_t pPion =(Float_t)gauss->Eval(timeth[0],0.,0.);
215 Float_t pKaon =(Float_t)gauss->Eval(timeth[1],0.,0.);
216 Float_t pProton =(Float_t)gauss->Eval(timeth[2],0.,0.);
217
218 Float_t dtpion=time-timeth[0];
219 Float_t dtkaon=time-timeth[1];
220 Float_t dtproton=time-timeth[2];
221 hdtpion->Fill(dtpion);
222 hdtkaon->Fill(dtkaon);
223 hdtproton->Fill(dtproton);
224
225 // normalizing probabilities
226 Float_t sumProb=pPion+pKaon+pProton;
227 pPion=pPion/sumProb;
228 pKaon=pKaon/sumProb;
229 pProton=pProton/sumProb;
230 hpion->Fill(pPion);
231 hkaon->Fill(pKaon);
232 hproton->Fill(pProton);
233 //cout << "pPion " << pPion << " pKaon " << pKaon << " pProton " << pProton << endl;
234 h3D->Fill(pPion,pKaon,pProton);
235 h3Ddt->Fill(dtpion,dtkaon,dtproton);
236 //cout << "pion " << timeth[0] << " kaon " << timeth[1] << " proton " << timeth[2] << endl;
237 }// End of cuts appling
238
239 delete gauss;
240
241 }// End of loop over particles
242
243 char outFileName[100];
244 strcpy(outFileName,"3DProbHisto.root");
245 TFile *houtfile = new TFile(outFileName,"recreate");
246 houtfile->cd();
247 h3D->Write(0,TObject::kOverwrite);
248 h3Ddt->Write(0,TObject::kOverwrite);
249 hpion->Write(0,TObject::kOverwrite);
250 hkaon->Write(0,TObject::kOverwrite);
251 hproton->Write(0,TObject::kOverwrite);
252 hdtpion->Write(0,TObject::kOverwrite);
253 hdtkaon->Write(0,TObject::kOverwrite);
254 hdtproton->Write(0,TObject::kOverwrite);
255
256 houtfile->Close();
257 houtfile->Write(0,TObject::kOverwrite);
258}
259
260
261//__________________________________________________________________
262Bool_t AliTOFProb::operator==( AliTOFProb const & tofrec)const
263{
264 // dummy version of Equal operator.
265 // requested by coding conventions
266 return kTRUE;
267
268}