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

//Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru)
#include "AliHBTLLWeightsPID.h"
#include "AliPDG.h"
#include "AliHBTPair.h"
#include "AliHBTParticle.h"
#include <TList.h>
#include <TRandom.h>                                                                     
#include <TMath.h>                                                                       


ClassImp(AliHBTLLWeightsPID)  
 
AliHBTLLWeightsPID* AliHBTLLWeightsPID::fWeightsPID=NULL; 

AliHBTLLWeightsPID::AliHBTLLWeightsPID()
{                                                                                           
//ctor
//initial parameters of model

  ptK   = new TH1F("ptK"," pt of K+- ",40,0,4);
  ptKefftpc   = new TH1F("ptEfficKpos"," pt of K+ after efficiency cut ",40,0,4);
  ptKefftpcboth   = new TH1F("ptEfficKneg"," pt of K- after efficiency cut ",40,0,4);

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

  effic2pol = new TF1("effic2pol4","pol4",100.,540.); 

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

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

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

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

}                                                                                             
                                            

AliHBTLLWeightsPID* AliHBTLLWeightsPID::Instance()
{                                                                                             
  if (fWeightsPID) {                                                                        
    return fWeightsPID;                                                                   
   } else {                                                                                  
   fWeightsPID = new AliHBTLLWeightsPID();                                                        
      return fWeightsPID;                                                                   
  }                                                                                         
}                                                                                             


Double_t AliHBTLLWeightsPID::GetWeightPID(const AliHBTPair* trackpair)
{

    AliHBTParticle *track1 = trackpair->Particle1();
    AliHBTParticle *track2 = trackpair->Particle2();
    
    Double_t Pt1=track1->Pt();
    Double_t eta1=track1->Eta();

    ptK->Fill(Pt1); 

          if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && Pt1 < 0.48) {
            efficTPC1 = effic1pol->Eval(Pt1*1000.);
          }else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && Pt1 < 0.54) {
            efficTPC1 = effic2pol->Eval(Pt1*1000.);
          }else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && Pt1 < 0.585) {
            efficTPC1 = effic3pol->Eval(Pt1*1000.);
          }else if(eta1 > -0.3 && eta1 < 0.3 && Pt1 < 0.6) {
            efficTPC1 = effic4pol->Eval(Pt1*1000.);
          }

          // TOF efficiency

          if(Pt1 > 0.2) {
           if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && Pt1 < 2.0) {
            efficTOF1 = effic1polTOF->Eval(Pt1);
           }else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && Pt1 < 2.4) {
            efficTOF1 = effic2polTOF->Eval(Pt1);
           }else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && Pt1 < 2.4) {
            efficTOF1 = effic3polTOF->Eval(Pt1);
           }else if(eta1 > -0.3 && eta1 < 0.3 && Pt1 < 2.6) {
            efficTOF1 = effic4polTOF->Eval(Pt1);
           }
          }

          Double_t rndmtpc=gRandom->Rndm();
          Double_t rndmtof=gRandom->Rndm();
          Double_t weightPID1=1.;
          if(efficTPC1 < rndmtpc && efficTOF1 < rndmtof) { weightPID1=0.;}
          ptKefftpc->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) {
            efficTPC1 = effic1pol->Eval(Pt2*1000.);
          }else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && Pt2 < 0.54) {
            efficTPC1 = effic2pol->Eval(Pt2*1000.);
          }else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && Pt2 < 0.585) {
            efficTPC1 = effic3pol->Eval(Pt2*1000.);
          }else if(eta2 > -0.3 && eta2 < 0.3 && Pt2 < 0.6) {
            efficTPC1 = effic4pol->Eval(Pt2*1000.);
          }

          // TOF efficiency

          if(Pt2 > 0.2) {
           if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && Pt2 < 2.0) {
            efficTOF1 = effic1polTOF->Eval(Pt2);
           }else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && Pt2 < 2.4) {
            efficTOF1 = effic2polTOF->Eval(Pt2);
           }else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && Pt2 < 2.4) {
            efficTOF1 = effic3polTOF->Eval(Pt2);
           }else if(eta2 > -0.3 && eta2 < 0.3 && Pt2 < 2.6) {
            efficTOF1 = effic4polTOF->Eval(Pt2);
           }
          }

          rndmtpc=gRandom->Rndm();
          rndmtof=gRandom->Rndm();

          Double_t weightPID2=1.;
          if(efficTPC1 < rndmtpc && efficTOF1 < rndmtof) { weightPID2=0.;}

          Double_t weightPID=weightPID1*weightPID2;
          ptKefftpcboth->Fill(Pt1,weightPID); 

          return weightPID;
}
Commit	Line	Data
a6e49985	1	//Author: Ludmila Malinina, JINR (malinina@sunhe.jinr.ru)
	2	#include "AliHBTLLWeightsPID.h"
	3	#include "AliPDG.h"
	4	#include "AliHBTPair.h"
	5	#include "AliHBTParticle.h"
	6	#include <TList.h>
	7	#include <TRandom.h>
	8	#include <TMath.h>
	9
	10
	11	ClassImp(AliHBTLLWeightsPID)
	12
	13	AliHBTLLWeightsPID* AliHBTLLWeightsPID::fWeightsPID=NULL;
	14
	15	AliHBTLLWeightsPID::AliHBTLLWeightsPID()
	16	{
	17	//ctor
	18	//initial parameters of model
	19
	20	ptK = new TH1F("ptK"," pt of K+- ",40,0,4);
	21	ptKefftpc = new TH1F("ptEfficKpos"," pt of K+ after efficiency cut ",40,0,4);
	22	ptKefftpcboth = new TH1F("ptEfficKneg"," pt of K- after efficiency cut ",40,0,4);
	23
	24	// efficiency functions (from pt) for K+- in the TPC
	25	effic1pol = new TF1("effic1pol4","pol4",100.,480.);
	26	// for -0.9 < eta < -0.7
	27	// for 0.7 < eta < 0.9
	28	effic1pol->SetParameter(0,-0.266362);
	29	effic1pol->SetParameter(1,0.00565461);
	30	effic1pol->SetParameter(2,-4.06686e-05);
	31	effic1pol->SetParameter(3,1.39387e-07);
	32	effic1pol->SetParameter(4,-1.59674e-10);
	33
	34	effic2pol = new TF1("effic2pol4","pol4",100.,540.);
	35
	36	effic2pol->SetParameter(0,-0.324881);
	37	effic2pol->SetParameter(1,0.00565381);
	38	effic2pol->SetParameter(2,-3.23633e-05);
	39	effic2pol->SetParameter(3,9.72523e-08);
	40	effic2pol->SetParameter(4,-1.01013e-10);
	41
	42	effic3pol = new TF1("effic3pol4","pol4",100.,585.);
	43	// for -0.5 < eta < -0.3
	44	// for 0.3 < eta < 0.5
	45	effic3pol->SetParameter(0,-0.306572);
	46	effic3pol->SetParameter(1,0.00557472);
	47	effic3pol->SetParameter(2,-3.33752e-05);
	48	effic3pol->SetParameter(3,9.83241e-08);
	49	effic3pol->SetParameter(4,-9.5827e-11);
	50
	51	effic4pol = new TF1("effic4pol4","pol4",100.,600.); // for -0.3 < eta < 0.3
	52	effic4pol->SetParameter(0,-0.168648);
	53	effic4pol->SetParameter(1,0.00252021);
	54	effic4pol->SetParameter(2,-1.09113e-05);
	55	effic4pol->SetParameter(3,3.34871e-08);
	56	effic4pol->SetParameter(4,-3.31691e-11);
	57
	58	// efficiency functions (from pt) for K+- in the TOF
	59	effic1polTOF = new TF1("effic1pol4TOF","pol4",0.2,2.0);
	60	// for -0.9 < eta < -0.7
	61	// for 0.7 < eta < 0.9
	62	effic1polTOF->SetParameter(0,-0.165885);
	63	effic1polTOF->SetParameter(1,0.717459);
	64	effic1polTOF->SetParameter(2,-0.457131);
65	effic1polTOF->SetParameter(3,0.284753);
66	effic1polTOF->SetParameter(4,-0.105215);
67
68	effic2polTOF = new TF1("effic2pol4TOF","pol4",0.2,2.4);
69	// for -0.7 < eta < -0.5
70	// for 0.5 < eta < 0.7
71
72	effic2polTOF->SetParameter(0,-0.165947);
73	effic2polTOF->SetParameter(1,0.702475);
74	effic2polTOF->SetParameter(2,-0.300313);
75	effic2polTOF->SetParameter(3,0.127047);
76	effic2polTOF->SetParameter(4,-0.0489395);
77
78
79	effic3polTOF = new TF1("effic3pol4TOF","pol4",0.2,2.4);
80	// for -0.5 < eta < -0.3
81	// for 0.3 < eta < 0.5
82	effic3polTOF->SetParameter(0,-0.339516);
83	effic3polTOF->SetParameter(1,1.56942);
84	effic3polTOF->SetParameter(2,-1.43132);
85	effic3polTOF->SetParameter(3,0.727148);
86	effic3polTOF->SetParameter(4,-0.158444);
87
88	effic4polTOF = new TF1("effic4pol4TOF","pol4",0.2,2.6); // for -0.3 < eta < 0.3
89	effic4polTOF->SetParameter(0,-0.243435);
90	effic4polTOF->SetParameter(1,1.00928);
91	effic4polTOF->SetParameter(2,-0.594597);
92	effic4polTOF->SetParameter(3,0.212601);
93	effic4polTOF->SetParameter(4,-0.0453419);
94
95	}
96
97
98	AliHBTLLWeightsPID* AliHBTLLWeightsPID::Instance()
99	{
100	if (fWeightsPID) {
101	return fWeightsPID;
102	} else {
103	fWeightsPID = new AliHBTLLWeightsPID();
104	return fWeightsPID;
105	}
106	}
107
108
109
110	Double_t AliHBTLLWeightsPID::GetWeightPID(const AliHBTPair* trackpair)
111	{
112
113	AliHBTParticle *track1 = trackpair->Particle1();
114	AliHBTParticle *track2 = trackpair->Particle2();
115
116	Double_t Pt1=track1->Pt();
117	Double_t eta1=track1->Eta();
118
119	ptK->Fill(Pt1);
120
121	if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && Pt1 < 0.48) {
122	efficTPC1 = effic1pol->Eval(Pt1*1000.);
123	}else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && Pt1 < 0.54) {
124	efficTPC1 = effic2pol->Eval(Pt1*1000.);
125	}else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && Pt1 < 0.585) {
126	efficTPC1 = effic3pol->Eval(Pt1*1000.);
127	}else if(eta1 > -0.3 && eta1 < 0.3 && Pt1 < 0.6) {
128	efficTPC1 = effic4pol->Eval(Pt1*1000.);
129	}
130
131	// TOF efficiency
132
133	if(Pt1 > 0.2) {
134	if(TMath::Abs(eta1) > 0.7 && TMath::Abs(eta1) < 0.9 && Pt1 < 2.0) {
135	efficTOF1 = effic1polTOF->Eval(Pt1);
136	}else if(TMath::Abs(eta1) > 0.5 && TMath::Abs(eta1) < 0.7 && Pt1 < 2.4) {
137	efficTOF1 = effic2polTOF->Eval(Pt1);
138	}else if(TMath::Abs(eta1) > 0.3 && TMath::Abs(eta1) < 0.5 && Pt1 < 2.4) {
139	efficTOF1 = effic3polTOF->Eval(Pt1);
140	}else if(eta1 > -0.3 && eta1 < 0.3 && Pt1 < 2.6) {
141	efficTOF1 = effic4polTOF->Eval(Pt1);
142	}
143	}
144
145	Double_t rndmtpc=gRandom->Rndm();
146	Double_t rndmtof=gRandom->Rndm();
147	Double_t weightPID1=1.;
148	if(efficTPC1 < rndmtpc && efficTOF1 < rndmtof) { weightPID1=0.;}
149	ptKefftpc->Fill(Pt1,weightPID1);
150
151
152	Double_t Pt2=track2->Pt();
153	Double_t eta2=track2->Eta();
154
155	// TPC efficiency
156
157	if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && Pt2 < 0.48) {
158	efficTPC1 = effic1pol->Eval(Pt2*1000.);
159	}else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && Pt2 < 0.54) {
160	efficTPC1 = effic2pol->Eval(Pt2*1000.);
161	}else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && Pt2 < 0.585) {
162	efficTPC1 = effic3pol->Eval(Pt2*1000.);
163	}else if(eta2 > -0.3 && eta2 < 0.3 && Pt2 < 0.6) {
164	efficTPC1 = effic4pol->Eval(Pt2*1000.);
165	}
166
167	// TOF efficiency
168
169	if(Pt2 > 0.2) {
170	if(TMath::Abs(eta2) > 0.7 && TMath::Abs(eta2) < 0.9 && Pt2 < 2.0) {
171	efficTOF1 = effic1polTOF->Eval(Pt2);
172	}else if(TMath::Abs(eta2) > 0.5 && TMath::Abs(eta2) < 0.7 && Pt2 < 2.4) {
173	efficTOF1 = effic2polTOF->Eval(Pt2);
174	}else if(TMath::Abs(eta2) > 0.3 && TMath::Abs(eta2) < 0.5 && Pt2 < 2.4) {
175	efficTOF1 = effic3polTOF->Eval(Pt2);
176	}else if(eta2 > -0.3 && eta2 < 0.3 && Pt2 < 2.6) {
177	efficTOF1 = effic4polTOF->Eval(Pt2);
178	}
179	}
180
181	rndmtpc=gRandom->Rndm();
182	rndmtof=gRandom->Rndm();
183
184	Double_t weightPID2=1.;
185	if(efficTPC1 < rndmtpc && efficTOF1 < rndmtof) { weightPID2=0.;}
186
187	Double_t weightPID=weightPID1*weightPID2;
188	ptKefftpcboth->Fill(Pt1,weightPID);
189
190	return weightPID;
191	}