[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 "AliAODJet.h"
#include "AliPDG.h"
#include "AliJetKineReaderHeader.h"
#include "AliJetKineReader.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliHeader.h"
#include "AliGenPythiaEventHeader.h"
#include "AliLog.h"

ClassImp(AliJetKineReader)

AliJetKineReader::AliJetKineReader():
    AliJetReader(),  
    fAliHeader(0),
    fMCEvent(0),
    fGenJets(0)
{
  // Default constructor
}

//____________________________________________________________________________

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

//____________________________________________________________________________

void AliJetKineReader::OpenInputFiles()
{
  // Opens the input file using the run loader
  // Not used anymore   
}

//____________________________________________________________________________

Bool_t AliJetKineReader::FillMomentumArray()
{
  // Fill momentum array for event
    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];
  
    ClearArray();
    // Get the stack
    AliStack* stack = fMCEvent->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 = fMCEvent->Header();
      
    
    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();


	if (((AliJetKineReaderHeader*)fReaderHeader)->ChargedOnly()) {
	    // Charged particles only
	    Float_t charge = 
		TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
	    if (charge == 0)               continue;
	} // End charged only

	
	// 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();
	TParticlePDG *pPDG = part->GetPDG();
	Float_t m = 0;
	if(!pPDG){
	  // this is very rare...
	  // Is avoided by AliPDG::AddParticlesToPdgDataBase();
	  // but this should be called only once... (how in proof?)
	  // Calucluate mass with unsmeared momentum values
	  m  = part->Energy()*part->Energy() - 
	    (px * px + py * py + pz * pz)/r/r; 
	  if(m>0)m = TMath::Sqrt(m);
	  else m = 0;
	  AliInfo(Form("Unknown Particle using %d calculated mass m = %3.3f",part->GetPdgCode(),m));

	}
	else{
	  m  = pPDG->Mass();
	}
	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++;
    } // track loop

// set the signal flags
    fSignalFlag.Set(goodTrack, sflag);
    fCutFlag.Set(goodTrack, cflag);
    AliInfo(Form(" Number of good tracks %d \n", 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;
  }

}

TClonesArray*  AliJetKineReader::GetGeneratedJets()
{
    //
    // Get generated jets from mc header
    //
    if (!fGenJets) {
	fGenJets = new TClonesArray("AliAODJets", 100);
    } else {
	fGenJets->Clear();
    }
    
    
    AliHeader* alih = GetAliHeader(); 
    if (alih == 0) return 0;
    AliGenEventHeader * genh = alih->GenEventHeader();
    if (genh == 0) return 0;

    Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets(); 
    for (Int_t i = 0; i < nj; i++) {
	Float_t p[4];
	((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
	new ((*fGenJets)[i]) AliAODJet(p[0], p[1], p[2], p[3]);
    }

    return fGenJets;
}

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