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 "AliTOFConstants.h"
  56 #include "AliTOFProb.h"
  57
  58 #include <stdlib.h>
  59 #include <Riostream.h>
  60 #include <Riostream.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", "", 120, -0.1, 1.1, 120, -0.1, 1.1, 120, -0.1, 1.1);
 143   TH3* h3Ddt=new TH3F("h3Ddt", "", 360, -120., 120., 360, -120., 120., 360, -120., 120.);
 144
 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);
 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     // 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);
 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;
 176     //Float_t twopi=2*TMath::Pi();
 177     //Float_t sqrtwopi=TMath::Sqrt(twopi);
 178     //Double_t p0=1./(p2*sqrtwopi);
 179     // choising a different normalization
 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(px*px+py*py+pz*pz);
 190     Float_t ptvtx=TMath::Sqrt(px*px+py*py);
 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 //__________________________________________________________________
 265 Bool_t AliTOFProb::operator==( AliTOFProb const & tofrec)const
 266 {
 267   // dummy version of Equal operator.
 268   // requested by coding conventions
 269   return kTRUE;
 270
 271 }