[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 "AliTOFProb.h"

#include <stdlib.h>
#include <Riostream.h>
#include <Riostream.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", "", 120, -0.1, 1.1, 120, -0.1, 1.1, 120, -0.1, 1.1);
  TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);

  TH1* hpion  =new TH1F("hpion", "", 120, -0.1, 1.1);
  TH1* hkaon  =new TH1F("hkaon", "", 120, -0.1, 1.1);
  TH1* hproton=new TH1F("hproton", "", 120, -0.1, 1.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());
    // unused for the time being in the context of probabilities
    // probably it will be included in AliESD approach
    //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);
    // choising a different normalization
    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"
72b127a8	55	#include "AliTOFProb.h"
	56
	57	#include <stdlib.h>
f8014e68	58	#include <Riostream.h>
f8014e68	59	#include <Riostream.h>
72b127a8	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	//____________________________________________________________________________
5c016a7b	133	void AliTOFProb::Exec(const Option_t * /dummyOpt/)
72b127a8	134	{
	135	//
	136	// Performs Prob for TOF detector
	137	//
	138	fTask=1;
	139	fhfile->cd();
	140	// create the histo
c2c1e558	141	TH3* h3D=new TH3F("h3D", "", 120, -0.1, 1.1, 120, -0.1, 1.1, 120, -0.1, 1.1);
72b127a8	142	TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);
72b127a8	143
c2c1e558	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);
72b127a8	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());
d8cf2f16	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);
72b127a8	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;
be167bf4	175	//Float_t twopi=2*TMath::Pi();
be167bf4	176	//Float_t sqrtwopi=TMath::Sqrt(twopi);
72b127a8	177	//Double_t p0=1./(p2*sqrtwopi);
be167bf4	178	// choising a different normalization
72b127a8	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(pxpx+pypy+pz*pz);
	189	Float_t ptvtx=TMath::Sqrt(pxpx+pypy);
	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	//__________________________________________________________________
5c016a7b	264	Bool_t AliTOFProb::operator==( AliTOFProb const & /tofrec/)const
72b127a8	265	{
	266	// dummy version of Equal operator.
	267	// requested by coding conventions
	268	return kTRUE;
	269
	270	}