uuhhhhaaa - what did I change? Added the run number and fixed coding conventions ;-)
[u/mrichter/AliRoot.git] / TOF / AliTOFProb.cxx
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 "AliTOFProb.h"
56
57 #include <stdlib.h>
58 #include <Riostream.h>
59 #include <Riostream.h>
60
61 ClassImp(AliTOFProb)
62
63 //____________________________________________________________________________ 
64   AliTOFProb::AliTOFProb():TTask("AliTOFProb","") 
65 {
66   // default ctor - set the pointer member vars to zero
67
68   fhfile    = 0;
69   fNtuple   = 0;
70   fgen      = 0;
71   foutfileName  = 0;
72 }
73            
74 //____________________________________________________________________________ 
75   AliTOFProb::AliTOFProb(char* headerFile):TTask("AliTOFProb","") 
76 {
77
78   fhfile = TFile::Open(headerFile); // connect file with ntuple
79   foutfileName=headerFile;
80   fNtuple   = 0;
81   fgen      = 0;
82
83   Init();
84   // add Task to //root/Tasks folder
85   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
86   roottasks->Add(this) ; 
87 }
88 //____________________________________________________________________________ 
89 void AliTOFProb::Init()
90 {
91
92   fNtuple= (TNtuple*)fhfile->Get("Ntuple"); // get ntuple from file
93   Int_t nvar = fNtuple->GetNvar(); cout <<"N of var.="<< nvar << endl;
94   fNtuple->GetEvent(0);
95
96 }
97
98 //____________________________________________________________________________ 
99   AliTOFProb::~AliTOFProb()
100 {
101   //
102   // dtor (free used memory)
103   //
104
105   if (fhfile)
106     {
107       delete fhfile;
108       fhfile = 0;
109     }
110
111
112   if (fNtuple)
113     {
114       delete fNtuple;
115       fNtuple = 0;
116     }
117
118   if (fgen)
119     {
120       delete fgen;
121       fgen = 0;
122     }
123
124   if (foutfileName)
125     {
126       delete foutfileName;
127       foutfileName = 0;
128     }
129 }
130
131
132 //____________________________________________________________________________
133 void AliTOFProb::Exec(const Option_t * /*dummyOpt*/) 
134
135   //
136   // Performs Prob for TOF detector
137   // 
138   fTask=1;
139   fhfile->cd();
140   // create the histo
141   TH3* h3D=new TH3F("h3D", "", 120, -0.1, 1.1, 120, -0.1, 1.1, 120, -0.1, 1.1);
142   TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);
143
144   TH1* hpion  =new TH1F("hpion", "", 120, -0.1, 1.1);
145   TH1* hkaon  =new TH1F("hkaon", "", 120, -0.1, 1.1);
146   TH1* hproton=new TH1F("hproton", "", 120, -0.1, 1.1);
147
148   TH1* hdtpion  =new TH1F("hdtpion", "", 600, -120., 120.);
149   TH1* hdtkaon  =new TH1F("hdtkaon", "", 600, -120., 120.);
150   TH1* hdtproton=new TH1F("hdtproton", "", 600, -120., 120.);
151
152   const Float_t timeResolution=0.12;
153
154   Int_t nparticles = (Int_t)fNtuple->GetEntries();
155   cout << " Number of nparticles =" << nparticles << endl;
156   if (nparticles <= 0) return;
157   
158   for (Int_t i=0; i < nparticles; i++) {
159     fNtuple->GetEvent(i);
160     //Int_t event=(Int_t)(fNtuple->GetLeaf("event")->GetValue());
161     Int_t pdgcode=(Int_t)(fNtuple->GetLeaf("ipart")->GetValue());
162     // unused for the time being in the context of probabilities
163     // probably it will be included in AliESD approach
164     //Float_t mass=fNtuple->GetLeaf("mext")->GetValue(0);
165     Int_t matc=(Int_t)(fNtuple->GetLeaf("matc")->GetValue(0));
166     Int_t imam=(Int_t)(fNtuple->GetLeaf("imam")->GetValue(0));
167     Float_t time=fNtuple->GetLeaf("text")->GetValue(0); // [ns]
168     // generate a gaussian around time with sigma timeresolution [ns]
169     Float_t lowTime=time-1000.*timeResolution;
170     Float_t upTime =time+1000.*timeResolution;
171     TF1* gauss=new TF1("tg","gaus",lowTime,upTime);
172
173     Double_t p1=time;
174     Double_t p2=timeResolution;
175     //Float_t twopi=2*TMath::Pi();
176     //Float_t sqrtwopi=TMath::Sqrt(twopi);
177     //Double_t p0=1./(p2*sqrtwopi);
178     // choising a different normalization
179     Double_t p0=1.;
180     gauss->SetParameters(p0,p1,p2);
181
182     Float_t px=fNtuple->GetLeaf("pxvtx")->GetValue(0);
183     Float_t py=fNtuple->GetLeaf("pyvtx")->GetValue(0);
184     Float_t pz=fNtuple->GetLeaf("pzvtx")->GetValue(0);
185     // track length [cm]
186     Float_t toflen=fNtuple->GetLeaf("leng")->GetValue(0);
187
188     Float_t pvtx=TMath::Sqrt(px*px+py*py+pz*pz);
189     Float_t ptvtx=TMath::Sqrt(px*px+py*py);
190
191     Bool_t isSelected=(imam == 0 && pz !=0 && TMath::ATan(TMath::Abs(ptvtx/pz))>TMath::Pi()*45./180.);
192     Bool_t trackWithTime=(matc==2 || matc==3 || matc==4);
193
194     Int_t abspdgcode=TMath::Abs(pdgcode);
195     switch(abspdgcode){
196     case 321:
197       break;
198     case 2212:
199       break;
200     case 11:
201       break;
202     default:
203       break;
204     }
205
206     if (isSelected && trackWithTime && toflen>370.)
207       {//only primary +/-45  
208         Float_t beta[3]={0.,0.,0.};
209         Float_t timeth[3]={0.,0.,0.};
210         Float_t massarray[3]={0.13957,0.493677,0.9382723};
211         for (Int_t j=0; j < 3; j++) {
212           Float_t dummy=(pvtx/massarray[j])*(pvtx/massarray[j]);
213           beta[j]=TMath::Sqrt(dummy/(1.+dummy));
214           timeth[j]=toflen/(29.9792458*beta[j]);
215         }
216         Float_t pPion   =(Float_t)gauss->Eval(timeth[0],0.,0.);
217         Float_t pKaon   =(Float_t)gauss->Eval(timeth[1],0.,0.);
218         Float_t pProton =(Float_t)gauss->Eval(timeth[2],0.,0.);
219
220         Float_t dtpion=time-timeth[0];
221         Float_t dtkaon=time-timeth[1];
222         Float_t dtproton=time-timeth[2];
223         hdtpion->Fill(dtpion);
224         hdtkaon->Fill(dtkaon);
225         hdtproton->Fill(dtproton);
226
227         // normalizing probabilities
228         Float_t sumProb=pPion+pKaon+pProton;
229         pPion=pPion/sumProb;
230         pKaon=pKaon/sumProb;
231         pProton=pProton/sumProb;
232         hpion->Fill(pPion);
233         hkaon->Fill(pKaon);
234         hproton->Fill(pProton);
235         //cout << "pPion " << pPion << " pKaon " << pKaon << " pProton " << pProton << endl;
236         h3D->Fill(pPion,pKaon,pProton);
237         h3Ddt->Fill(dtpion,dtkaon,dtproton);
238         //cout << "pion " << timeth[0] << " kaon " << timeth[1] << " proton " << timeth[2] << endl;
239       }// End of cuts appling
240
241     delete gauss;
242
243   }// End of loop over particles
244
245   char outFileName[100];
246   strcpy(outFileName,"3DProbHisto.root");
247   TFile *houtfile = new TFile(outFileName,"recreate");
248   houtfile->cd();
249   h3D->Write(0,TObject::kOverwrite);
250   h3Ddt->Write(0,TObject::kOverwrite);
251   hpion->Write(0,TObject::kOverwrite);
252   hkaon->Write(0,TObject::kOverwrite);
253   hproton->Write(0,TObject::kOverwrite);
254   hdtpion->Write(0,TObject::kOverwrite);
255   hdtkaon->Write(0,TObject::kOverwrite);
256   hdtproton->Write(0,TObject::kOverwrite);
257
258   houtfile->Close();
259   houtfile->Write(0,TObject::kOverwrite);
260 }
261
262
263 //__________________________________________________________________
264 Bool_t AliTOFProb::operator==( AliTOFProb const & /*tofrec*/)const
265 {
266   // dummy version of Equal operator.
267   // requested by coding conventions
268   return kTRUE;
269
270 }