[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 <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"
	55	#include "AliTOFConstants.h"
	56	#include "AliTOFProb.h"
	57
	58	#include <stdlib.h>
f8014e68	59	#include <Riostream.h>
f8014e68	60	#include <Riostream.h>
72b127a8	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	//____________________________________________________________________________
5c016a7b	134	void AliTOFProb::Exec(const Option_t * /dummyOpt/)
72b127a8	135	{
	136	//
	137	// Performs Prob for TOF detector
	138	//
	139	fTask=1;
	140	fhfile->cd();
	141	// create the histo
c2c1e558	142	TH3* h3D=new TH3F("h3D", "", 120, -0.1, 1.1, 120, -0.1, 1.1, 120, -0.1, 1.1);
72b127a8	143	TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);
72b127a8	144
c2c1e558	145	TH1* hpion =new TH1F("hpion", "", 120, -0.1, 1.1);
	146	TH1* hkaon =new TH1F("hkaon", "", 120, -0.1, 1.1);
	147	TH1* hproton=new TH1F("hproton", "", 120, -0.1, 1.1);
72b127a8	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());
d8cf2f16	163	// unused for the time being in the context of probabilities
	164	// probably it will be included in AliESD approach
	165	//Float_t mass=fNtuple->GetLeaf("mext")->GetValue(0);
72b127a8	166	Int_t matc=(Int_t)(fNtuple->GetLeaf("matc")->GetValue(0));
	167	Int_t imam=(Int_t)(fNtuple->GetLeaf("imam")->GetValue(0));
	168	Float_t time=fNtuple->GetLeaf("text")->GetValue(0); // [ns]
	169	// generate a gaussian around time with sigma timeresolution [ns]
	170	Float_t lowTime=time-1000.*timeResolution;
	171	Float_t upTime =time+1000.*timeResolution;
	172	TF1* gauss=new TF1("tg","gaus",lowTime,upTime);
	173
	174	Double_t p1=time;
	175	Double_t p2=timeResolution;
be167bf4	176	//Float_t twopi=2*TMath::Pi();
be167bf4	177	//Float_t sqrtwopi=TMath::Sqrt(twopi);
72b127a8	178	//Double_t p0=1./(p2*sqrtwopi);
be167bf4	179	// choising a different normalization
72b127a8	180	Double_t p0=1.;
	181	gauss->SetParameters(p0,p1,p2);
	182
	183	Float_t px=fNtuple->GetLeaf("pxvtx")->GetValue(0);
	184	Float_t py=fNtuple->GetLeaf("pyvtx")->GetValue(0);
	185	Float_t pz=fNtuple->GetLeaf("pzvtx")->GetValue(0);
	186	// track length [cm]
	187	Float_t toflen=fNtuple->GetLeaf("leng")->GetValue(0);
	188
	189	Float_t pvtx=TMath::Sqrt(pxpx+pypy+pz*pz);
	190	Float_t ptvtx=TMath::Sqrt(pxpx+pypy);
	191
	192	Bool_t isSelected=(imam == 0 && pz !=0 && TMath::ATan(TMath::Abs(ptvtx/pz))>TMath::Pi()*45./180.);
	193	Bool_t trackWithTime=(matc==2 \|\| matc==3 \|\| matc==4);
	194
	195	Int_t abspdgcode=TMath::Abs(pdgcode);
	196	switch(abspdgcode){
	197	case 321:
	198	break;
	199	case 2212:
	200	break;
	201	case 11:
	202	break;
	203	default:
	204	break;
	205	}
	206
	207	if (isSelected && trackWithTime && toflen>370.)
	208	{//only primary +/-45
	209	Float_t beta[3]={0.,0.,0.};
	210	Float_t timeth[3]={0.,0.,0.};
	211	Float_t massarray[3]={0.13957,0.493677,0.9382723};
	212	for (Int_t j=0; j < 3; j++) {
	213	Float_t dummy=(pvtx/massarray[j])*(pvtx/massarray[j]);
	214	beta[j]=TMath::Sqrt(dummy/(1.+dummy));
	215	timeth[j]=toflen/(29.9792458*beta[j]);
	216	}
	217	Float_t pPion =(Float_t)gauss->Eval(timeth[0],0.,0.);
	218	Float_t pKaon =(Float_t)gauss->Eval(timeth[1],0.,0.);
	219	Float_t pProton =(Float_t)gauss->Eval(timeth[2],0.,0.);
	220
	221	Float_t dtpion=time-timeth[0];
	222	Float_t dtkaon=time-timeth[1];
	223	Float_t dtproton=time-timeth[2];
	224	hdtpion->Fill(dtpion);
	225	hdtkaon->Fill(dtkaon);
	226	hdtproton->Fill(dtproton);
	227
	228	// normalizing probabilities
	229	Float_t sumProb=pPion+pKaon+pProton;
	230	pPion=pPion/sumProb;
	231	pKaon=pKaon/sumProb;
	232	pProton=pProton/sumProb;
	233	hpion->Fill(pPion);
	234	hkaon->Fill(pKaon);
	235	hproton->Fill(pProton);
	236	//cout << "pPion " << pPion << " pKaon " << pKaon << " pProton " << pProton << endl;
	237	h3D->Fill(pPion,pKaon,pProton);
	238	h3Ddt->Fill(dtpion,dtkaon,dtproton);
	239	//cout << "pion " << timeth[0] << " kaon " << timeth[1] << " proton " << timeth[2] << endl;
	240	}// End of cuts appling
	241
	242	delete gauss;
	243
244	}// End of loop over particles
245
246	char outFileName[100];
247	strcpy(outFileName,"3DProbHisto.root");
248	TFile *houtfile = new TFile(outFileName,"recreate");
249	houtfile->cd();
250	h3D->Write(0,TObject::kOverwrite);
251	h3Ddt->Write(0,TObject::kOverwrite);
252	hpion->Write(0,TObject::kOverwrite);
253	hkaon->Write(0,TObject::kOverwrite);
254	hproton->Write(0,TObject::kOverwrite);
255	hdtpion->Write(0,TObject::kOverwrite);
256	hdtkaon->Write(0,TObject::kOverwrite);
257	hdtproton->Write(0,TObject::kOverwrite);
258
259	houtfile->Close();
260	houtfile->Write(0,TObject::kOverwrite);
261	}
262
263
264	//__________________________________________________________________
5c016a7b	265	Bool_t AliTOFProb::operator==( AliTOFProb const & /tofrec/)const
72b127a8	266	{
	267	// dummy version of Equal operator.
	268	// requested by coding conventions
	269	return kTRUE;
	270
	271	}