[u/mrichter/AliRoot.git] / JETAN / AliJetKineReader.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.                  *
 **************************************************************************/

//-------------------------------------------------------------------------
// Jet kinematics reader 
// MC reader for jet analysis
// Author: Andreas Morsch (andreas.morsch@cern.ch)
//-------------------------------------------------------------------------

// From root ...
#include <TClonesArray.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TParticlePDG.h>
#include <TLorentzVector.h>
#include <TSystem.h>
#include <TDatabasePDG.h>
#include <TRandom.h>
// From AliRoot ...
#include "AliJet.h"
#include "AliJetKineReaderHeader.h"
#include "AliJetKineReader.h"
#include "AliRunLoader.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliGenPythiaEventHeader.h"

ClassImp(AliJetKineReader)

AliJetKineReader::AliJetKineReader():
    AliJetReader(),  
    fRunLoader(0x0),
    fAliHeader(0x0)
{
  // Default constructor
}

//____________________________________________________________________________

AliJetKineReader::~AliJetKineReader()
{
  // Destructor
  delete fAliHeader;
}

//____________________________________________________________________________

void AliJetKineReader::OpenInputFiles()
{
  // Opens the input file using the run loader
  const char* dirName = fReaderHeader->GetDirectory();
  char path[256];
  sprintf(path, "%s/galice.root",dirName);
  fRunLoader = AliRunLoader::Open(path);
  fRunLoader->LoadKinematics();
  fRunLoader->LoadHeader(); 
  
  // get number of events
  Int_t nMax = fRunLoader->GetNumberOfEvents();
  printf("\nTotal number of events = %d", nMax);
  
  // set number of events in header
  if (fReaderHeader->GetLastEvent() == -1)
    fReaderHeader->SetLastEvent(nMax);
  else {
    Int_t nUsr = fReaderHeader->GetLastEvent();
    fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
  }
}

//____________________________________________________________________________

Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
{
  // Fill momentum array for event
  char path[256];
  const char* dirName   = fReaderHeader->GetDirectory();
  const char* bgdirName = fReaderHeader->GetBgDirectory();
  Int_t goodTrack = 0;
  // Clear array
  // temporary storage of signal and cut flags
  Int_t* sflag  = new Int_t[50000];
  Int_t* cflag  = new Int_t[50000];
  Int_t  spb = 0;
  
  ClearArray();
  for (Int_t i = 0; i < 2; i++) {
      if (i == 1) {
// Pythia Events
	  if (fRunLoader) delete fRunLoader;
	  sprintf(path, "%s/galice.root", dirName);
	  fRunLoader = AliRunLoader::Open(path);
	  fRunLoader->LoadKinematics();
	  fRunLoader->LoadHeader(); 
	  // Get event from runloader
	  fRunLoader->GetEvent(event);
      } else {
// Background events
	  if ((spb = fReaderHeader->GetSignalPerBg()) > 0) {
	      if (fRunLoader) delete fRunLoader;
	      sprintf(path, "%s/galice.root", bgdirName);
	      fRunLoader = AliRunLoader::Open(path);
	      fRunLoader->LoadKinematics();
	      fRunLoader->LoadHeader(); 
	      // Get event from runloader
	      fRunLoader->GetEvent(event / spb);
	  } else {
	      continue;
	  }
      }

      // Get the stack
      AliStack* stack = fRunLoader->Stack();
      // Number of primaries
      Int_t nt = stack->GetNprimary();
      
      // Get cuts set by user and header
      Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
      Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
      Float_t etaMax = fReaderHeader->GetFiducialEtaMax();  
      fAliHeader = fRunLoader->GetHeader();
      
      
      TLorentzVector p4;
      // Loop over particles
      for (Int_t it = 0; it < nt; it++) {
	  TParticle *part = stack->Particle(it);
	  Int_t   status  = part->GetStatusCode();
	  Int_t   pdg     = TMath::Abs(part->GetPdgCode());
	  Float_t pt      = part->Pt(); 
	  
	  // Skip non-final state particles, neutrinos 
	  // and particles with pt < pt_min 
	  if ((status != 1)            
	      || (pdg == 12 || pdg == 14 || pdg == 16)) continue; 
	  
	  Float_t p       = part->P();
	  Float_t p0      = p;
	  Float_t eta     = part->Eta();
	  Float_t phi     = part->Phi();
	  
	  // Fast simulation of TPC if requested
	  if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
	      // Charged particles only
	      Float_t charge = 
		  TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
	      if (charge == 0)               continue;
	      // Simulate efficiency
	      if (!Efficiency(p0, eta, phi)) continue;
	      // Simulate resolution
	      p = SmearMomentum(4, p0);
	  } // End fast TPC
	  
	  // Fast simulation of EMCAL if requested
	  if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
	      Float_t charge = 
		  TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
	      // Charged particles only
	      if (charge != 0){
		  // Simulate efficiency
		  if (!Efficiency(p0, eta, phi)) continue;
		  // Simulate resolution
		  p = SmearMomentum(4, p0);
	      } // end "if" charged particles
	      // Neutral particles (exclude K0L, n, nbar)
	      if (pdg == kNeutron || pdg == kK0Long) continue;
	  } // End fast EMCAL
    	  
	  // Fill momentum array
	  Float_t r  = p/p0;
	  Float_t px = r * part->Px(); 
	  Float_t py = r * part->Py(); 
	  Float_t pz = r * part->Pz();
	  Float_t m  = part->GetMass();
	  Float_t e  = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
	  p4 = TLorentzVector(px, py, pz, e);
	  if ( (p4.Eta()>etaMax) || (p4.Eta()<etaMin)) continue;
	  new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
	  
	  // all tracks are signal ... ?
	  sflag[goodTrack]=1;
	  cflag[goodTrack]=0;
	  if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
	  goodTrack++;
      }
  }

  // set the signal flags
  fSignalFlag.Set(goodTrack, sflag);
  fCutFlag.Set(goodTrack, cflag);
  printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
  
  delete[] sflag;
  delete[] cflag;
  
  return kTRUE;
}


Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
{
  //  The relative momentum error, i.e. 
  //  (delta p)/p = sqrt (a**2 + (b*p)**2) * 10**-2,
  //  where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
  //  (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
  //
  //  If we include information from TRD, b will be a factor 2/3 smaller.
  //
  //  ind = 1: high multiplicity
  //  ind = 2: low  multiplicity
  //  ind = 3: high multiplicity + TRD
  //  ind = 4: low  multiplicity + TRD
  
  Float_t pSmeared;
  Float_t a = 0.75;
  Float_t b = 0.12;
  
  if (ind == 1) b = 0.12;
  if (ind == 2) b = 0.08;
  if (ind == 3) b = 0.12;    
  if (ind == 4) b = 0.08;    
  
  Float_t sigma =  p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
  pSmeared = p + gRandom->Gaus(0., sigma);
  return pSmeared;
}


Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /*eta*/, Float_t phi)
{
  // Fast simulation of geometrical acceptance and tracking efficiency
 
  //  Tracking
  Float_t eff = 0.99;
  if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
  // Geometry
  if (p > 0.5) {
    phi *= 180. / TMath::Pi();
    // Sector number 0 - 17
    Int_t isec  = Int_t(phi / 20.);
    // Sector centre
    Float_t phi0 = isec * 20. + 10.;
    Float_t phir = TMath::Abs(phi-phi0);
    // 2 deg of dead space
    if (phir > 9.) eff = 0.;
  }

  if (gRandom->Rndm() > eff) {
    return kFALSE;
  } else {
    return kTRUE;
  }

}

Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
{
    // Get generated jets from mc header
    AliHeader* alih = GetAliHeader(); 
    if (alih == 0) return kFALSE;
    AliGenEventHeader * genh = alih->GenEventHeader();
    if (genh == 0) return kFALSE;
    Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
    Int_t* m = new Int_t[nj];
    Int_t* k = new Int_t[nj];
    for (Int_t i=0; i< nj; i++) {
	Float_t p[4];
	((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
	genJets->AddJet(p[0],p[1],p[2],p[3]);
	m[i]=1;
	k[i]=i;
    }
    genJets->SetNinput(nj);
    genJets->SetMultiplicities(m);
    genJets->SetInJet(k);
    delete[] k;
    delete[] m;
    return kTRUE;
}
Commit	Line	Data
99e5fe42	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
83a444b1	16	//-------------------------------------------------------------------------
99e5fe42	17	// Jet kinematics reader
99e5fe42	18	// MC reader for jet analysis
83a444b1	19	// Author: Andreas Morsch (andreas.morsch@cern.ch)
83a444b1	20	//-------------------------------------------------------------------------
99e5fe42	21
	22	// From root ...
	23	#include <TClonesArray.h>
	24	#include <TPDGCode.h>
	25	#include <TParticle.h>
	26	#include <TParticlePDG.h>
99e5fe42	27	#include <TLorentzVector.h>
	28	#include <TSystem.h>
	29	#include <TDatabasePDG.h>
	30	#include <TRandom.h>
	31	// From AliRoot ...
0ffa8579	32	#include "AliJet.h"
99e5fe42	33	#include "AliJetKineReaderHeader.h"
	34	#include "AliJetKineReader.h"
	35	#include "AliRunLoader.h"
	36	#include "AliStack.h"
0ffa8579	37	#include "AliHeader.h"
0ffa8579	38	#include "AliGenPythiaEventHeader.h"
99e5fe42	39
	40	ClassImp(AliJetKineReader)
	41
1b7d5d7e	42	AliJetKineReader::AliJetKineReader():
b45b0c92	43	AliJetReader(),
	44	fRunLoader(0x0),
	45	fAliHeader(0x0)
99e5fe42	46	{
1b7d5d7e	47	// Default constructor
99e5fe42	48	}
	49
	50	//____________________________________________________________________________
	51
	52	AliJetKineReader::~AliJetKineReader()
	53	{
	54	// Destructor
0ffa8579	55	delete fAliHeader;
99e5fe42	56	}
	57
	58	//____________________________________________________________________________
	59
	60	void AliJetKineReader::OpenInputFiles()
	61	{
83a444b1	62	// Opens the input file using the run loader
	63	const char* dirName = fReaderHeader->GetDirectory();
	64	char path[256];
	65	sprintf(path, "%s/galice.root",dirName);
	66	fRunLoader = AliRunLoader::Open(path);
	67	fRunLoader->LoadKinematics();
	68	fRunLoader->LoadHeader();
	69
	70	// get number of events
	71	Int_t nMax = fRunLoader->GetNumberOfEvents();
	72	printf("\nTotal number of events = %d", nMax);
	73
99e5fe42	74	// set number of events in header
83a444b1	75	if (fReaderHeader->GetLastEvent() == -1)
	76	fReaderHeader->SetLastEvent(nMax);
	77	else {
	78	Int_t nUsr = fReaderHeader->GetLastEvent();
	79	fReaderHeader->SetLastEvent(TMath::Min(nMax, nUsr));
	80	}
99e5fe42	81	}
	82
	83	//____________________________________________________________________________
	84
0ffa8579	85	Bool_t AliJetKineReader::FillMomentumArray(Int_t event)
99e5fe42	86	{
83a444b1	87	// Fill momentum array for event
4faf4c3d	88	char path[256];
	89	const char* dirName = fReaderHeader->GetDirectory();
	90	const char* bgdirName = fReaderHeader->GetBgDirectory();
83a444b1	91	Int_t goodTrack = 0;
83a444b1	92	// Clear array
83a444b1	93	// temporary storage of signal and cut flags
4faf4c3d	94	Int_t* sflag = new Int_t[50000];
	95	Int_t* cflag = new Int_t[50000];
	96	Int_t spb = 0;
83a444b1	97
4faf4c3d	98	ClearArray();
	99	for (Int_t i = 0; i < 2; i++) {
	100	if (i == 1) {
	101	// Pythia Events
	102	if (fRunLoader) delete fRunLoader;
	103	sprintf(path, "%s/galice.root", dirName);
	104	fRunLoader = AliRunLoader::Open(path);
	105	fRunLoader->LoadKinematics();
	106	fRunLoader->LoadHeader();
	107	// Get event from runloader
	108	fRunLoader->GetEvent(event);
	109	} else {
	110	// Background events
	111	if ((spb = fReaderHeader->GetSignalPerBg()) > 0) {
	112	if (fRunLoader) delete fRunLoader;
	113	sprintf(path, "%s/galice.root", bgdirName);
	114	fRunLoader = AliRunLoader::Open(path);
	115	fRunLoader->LoadKinematics();
	116	fRunLoader->LoadHeader();
	117	// Get event from runloader
	118	fRunLoader->GetEvent(event / spb);
	119	} else {
	120	continue;
	121	}
	122	}
	123
	124	// Get the stack
	125	AliStack* stack = fRunLoader->Stack();
	126	// Number of primaries
	127	Int_t nt = stack->GetNprimary();
	128
	129	// Get cuts set by user and header
	130	Double_t ptMin = ((AliJetKineReaderHeader*) fReaderHeader)->GetPtCut();
	131	Float_t etaMin = fReaderHeader->GetFiducialEtaMin();
	132	Float_t etaMax = fReaderHeader->GetFiducialEtaMax();
	133	fAliHeader = fRunLoader->GetHeader();
	134
	135
	136	TLorentzVector p4;
	137	// Loop over particles
	138	for (Int_t it = 0; it < nt; it++) {
	139	TParticle *part = stack->Particle(it);
	140	Int_t status = part->GetStatusCode();
	141	Int_t pdg = TMath::Abs(part->GetPdgCode());
	142	Float_t pt = part->Pt();
	143
	144	// Skip non-final state particles, neutrinos
	145	// and particles with pt < pt_min
	146	if ((status != 1)
	147	\|\| (pdg == 12 \|\| pdg == 14 \|\| pdg == 16)) continue;
	148
	149	Float_t p = part->P();
	150	Float_t p0 = p;
	151	Float_t eta = part->Eta();
	152	Float_t phi = part->Phi();
	153
	154	// Fast simulation of TPC if requested
	155	if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimTPC()) {
	156	// Charged particles only
	157	Float_t charge =
	158	TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
	159	if (charge == 0) continue;
	160	// Simulate efficiency
	161	if (!Efficiency(p0, eta, phi)) continue;
162	// Simulate resolution
163	p = SmearMomentum(4, p0);
164	} // End fast TPC
165
166	// Fast simulation of EMCAL if requested
167	if (((AliJetKineReaderHeader*)fReaderHeader)->FastSimEMCAL()) {
168	Float_t charge =
169	TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
170	// Charged particles only
171	if (charge != 0){
172	// Simulate efficiency
173	if (!Efficiency(p0, eta, phi)) continue;
174	// Simulate resolution
175	p = SmearMomentum(4, p0);
176	} // end "if" charged particles
177	// Neutral particles (exclude K0L, n, nbar)
178	if (pdg == kNeutron \|\| pdg == kK0Long) continue;
179	} // End fast EMCAL
180
181	// Fill momentum array
182	Float_t r = p/p0;
183	Float_t px = r * part->Px();
184	Float_t py = r * part->Py();
185	Float_t pz = r * part->Pz();
186	Float_t m = part->GetMass();
187	Float_t e = TMath::Sqrt(px * px + py * py + pz * pz + m * m);
188	p4 = TLorentzVector(px, py, pz, e);
189	if ( (p4.Eta()>etaMax) \|\| (p4.Eta()<etaMin)) continue;
190	new ((*fMomentumArray)[goodTrack]) TLorentzVector(p4);
191
192	// all tracks are signal ... ?
193	sflag[goodTrack]=1;
194	cflag[goodTrack]=0;
195	if (pt > ptMin) cflag[goodTrack]=1; // track surviving pt cut
196	goodTrack++;
197	}
99e5fe42	198	}
83a444b1	199
83a444b1	200	// set the signal flags
63a9dacf	201	fSignalFlag.Set(goodTrack, sflag);
63a9dacf	202	fCutFlag.Set(goodTrack, cflag);
83a444b1	203	printf("\nIn event %d, number of good tracks %d \n", event, goodTrack);
63a9dacf	204
	205	delete[] sflag;
	206	delete[] cflag;
	207
0ffa8579	208	return kTRUE;
99e5fe42	209	}
	210
	211
48bb52d8	212	Float_t AliJetKineReader::SmearMomentum(Int_t ind, Float_t p)
48bb52d8	213	{
83a444b1	214	// The relative momentum error, i.e.
	215	// (delta p)/p = sqrt (a*2 + (bp)*2) 10**-2,
	216	// where typically a = 0.75 and b = 0.16 - 0.24 depending on multiplicity
	217	// (the lower value is for dn/d(eta) about 2000, and the higher one for 8000)
	218	//
	219	// If we include information from TRD, b will be a factor 2/3 smaller.
	220	//
	221	// ind = 1: high multiplicity
	222	// ind = 2: low multiplicity
	223	// ind = 3: high multiplicity + TRD
	224	// ind = 4: low multiplicity + TRD
	225
	226	Float_t pSmeared;
	227	Float_t a = 0.75;
	228	Float_t b = 0.12;
	229
	230	if (ind == 1) b = 0.12;
	231	if (ind == 2) b = 0.08;
	232	if (ind == 3) b = 0.12;
	233	if (ind == 4) b = 0.08;
	234
	235	Float_t sigma = p * TMath::Sqrt(a * a + b * b * p * p)*0.01;
	236	pSmeared = p + gRandom->Gaus(0., sigma);
	237	return pSmeared;
48bb52d8	238	}
83a444b1	239
83a444b1	240
48bb52d8	241	Bool_t AliJetKineReader::Efficiency(Float_t p, Float_t /eta/, Float_t phi)
48bb52d8	242	{
83a444b1	243	// Fast simulation of geometrical acceptance and tracking efficiency
	244
	245	// Tracking
	246	Float_t eff = 0.99;
	247	if (p < 0.5) eff -= (0.5-p)*0.2/0.3;
	248	// Geometry
	249	if (p > 0.5) {
	250	phi *= 180. / TMath::Pi();
	251	// Sector number 0 - 17
	252	Int_t isec = Int_t(phi / 20.);
	253	// Sector centre
	254	Float_t phi0 = isec * 20. + 10.;
	255	Float_t phir = TMath::Abs(phi-phi0);
	256	// 2 deg of dead space
	257	if (phir > 9.) eff = 0.;
	258	}
	259
	260	if (gRandom->Rndm() > eff) {
	261	return kFALSE;
	262	} else {
	263	return kTRUE;
	264	}
	265
48bb52d8	266	}
48bb52d8	267
0ffa8579	268	Bool_t AliJetKineReader::GetGenJets(AliJet* genJets)
	269	{
	270	// Get generated jets from mc header
	271	AliHeader* alih = GetAliHeader();
	272	if (alih == 0) return kFALSE;
	273	AliGenEventHeader * genh = alih->GenEventHeader();
	274	if (genh == 0) return kFALSE;
	275	Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
	276	Int_t* m = new Int_t[nj];
	277	Int_t* k = new Int_t[nj];
	278	for (Int_t i=0; i< nj; i++) {
	279	Float_t p[4];
	280	((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
	281	genJets->AddJet(p[0],p[1],p[2],p[3]);
	282	m[i]=1;
	283	k[i]=i;
	284	}
	285	genJets->SetNinput(nj);
	286	genJets->SetMultiplicities(m);
	287	genJets->SetInJet(k);
63a9dacf	288	delete[] k;
63a9dacf	289	delete[] m;
0ffa8579	290	return kTRUE;
0ffa8579	291	}