[u/mrichter/AliRoot.git] / HBTAN / AliHBTWeightsPID.cxx

/* $Id $ */

//-----------------------------------------------------------
//This class introduces the weights calculated according 
//with functions of efficiency of identification (TPC+TOF) 
//(calculated by B.V. Batyunia).
//Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru)
//-----------------------------------------------------------

#include "AliHBTWeightsPID.h"
#include "AliHBTPair.h"
#include <TRandom.h>
#include <TMath.h>
#include <TH1.h>
#include <TF1.h>

ClassImp(AliHBTWeightsPID)  
 
AliHBTWeightsPID* AliHBTWeightsPID::fgWeightsPID=NULL; 

AliHBTWeightsPID::AliHBTWeightsPID() :
  fEfficTPC1(0),
  fEfficTPC2(0),
  fEfficTOF1(0),
  fEfficTOF2(0),
  fPtK(0),
  fPtKefftpc(0),
  fPtKefftpcboth(0),
  fEffic1pol(0),
  fEffic2pol(0),
  fEffic3pol(0),
  fEffic4pol(0),
  fEffic1polTOF(0),
  fEffic2polTOF(0),
  fEffic3polTOF(0),
  fEffic4polTOF(0)
{
  //ctor
  //initial parameters of model

  fPtK   = new TH1F("fPtK"," pt of K+- ",40,0,4);
  fPtKefftpc   = new TH1F("fPtEfficKpos"," pt of K+ after efficiency cut ",40,0,4);
  fPtKefftpcboth   = new TH1F("fPtEfficKneg"," pt of K- after efficiency cut ",40,0,4);

  // efficiency functions (from pt) for K+- in the TPC
  fEffic1pol = new TF1("fEffic1pol4","pol4",100.,480.); 
  // for -0.9 < eta < -0.7
  // for  0.7 < eta <  0.9
  fEffic1pol->SetParameter(0,-0.266362);
  fEffic1pol->SetParameter(1,0.00565461);
  fEffic1pol->SetParameter(2,-4.06686e-05);
  fEffic1pol->SetParameter(3,1.39387e-07);
  fEffic1pol->SetParameter(4,-1.59674e-10);

  fEffic2pol = new TF1("fEffic2pol4","pol4",100.,540.); 

  fEffic2pol->SetParameter(0,-0.324881);
  fEffic2pol->SetParameter(1,0.00565381);
  fEffic2pol->SetParameter(2,-3.23633e-05);
  fEffic2pol->SetParameter(3,9.72523e-08);
  fEffic2pol->SetParameter(4,-1.01013e-10);
                                                                               
  fEffic3pol = new TF1("fEffic3pol4","pol4",100.,585.); 
  // for -0.5 < eta < -0.3
  // for  0.3 < eta <  0.5
  fEffic3pol->SetParameter(0,-0.306572);
  fEffic3pol->SetParameter(1,0.00557472);
  fEffic3pol->SetParameter(2,-3.33752e-05);
  fEffic3pol->SetParameter(3,9.83241e-08);
  fEffic3pol->SetParameter(4,-9.5827e-11);
  fEffic4pol = new TF1("fEffic4pol4","pol4",100.,600.);
  // for -0.3 < eta < 0.3
  fEffic4pol->SetParameter(0,-0.168648);
  fEffic4pol->SetParameter(1,0.00252021);
  fEffic4pol->SetParameter(2,-1.09113e-05);
  fEffic4pol->SetParameter(3,3.34871e-08);
  fEffic4pol->SetParameter(4,-3.31691e-11);

  // efficiency functions (from pt) for K+- in the TOF
  fEffic1polTOF = new TF1("fEffic1pol4TOF","pol4",0.2,2.0); 
  // for -0.9 < eta < -0.7
  // for  0.7 < eta <  0.9
  fEffic1polTOF->SetParameter(0,-0.165885);
  fEffic1polTOF->SetParameter(1,0.717459);
  fEffic1polTOF->SetParameter(2,-0.457131);
  fEffic1polTOF->SetParameter(3,0.284753);
  fEffic1polTOF->SetParameter(4,-0.105215);
                                                                               
  fEffic2polTOF = new TF1("fEffic2pol4TOF","pol4",0.2,2.4); 
  // for -0.7 < eta < -0.5
  // for  0.5 < eta <  0.7

  fEffic2polTOF->SetParameter(0,-0.165947);
  fEffic2polTOF->SetParameter(1,0.702475);
  fEffic2polTOF->SetParameter(2,-0.300313);
  fEffic2polTOF->SetParameter(3,0.127047);
  fEffic2polTOF->SetParameter(4,-0.0489395);
                                                                               

  fEffic3polTOF = new TF1("fEffic3pol4TOF","pol4",0.2,2.4); 
  // for -0.5 < eta < -0.3
  // for  0.3 < eta <  0.5
  fEffic3polTOF->SetParameter(0,-0.339516);
  fEffic3polTOF->SetParameter(1,1.56942);
  fEffic3polTOF->SetParameter(2,-1.43132);
  fEffic3polTOF->SetParameter(3,0.727148);
  fEffic3polTOF->SetParameter(4,-0.158444);
  fEffic4polTOF = new TF1("fEffic4pol4TOF","pol4",0.2,2.6);
  // for -0.3 < eta < 0.3
  fEffic4polTOF->SetParameter(0,-0.243435);
  fEffic4polTOF->SetParameter(1,1.00928);
  fEffic4polTOF->SetParameter(2,-0.594597);
  fEffic4polTOF->SetParameter(3,0.212601);
  fEffic4polTOF->SetParameter(4,-0.0453419);

}                                            

AliHBTWeightsPID* AliHBTWeightsPID::Instance()
{
  //Creates an instance of the class 
  //or returns pointer to already existing one
  if (fgWeightsPID) {
    return fgWeightsPID;
  } else {
    fgWeightsPID = new AliHBTWeightsPID();
    return fgWeightsPID;
  }
}


Double_t AliHBTWeightsPID::GetWeightPID(const AliHBTPair* trackpair)
{
  //Calculates the weight of "trackpair"
  AliVAODParticle *track1 = trackpair->Particle1();
  AliVAODParticle *track2 = trackpair->Particle2();
    
  Double_t pt1=track1->Pt();
  Double_t eta1=track1->Eta();

  fPtK->Fill(pt1); 

  if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && pt1 < 0.48) {
    fEfficTPC1 = fEffic1pol->Eval(pt1*1000.);
  }else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && pt1 < 0.54) {
    fEfficTPC1 = fEffic2pol->Eval(pt1*1000.);
  }else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && pt1 < 0.585) {
    fEfficTPC1 = fEffic3pol->Eval(pt1*1000.);
  }else if(eta1 > -0.3 && eta1 < 0.3 && pt1 < 0.6) {
    fEfficTPC1 = fEffic4pol->Eval(pt1*1000.);
  }
  
  // TOF efficiency
  
  if(pt1 > 0.2) {
    if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && pt1 < 2.0) {
      fEfficTOF1 = fEffic1polTOF->Eval(pt1);
    }else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && pt1 < 2.4) {
      fEfficTOF1 = fEffic2polTOF->Eval(pt1);
    }else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && pt1 < 2.4) {
      fEfficTOF1 = fEffic3polTOF->Eval(pt1);
    }else if(eta1 > -0.3 && eta1 < 0.3 && pt1 < 2.6) {
      fEfficTOF1 = fEffic4polTOF->Eval(pt1);
    }
  }
  
  Double_t rndmtpc=gRandom->Rndm();
  Double_t rndmtof=gRandom->Rndm();
  Double_t weightPID1=1.;
  if(fEfficTPC1 < rndmtpc && fEfficTOF1 < rndmtof) { weightPID1=0.;}
  fPtKefftpc->Fill(pt1,weightPID1); 
  
  
  Double_t pt2=track2->Pt();
  Double_t eta2=track2->Eta();
  
  // TPC efficiency

  if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && pt2 < 0.48) {
    fEfficTPC1 = fEffic1pol->Eval(pt2*1000.);
  }else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && pt2 < 0.54) {
    fEfficTPC1 = fEffic2pol->Eval(pt2*1000.);
  }else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && pt2 < 0.585) {
    fEfficTPC1 = fEffic3pol->Eval(pt2*1000.);
  }else if(eta2 > -0.3 && eta2 < 0.3 && pt2 < 0.6) {
    fEfficTPC1 = fEffic4pol->Eval(pt2*1000.);
  }
  
  // TOF efficiency
  
  if(pt2 > 0.2) {
    if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && pt2 < 2.0) {
      fEfficTOF1 = fEffic1polTOF->Eval(pt2);
    }else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && pt2 < 2.4) {
      fEfficTOF1 = fEffic2polTOF->Eval(pt2);
    }else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && pt2 < 2.4) {
      fEfficTOF1 = fEffic3polTOF->Eval(pt2);
    }else if(eta2 > -0.3 && eta2 < 0.3 && pt2 < 2.6) {
      fEfficTOF1 = fEffic4polTOF->Eval(pt2);
    }
  }
  
  rndmtpc=gRandom->Rndm();
  rndmtof=gRandom->Rndm();
  
  Double_t weightPID2=1.;
  if(fEfficTPC1 < rndmtpc && fEfficTOF1 < rndmtof) { weightPID2=0.;}
  
  Double_t weightPID=weightPID1*weightPID2;
  fPtKefftpcboth->Fill(pt1,weightPID); 
  
  return weightPID;
}
Commit	Line	Data
99927f15	1	/* $Id $ */
	2
	3	//-----------------------------------------------------------
	4	//This class introduces the weights calculated according
	5	//with functions of efficiency of identification (TPC+TOF)
	6	//(calculated by B.V. Batyunia).
a6e49985	7	//Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru)
99927f15	8	//-----------------------------------------------------------
99927f15	9
dd82cadc	10	#include "AliHBTWeightsPID.h"
a6e49985	11	#include "AliHBTPair.h"
99927f15	12	#include <TRandom.h>
	13	#include <TMath.h>
	14	#include <TH1.h>
	15	#include <TF1.h>
a6e49985	16
dd82cadc	17	ClassImp(AliHBTWeightsPID)
a6e49985	18
dd82cadc	19	AliHBTWeightsPID* AliHBTWeightsPID::fgWeightsPID=NULL;
a6e49985	20
4b1c9620	21	AliHBTWeightsPID::AliHBTWeightsPID() :
	22	fEfficTPC1(0),
	23	fEfficTPC2(0),
	24	fEfficTOF1(0),
	25	fEfficTOF2(0),
	26	fPtK(0),
	27	fPtKefftpc(0),
	28	fPtKefftpcboth(0),
	29	fEffic1pol(0),
	30	fEffic2pol(0),
	31	fEffic3pol(0),
	32	fEffic4pol(0),
	33	fEffic1polTOF(0),
	34	fEffic2polTOF(0),
	35	fEffic3polTOF(0),
	36	fEffic4polTOF(0)
99927f15	37	{
	38	//ctor
	39	//initial parameters of model
	40
	41	fPtK = new TH1F("fPtK"," pt of K+- ",40,0,4);
	42	fPtKefftpc = new TH1F("fPtEfficKpos"," pt of K+ after efficiency cut ",40,0,4);
	43	fPtKefftpcboth = new TH1F("fPtEfficKneg"," pt of K- after efficiency cut ",40,0,4);
	44
	45	// efficiency functions (from pt) for K+- in the TPC
	46	fEffic1pol = new TF1("fEffic1pol4","pol4",100.,480.);
	47	// for -0.9 < eta < -0.7
	48	// for 0.7 < eta < 0.9
	49	fEffic1pol->SetParameter(0,-0.266362);
	50	fEffic1pol->SetParameter(1,0.00565461);
	51	fEffic1pol->SetParameter(2,-4.06686e-05);
	52	fEffic1pol->SetParameter(3,1.39387e-07);
	53	fEffic1pol->SetParameter(4,-1.59674e-10);
	54
	55	fEffic2pol = new TF1("fEffic2pol4","pol4",100.,540.);
	56
	57	fEffic2pol->SetParameter(0,-0.324881);
	58	fEffic2pol->SetParameter(1,0.00565381);
	59	fEffic2pol->SetParameter(2,-3.23633e-05);
	60	fEffic2pol->SetParameter(3,9.72523e-08);
	61	fEffic2pol->SetParameter(4,-1.01013e-10);
a6e49985	62
99927f15	63	fEffic3pol = new TF1("fEffic3pol4","pol4",100.,585.);
	64	// for -0.5 < eta < -0.3
	65	// for 0.3 < eta < 0.5
	66	fEffic3pol->SetParameter(0,-0.306572);
	67	fEffic3pol->SetParameter(1,0.00557472);
	68	fEffic3pol->SetParameter(2,-3.33752e-05);
	69	fEffic3pol->SetParameter(3,9.83241e-08);
	70	fEffic3pol->SetParameter(4,-9.5827e-11);
	71	fEffic4pol = new TF1("fEffic4pol4","pol4",100.,600.);
	72	// for -0.3 < eta < 0.3
	73	fEffic4pol->SetParameter(0,-0.168648);
	74	fEffic4pol->SetParameter(1,0.00252021);
	75	fEffic4pol->SetParameter(2,-1.09113e-05);
	76	fEffic4pol->SetParameter(3,3.34871e-08);
	77	fEffic4pol->SetParameter(4,-3.31691e-11);
a6e49985	78
a6e49985	79	// efficiency functions (from pt) for K+- in the TOF
99927f15	80	fEffic1polTOF = new TF1("fEffic1pol4TOF","pol4",0.2,2.0);
	81	// for -0.9 < eta < -0.7
	82	// for 0.7 < eta < 0.9
	83	fEffic1polTOF->SetParameter(0,-0.165885);
	84	fEffic1polTOF->SetParameter(1,0.717459);
	85	fEffic1polTOF->SetParameter(2,-0.457131);
	86	fEffic1polTOF->SetParameter(3,0.284753);
	87	fEffic1polTOF->SetParameter(4,-0.105215);
a6e49985	88
99927f15	89	fEffic2polTOF = new TF1("fEffic2pol4TOF","pol4",0.2,2.4);
	90	// for -0.7 < eta < -0.5
	91	// for 0.5 < eta < 0.7
	92
	93	fEffic2polTOF->SetParameter(0,-0.165947);
	94	fEffic2polTOF->SetParameter(1,0.702475);
	95	fEffic2polTOF->SetParameter(2,-0.300313);
	96	fEffic2polTOF->SetParameter(3,0.127047);
	97	fEffic2polTOF->SetParameter(4,-0.0489395);
a6e49985	98
a6e49985	99
99927f15	100	fEffic3polTOF = new TF1("fEffic3pol4TOF","pol4",0.2,2.4);
	101	// for -0.5 < eta < -0.3
	102	// for 0.3 < eta < 0.5
	103	fEffic3polTOF->SetParameter(0,-0.339516);
	104	fEffic3polTOF->SetParameter(1,1.56942);
	105	fEffic3polTOF->SetParameter(2,-1.43132);
	106	fEffic3polTOF->SetParameter(3,0.727148);
	107	fEffic3polTOF->SetParameter(4,-0.158444);
	108	fEffic4polTOF = new TF1("fEffic4pol4TOF","pol4",0.2,2.6);
	109	// for -0.3 < eta < 0.3
	110	fEffic4polTOF->SetParameter(0,-0.243435);
	111	fEffic4polTOF->SetParameter(1,1.00928);
	112	fEffic4polTOF->SetParameter(2,-0.594597);
	113	fEffic4polTOF->SetParameter(3,0.212601);
	114	fEffic4polTOF->SetParameter(4,-0.0453419);
	115
	116	}
a6e49985	117
dd82cadc	118	AliHBTWeightsPID* AliHBTWeightsPID::Instance()
99927f15	119	{
	120	//Creates an instance of the class
	121	//or returns pointer to already existing one
	122	if (fgWeightsPID) {
	123	return fgWeightsPID;
	124	} else {
dd82cadc	125	fgWeightsPID = new AliHBTWeightsPID();
99927f15	126	return fgWeightsPID;
	127	}
	128	}
a6e49985	129
a6e49985	130
dd82cadc	131	Double_t AliHBTWeightsPID::GetWeightPID(const AliHBTPair* trackpair)
a6e49985	132	{
99927f15	133	//Calculates the weight of "trackpair"
78d7c6d3	134	AliVAODParticle *track1 = trackpair->Particle1();
78d7c6d3	135	AliVAODParticle *track2 = trackpair->Particle2();
a6e49985	136
99927f15	137	Double_t pt1=track1->Pt();
	138	Double_t eta1=track1->Eta();
	139
	140	fPtK->Fill(pt1);
	141
	142	if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && pt1 < 0.48) {
	143	fEfficTPC1 = fEffic1pol->Eval(pt1*1000.);
	144	}else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && pt1 < 0.54) {
	145	fEfficTPC1 = fEffic2pol->Eval(pt1*1000.);
	146	}else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && pt1 < 0.585) {
	147	fEfficTPC1 = fEffic3pol->Eval(pt1*1000.);
	148	}else if(eta1 > -0.3 && eta1 < 0.3 && pt1 < 0.6) {
	149	fEfficTPC1 = fEffic4pol->Eval(pt1*1000.);
	150	}
	151
	152	// TOF efficiency
	153
	154	if(pt1 > 0.2) {
	155	if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && pt1 < 2.0) {
	156	fEfficTOF1 = fEffic1polTOF->Eval(pt1);
	157	}else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && pt1 < 2.4) {
	158	fEfficTOF1 = fEffic2polTOF->Eval(pt1);
	159	}else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && pt1 < 2.4) {
	160	fEfficTOF1 = fEffic3polTOF->Eval(pt1);
	161	}else if(eta1 > -0.3 && eta1 < 0.3 && pt1 < 2.6) {
	162	fEfficTOF1 = fEffic4polTOF->Eval(pt1);
	163	}
	164	}
	165
	166	Double_t rndmtpc=gRandom->Rndm();
	167	Double_t rndmtof=gRandom->Rndm();
	168	Double_t weightPID1=1.;
	169	if(fEfficTPC1 < rndmtpc && fEfficTOF1 < rndmtof) { weightPID1=0.;}
	170	fPtKefftpc->Fill(pt1,weightPID1);
	171
	172
	173	Double_t pt2=track2->Pt();
	174	Double_t eta2=track2->Eta();
	175
	176	// TPC efficiency
	177
	178	if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && pt2 < 0.48) {
	179	fEfficTPC1 = fEffic1pol->Eval(pt2*1000.);
	180	}else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && pt2 < 0.54) {
	181	fEfficTPC1 = fEffic2pol->Eval(pt2*1000.);
	182	}else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && pt2 < 0.585) {
	183	fEfficTPC1 = fEffic3pol->Eval(pt2*1000.);
	184	}else if(eta2 > -0.3 && eta2 < 0.3 && pt2 < 0.6) {
	185	fEfficTPC1 = fEffic4pol->Eval(pt2*1000.);
	186	}
	187
	188	// TOF efficiency
	189
	190	if(pt2 > 0.2) {
	191	if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && pt2 < 2.0) {
	192	fEfficTOF1 = fEffic1polTOF->Eval(pt2);
	193	}else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && pt2 < 2.4) {
	194	fEfficTOF1 = fEffic2polTOF->Eval(pt2);
	195	}else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && pt2 < 2.4) {
	196	fEfficTOF1 = fEffic3polTOF->Eval(pt2);
	197	}else if(eta2 > -0.3 && eta2 < 0.3 && pt2 < 2.6) {
	198	fEfficTOF1 = fEffic4polTOF->Eval(pt2);
	199	}
	200	}
201
202	rndmtpc=gRandom->Rndm();
203	rndmtof=gRandom->Rndm();
204
205	Double_t weightPID2=1.;
206	if(fEfficTPC1 < rndmtpc && fEfficTOF1 < rndmtof) { weightPID2=0.;}
207
208	Double_t weightPID=weightPID1*weightPID2;
209	fPtKefftpcboth->Fill(pt1,weightPID);
210
211	return weightPID;
a6e49985	212	}