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 }