[u/mrichter/AliRoot.git] / TOF / AliTOFProb.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//_________________________________________________________________________
// TTask class for TOF Probabilities (ALICE-Offline week June 2002).
// Use case: start root and execute the following macro
// ntuple.root is assumed to be the ntuple with the reconstructed
// data
/*
{
  // Dynamically link some shared libs
  if (gClassTable->GetID("AliRun") < 0) {
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
  }
  // create an instance of the AliTOFProb class
  // You have to pass the ntuple
 AliTOFProb* tofprob=new AliTOFProb("ntuple.root");
 tofprob->Exec("");
}
*/
// 
//-- Author: F. Pierella | pierella@bo.infn.it
//////////////////////////////////////////////////////////////////////////////

#include "TStyle.h"
#include "TTask.h"
#include "TTree.h"
#include "TSystem.h"
#include "TFile.h"
#include "TH3.h"
//#include "TH3.h"
#include "TCanvas.h"
#include "TFile.h"
#include <TF1.h>
#include <TF2.h>
#include "TROOT.h"
#include "TFolder.h"
#include "TNtuple.h"
#include "TLeaf.h"

#include "AliConst.h"
#include "AliTOFConstants.h"
#include "AliTOFProb.h"

#include <stdlib.h>
#include <iostream.h>
#include <fstream.h>

ClassImp(AliTOFProb)

//____________________________________________________________________________ 
  AliTOFProb::AliTOFProb():TTask("AliTOFProb","") 
{
  // default ctor - set the pointer member vars to zero

  fhfile    = 0;
  fNtuple   = 0;
  fgen      = 0;
  foutfileName  = 0;
}
           
//____________________________________________________________________________ 
  AliTOFProb::AliTOFProb(char* headerFile):TTask("AliTOFProb","") 
{

  fhfile = TFile::Open(headerFile); // connect file with ntuple
  foutfileName=headerFile;
  fNtuple   = 0;
  fgen      = 0;

  Init();
  // add Task to //root/Tasks folder
  TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
  roottasks->Add(this) ; 
}
//____________________________________________________________________________ 
void AliTOFProb::Init()
{

  fNtuple= (TNtuple*)fhfile->Get("Ntuple"); // get ntuple from file
  Int_t nvar = fNtuple->GetNvar(); cout <<"N of var.="<< nvar << endl;
  fNtuple->GetEvent(0);

}

//____________________________________________________________________________ 
  AliTOFProb::~AliTOFProb()
{
  //
  // dtor (free used memory)
  //

  if (fhfile)
    {
      delete fhfile;
      fhfile = 0;
    }


  if (fNtuple)
    {
      delete fNtuple;
      fNtuple = 0;
    }

  if (fgen)
    {
      delete fgen;
      fgen = 0;
    }

  if (foutfileName)
    {
      delete foutfileName;
      foutfileName = 0;
    }
}


//____________________________________________________________________________
void AliTOFProb::Exec(const Option_t *dummyOpt) 
{ 
  //
  // Performs Prob for TOF detector
  // 
  fTask=1;
  fhfile->cd();
  // create the histo
  TH3* h3D=new TH3F("h3D", "", 100, 0., 1., 100, 0., 1., 100, 0., 1.);
  TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);

  TH1* hpion  =new TH1F("hpion", "", 100, 0., 1.);
  TH1* hkaon  =new TH1F("hkaon", "", 100, 0., 1.);
  TH1* hproton=new TH1F("hproton", "", 100, 0., 1.);

  TH1* hdtpion  =new TH1F("hdtpion", "", 600, -120., 120.);
  TH1* hdtkaon  =new TH1F("hdtkaon", "", 600, -120., 120.);
  TH1* hdtproton=new TH1F("hdtproton", "", 600, -120., 120.);

  const Float_t timeResolution=0.12;

  Int_t nparticles = (Int_t)fNtuple->GetEntries();
  cout << " Number of nparticles =" << nparticles << endl;
  if (nparticles <= 0) return;
  
  for (Int_t i=0; i < nparticles; i++) {
    fNtuple->GetEvent(i);
    //Int_t event=(Int_t)(fNtuple->GetLeaf("event")->GetValue());
    Int_t pdgcode=(Int_t)(fNtuple->GetLeaf("ipart")->GetValue());
    Float_t mass=fNtuple->GetLeaf("mext")->GetValue(0);
    Int_t matc=(Int_t)(fNtuple->GetLeaf("matc")->GetValue(0));
    Int_t imam=(Int_t)(fNtuple->GetLeaf("imam")->GetValue(0));
    Float_t time=fNtuple->GetLeaf("text")->GetValue(0); // [ns]
    // generate a gaussian around time with sigma timeresolution [ns]
    Float_t lowTime=time-1000.*timeResolution;
    Float_t upTime =time+1000.*timeResolution;
    TF1* gauss=new TF1("tg","gaus",lowTime,upTime);

    Double_t p1=time;
    Double_t p2=timeResolution;
    Float_t twopi=2*TMath::Pi();
    Float_t sqrtwopi=TMath::Sqrt(twopi);
    //Double_t p0=1./(p2*sqrtwopi);
    Double_t p0=1.;
    gauss->SetParameters(p0,p1,p2);

    Float_t px=fNtuple->GetLeaf("pxvtx")->GetValue(0);
    Float_t py=fNtuple->GetLeaf("pyvtx")->GetValue(0);
    Float_t pz=fNtuple->GetLeaf("pzvtx")->GetValue(0);
    // track length [cm]
    Float_t toflen=fNtuple->GetLeaf("leng")->GetValue(0);

    Float_t pvtx=TMath::Sqrt(px*px+py*py+pz*pz);
    Float_t ptvtx=TMath::Sqrt(px*px+py*py);

    Bool_t isSelected=(imam == 0 && pz !=0 && TMath::ATan(TMath::Abs(ptvtx/pz))>TMath::Pi()*45./180.);
    Bool_t trackWithTime=(matc==2 || matc==3 || matc==4);

    Int_t abspdgcode=TMath::Abs(pdgcode);
    switch(abspdgcode){
    case 321:
      break;
    case 2212:
      break;
    case 11:
      break;
    default:
      break;
    }

    if (isSelected && trackWithTime && toflen>370.)
      {//only primary +/-45  
	Float_t beta[3]={0.,0.,0.};
	Float_t timeth[3]={0.,0.,0.};
	Float_t massarray[3]={0.13957,0.493677,0.9382723};
	for (Int_t j=0; j < 3; j++) {
	  Float_t dummy=(pvtx/massarray[j])*(pvtx/massarray[j]);
	  beta[j]=TMath::Sqrt(dummy/(1.+dummy));
	  timeth[j]=toflen/(29.9792458*beta[j]);
	}
	Float_t pPion   =(Float_t)gauss->Eval(timeth[0],0.,0.);
	Float_t pKaon   =(Float_t)gauss->Eval(timeth[1],0.,0.);
	Float_t pProton =(Float_t)gauss->Eval(timeth[2],0.,0.);

	Float_t dtpion=time-timeth[0];
	Float_t dtkaon=time-timeth[1];
	Float_t dtproton=time-timeth[2];
	hdtpion->Fill(dtpion);
	hdtkaon->Fill(dtkaon);
	hdtproton->Fill(dtproton);

	// normalizing probabilities
	Float_t sumProb=pPion+pKaon+pProton;
	pPion=pPion/sumProb;
	pKaon=pKaon/sumProb;
	pProton=pProton/sumProb;
	hpion->Fill(pPion);
	hkaon->Fill(pKaon);
	hproton->Fill(pProton);
	//cout << "pPion " << pPion << " pKaon " << pKaon << " pProton " << pProton << endl;
	h3D->Fill(pPion,pKaon,pProton);
	h3Ddt->Fill(dtpion,dtkaon,dtproton);
	//cout << "pion " << timeth[0] << " kaon " << timeth[1] << " proton " << timeth[2] << endl;
      }// End of cuts appling

    delete gauss;

  }// End of loop over particles

  char outFileName[100];
  strcpy(outFileName,"3DProbHisto.root");
  TFile *houtfile = new TFile(outFileName,"recreate");
  houtfile->cd();
  h3D->Write(0,TObject::kOverwrite);
  h3Ddt->Write(0,TObject::kOverwrite);
  hpion->Write(0,TObject::kOverwrite);
  hkaon->Write(0,TObject::kOverwrite);
  hproton->Write(0,TObject::kOverwrite);
  hdtpion->Write(0,TObject::kOverwrite);
  hdtkaon->Write(0,TObject::kOverwrite);
  hdtproton->Write(0,TObject::kOverwrite);

  houtfile->Close();
  houtfile->Write(0,TObject::kOverwrite);
}


//__________________________________________________________________
Bool_t AliTOFProb::operator==( AliTOFProb const & tofrec)const
{
  // dummy version of Equal operator.
  // requested by coding conventions
  return kTRUE;

}
Commit	Line	Data
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
	62	ClassImp(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	//____________________________________________________________________________
90	void 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	//____________________________________________________________________________
134	void 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(pxpx+pypy+pz*pz);
187	Float_t ptvtx=TMath::Sqrt(pxpx+pypy);
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	//__________________________________________________________________
262	Bool_t AliTOFProb::operator==( AliTOFProb const & tofrec)const
263	{
264	// dummy version of Equal operator.
265	// requested by coding conventions
266	return kTRUE;
267
268	}