File indexing completed on 2024-09-27 07:03:27
0001
0002
0003
0004
0005
0006
0007
0008 #include "Kinematics.h"
0009 ClassImp(Kinematics)
0010
0011 Kinematics::Kinematics(
0012 Double_t enEleBeam,
0013 Double_t enIonBeam,
0014 Double_t crossAng
0015 )
0016 {
0017 srand(time(NULL));
0018
0019
0020 #ifdef SIDIS_MLPRED
0021 efnpackage = py::module_::import("testEFlowimport");
0022 pfnimport = efnpackage.attr("eflowPredict");
0023 #endif
0024
0025 IonMass = ProtonMass();
0026
0027
0028 crossAng *= 1e-3;
0029 crossAng = -1*TMath::Abs(crossAng);
0030
0031
0032
0033 mainFrame = fHeadOn;
0034
0035 qComponentsMethod = qQuadratic;
0036
0037
0038
0039
0040
0041
0042
0043 Double_t momEleBeam = EMtoP(enEleBeam,ElectronMass());
0044 Double_t momIonBeam = EMtoP(enIonBeam,IonMass);
0045 vecEleBeam.SetPxPyPzE(
0046 0,
0047 0,
0048 -momEleBeam,
0049 enEleBeam
0050 );
0051 vecIonBeam.SetPxPyPzE(
0052 momIonBeam * TMath::Sin(crossAng),
0053 0,
0054 momIonBeam * TMath::Cos(crossAng),
0055 enIonBeam
0056 );
0057 s = (vecEleBeam+vecIonBeam).M2();
0058
0059
0060
0061 BvecBoost = vecEleBeam + vecIonBeam;
0062 Bboost = -1*BvecBoost.BoostVector();
0063
0064 OvecBoost.SetXYZT( 0.0, 0.0, BvecBoost[2], BvecBoost[3] );
0065 Oboost = OvecBoost.BoostVector();
0066
0067 this->BoostToBeamComFrame(vecEleBeam,BvecEleBeam);
0068 this->BoostToBeamComFrame(vecIonBeam,BvecIonBeam);
0069
0070 rotAboutY = -TMath::ATan2( BvecIonBeam.Px(), BvecIonBeam.Pz() );
0071 BvecEleBeam.RotateY(rotAboutY);
0072 BvecIonBeam.RotateY(rotAboutY);
0073
0074 rotAboutX = TMath::ATan2( BvecIonBeam.Py(), BvecIonBeam.Pz() );
0075
0076
0077 tSpin = 1;
0078 lSpin = 1;
0079
0080
0081 polT = 0.80;
0082 polL = 0.;
0083 polBeam = 0.;
0084
0085
0086 RNG = std::make_unique<TRandomMixMax>(91874);
0087
0088
0089 countPIDsmeared=countPIDtrue=0;
0090 };
0091
0092
0093
0094
0095
0096
0097
0098 void Kinematics::GetQWNu_quadratic(){
0099
0100 double f,px,py,pz,pE;
0101 switch(mainFrame) {
0102 case fLab:
0103 f = y*(vecIonBeam.Dot(vecEleBeam));
0104 pz = vecIonBeam.Pz();
0105 py = vecIonBeam.Py();
0106 px = vecIonBeam.Px();
0107 pE = vecIonBeam.E();
0108 break;
0109 case fHeadOn:
0110 f = y*(HvecIonBeam.Dot(HvecEleBeam));
0111 pz = HvecIonBeam.Pz();
0112 py = HvecIonBeam.Py();
0113 px = HvecIonBeam.Px();
0114 pE = HvecIonBeam.E();
0115 break;
0116 };
0117 double hx = Pxh - px;
0118 double hy = Pyh - py;
0119
0120 double a = 1.0 - (pE*pE)/(pz*pz);
0121 double b = (2*pE/(pz*pz))*(px*hx + py*hy + f);
0122 double c = Q2 - hx*hx - hy*hy - (1/(pz*pz))*pow( (f+px*hx+py*hy) ,2.0);
0123 double disc = b*b - 4*a*c;
0124
0125 double qz1, qz2, qE1, qE2, qE, qz;
0126 if(disc>=0 && TMath::Abs(a)>1e-6) {
0127 qE1 = (-1*b+sqrt(b*b-4*a*c))/(2*a);
0128 qE2 = (-1*b-sqrt(b*b-4*a*c))/(2*a);
0129 qz1 = (-1*f + pE*qE1 - px*hx - py*hy)/(pz);
0130 qz2 = (-1*f + pE*qE2 - px*hx - py*hy)/(pz);
0131
0132 if(fabs(qE1) < fabs(qE2)){
0133 qE = qE1;
0134 qz = qz1;
0135 }
0136 else{
0137 qE = qE2;
0138 qz = qz2;
0139 }
0140
0141 vecQ.SetPxPyPzE(hx, hy, qz, qE);
0142 if(mainFrame==fHeadOn)
0143 this->TransformBackToLabFrame(vecQ,vecQ);
0144 vecW = vecIonBeam + vecQ;
0145 W = vecW.M();
0146 Nu = vecIonBeam.Dot(vecQ)/IonMass;
0147
0148 } else {
0149
0150 cerr << "ERROR: in Kinematics::GetQWNu_quadratic, ";
0151 if(isnan(disc)) cerr << "discriminant is NaN";
0152 else if(disc<0) cerr << "negative discriminant";
0153 else if(TMath::Abs(a)<1e-6) cerr << "zero denominator";
0154 else cerr << "unknown reason";
0155 cerr << "; skipping event" << endl;
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168 reconOK = false;
0169 }
0170 };
0171
0172
0173
0174
0175
0176 void Kinematics::GetQWNu_hadronic(){
0177 switch(mainFrame) {
0178 case fLab:
0179 vecQ = hadronSumVec - vecIonBeam;
0180 vecW = hadronSumVec;
0181 break;
0182 case fHeadOn:
0183 vecQ = hadronSumVec - HvecIonBeam;
0184 vecW = hadronSumVec;
0185 this->TransformBackToLabFrame(vecQ,vecQ);
0186 this->TransformBackToLabFrame(vecW,vecW);
0187 break;
0188 }
0189 W = vecW.M();
0190 Nu = vecIonBeam.Dot(vecQ)/IonMass;
0191 };
0192
0193
0194
0195
0196
0197 void Kinematics::GetQWNu_electronic(){
0198 vecQ = vecEleBeam - vecElectron;
0199 vecW = vecEleBeam + vecIonBeam - vecElectron;
0200 W = vecW.M();
0201 Nu = vecIonBeam.Dot(vecQ) / IonMass;
0202 };
0203
0204 void Kinematics::GetQWNu_ML(){
0205 hfsinfo.clear();
0206 float pidadj = 0;
0207 if(nHFS >= 2){
0208 std::vector<float> partHold;
0209 for(int i = 0; i < nHFS; i++){
0210 double pidsgn=(hfspid[i]/abs(hfspid[i]));
0211 if(abs(hfspid[i])==211) pidadj = 0.4*pidsgn;
0212 if(abs(hfspid[i])==22) pidadj = 0.2*pidsgn;
0213 if(abs(hfspid[i])==321) pidadj = 0.6*pidsgn;
0214 if(abs(hfspid[i])==2212) pidadj = 0.8*pidsgn;
0215 if(abs(hfspid[i])==11) pidadj = 1.0*pidsgn;
0216 partHold.push_back(hfspx[i]);
0217 partHold.push_back(hfspy[i]);
0218 partHold.push_back(hfspz[i]);
0219 partHold.push_back(hfsE[i]);
0220 partHold.push_back(pidadj);
0221 hfsinfo.push_back(partHold);
0222 partHold.clear();
0223 }
0224 double Q2ele, Q2DA, Q2JB;
0225 double xele, xDA, xJB;
0226 TLorentzVector vecQEle;
0227 globalinfo.clear();
0228 this->CalculateDISbyElectron();
0229 vecQEle.SetPxPyPzE(vecQ.Px(), vecQ.Py(), vecQ.Pz(), vecQ.E());
0230 Q2ele = Q2;
0231 xele = x;
0232 this->CalculateDISbyDA();
0233 Q2DA = Q2;
0234 xDA = x;
0235 this->CalculateDISbyJB();
0236 Q2JB = Q2;
0237 xJB = x;
0238 if( Q2DA > 0 && Q2DA < 1e4){
0239 globalinfo.push_back(log10(Q2DA));
0240 }
0241 else{
0242 globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0243 }
0244 if( Q2ele > 0 && Q2ele < 1e4){
0245 globalinfo.push_back(log10(Q2ele));
0246 }
0247 else{
0248 globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0249 }
0250 if( Q2JB > 0 && Q2JB < 1e4){
0251 globalinfo.push_back(log10(Q2JB));
0252 }
0253 else{
0254 globalinfo.push_back(log10((float) (rand()) / (float) (RAND_MAX/10000.0)));
0255 }
0256 if(xDA>0 && xDA < 1){
0257 globalinfo.push_back(-1*log10(xDA));
0258 }
0259 else{
0260 globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0) ));
0261 }
0262 if(xele>0 && xele < 1){
0263 globalinfo.push_back(-1*log10(xele));
0264 }
0265 else{
0266 globalinfo.push_back(-1*log10( (float) (rand()) / (float) (RAND_MAX/1.0) ));
0267 }
0268 if(xJB>0 && xJB < 1){
0269 globalinfo.push_back(-1*log10(xJB));
0270 }
0271 else{
0272 globalinfo.push_back( -1*log10((float) (rand()) / (float) (RAND_MAX/1.0) ) );
0273 }
0274 globalinfo.push_back(vecQEle.Px());
0275 globalinfo.push_back(vecQEle.Py());
0276 globalinfo.push_back(vecQEle.Pz());
0277 globalinfo.push_back(vecQEle.E());
0278 #ifdef SIDIS_MLPRED
0279 py::object nnoutput = pfnimport(hfsinfo, globalinfo);
0280 std::vector<float> nnvecq = nnoutput.cast<std::vector<float>>();
0281 vecQ.SetPxPyPzE(nnvecq[0],nnvecq[1],nnvecq[2],nnvecq[3]);
0282 #endif
0283 }
0284 else{
0285 this->CalculateDISbyElectron();
0286 }
0287 vecW = vecEleBeam + vecIonBeam - vecElectron;
0288 W = vecW.M();
0289 Nu = vecIonBeam.Dot(vecQ) / IonMass;
0290
0291 }
0292
0293
0294
0295
0296
0297 Bool_t Kinematics::CalculateDIS(TString recmethod){
0298
0299 reconOK = true;
0300
0301
0302
0303 this->TransformToHeadOnFrame(vecEleBeam,HvecEleBeam);
0304 this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0305 this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0306
0307
0308 if (recmethod.CompareTo( "Ele", TString::kIgnoreCase)==0) { this->CalculateDISbyElectron(); }
0309 else if(recmethod.CompareTo( "DA", TString::kIgnoreCase)==0) { this->CalculateDISbyDA(); }
0310 else if(recmethod.CompareTo( "JB", TString::kIgnoreCase)==0) { this->CalculateDISbyJB(); }
0311 else if(recmethod.CompareTo( "Mixed", TString::kIgnoreCase)==0) { this->CalculateDISbyMixed(); }
0312 else if(recmethod.CompareTo( "Sigma", TString::kIgnoreCase)==0) { this->CalculateDISbySigma(); }
0313 else if(recmethod.CompareTo( "eSigma", TString::kIgnoreCase)==0) { this->CalculateDISbyeSigma(); }
0314 else if(recmethod.CompareTo( "ML", TString::kIgnoreCase)==0) { this->CalculateDISbyML(); }
0315 else {
0316 cerr << "ERROR: unknown reconstruction method" << endl;
0317 return false;
0318 };
0319
0320
0321
0322 CvecBoost = vecQ + vecIonBeam;
0323 Cboost = -1*CvecBoost.BoostVector();
0324
0325 IvecBoost = vecIonBeam;
0326 Iboost = -1*IvecBoost.BoostVector();
0327
0328
0329
0330
0331 gamma = 2*ProtonMass()*x / TMath::Sqrt(Q2);
0332 epsilon = ( 1 - y - TMath::Power(gamma*y,2)/4 ) /
0333 ( 1 - y + y*y/2 + TMath::Power(gamma*y,2)/4 );
0334
0335 depolA = y*y / (2 - 2*epsilon);
0336 depolB = depolA * epsilon;
0337 depolC = depolA * TMath::Sqrt(1-epsilon*epsilon);
0338 depolV = depolA * TMath::Sqrt(2*epsilon*(1+epsilon));
0339 depolW = depolA * TMath::Sqrt(2*epsilon*(1-epsilon));
0340
0341 if(depolA==0) depolP1=depolP2=depolP3=depolP4=0;
0342 else {
0343 depolP1 = depolB / depolA;
0344 depolP2 = depolC / depolA;
0345 depolP3 = depolV / depolA;
0346 depolP4 = depolW / depolA;
0347 };
0348
0349 return reconOK;
0350 };
0351
0352
0353
0354 void Kinematics::CalculateDISbyElectron() {
0355 this->GetQWNu_electronic();
0356 Q2 = -1 * vecQ.M2();
0357 x = Q2 / ( 2 * vecQ.Dot(vecIonBeam) );
0358 y = vecIonBeam.Dot(vecQ) / vecIonBeam.Dot(vecEleBeam);
0359 };
0360
0361
0362 void Kinematics::CalculateDISbyJB(){
0363 switch(mainFrame) {
0364 case fLab: y = sigmah/(2*vecEleBeam.E()); break;
0365 case fHeadOn: y = sigmah/(2*HvecEleBeam.E()); break;
0366 };
0367 Q2 = (Pxh*Pxh + Pyh*Pyh)/(1-y);
0368 x = Q2/(s*y);
0369 switch(qComponentsMethod) {
0370 case qQuadratic: this->GetQWNu_quadratic(); break;
0371 case qHadronic: this->GetQWNu_hadronic(); break;
0372 case qElectronic: this->GetQWNu_electronic(); break;
0373 }
0374 };
0375
0376
0377
0378 void Kinematics::CalculateDISbyDA(){
0379 float thetah = acos( (Pxh*Pxh+Pyh*Pyh - sigmah*sigmah)/(Pxh*Pxh+Pyh*Pyh+sigmah*sigmah) );
0380 float thetae;
0381 switch(mainFrame) {
0382 case fLab:
0383 thetae = vecElectron.Theta();
0384 Q2 = 4.0*vecEleBeam.E()*vecEleBeam.E()*sin(thetah)*(1+cos(thetae))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0385 break;
0386 case fHeadOn:
0387 thetae = HvecElectron.Theta();
0388 Q2 = 4.0*HvecEleBeam.E()*HvecEleBeam.E()*sin(thetah)*(1+cos(thetae))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0389 break;
0390 };
0391 y = (sin(thetae)*(1-cos(thetah)))/(sin(thetah)+sin(thetae)-sin(thetah+thetae));
0392 x = Q2/(s*y);
0393 switch(qComponentsMethod) {
0394 case qQuadratic: this->GetQWNu_quadratic(); break;
0395 case qHadronic: this->GetQWNu_hadronic(); break;
0396 case qElectronic: this->GetQWNu_electronic(); break;
0397 }
0398 };
0399
0400
0401
0402 void Kinematics::CalculateDISbyMixed(){
0403 this->GetQWNu_electronic();
0404 Q2 = -1*vecQ.M2();
0405 switch(mainFrame) {
0406 case fLab: y = sigmah/(2*vecEleBeam.E()); break;
0407 case fHeadOn: y = sigmah/(2*HvecEleBeam.E()); break;
0408 }
0409 x = Q2/(s*y);
0410 };
0411
0412
0413
0414 void Kinematics::CalculateDISbySigma(){
0415 switch(mainFrame) {
0416 case fLab:
0417 y = sigmah/(sigmah + vecElectron.E()*(1-cos(vecElectron.Theta())));
0418 Q2 = (vecElectron.Px()*vecElectron.Px() + vecElectron.Py()*vecElectron.Py())/(1-y);
0419 break;
0420 case fHeadOn:
0421 y = sigmah/(sigmah + HvecElectron.E()*(1-cos(HvecElectron.Theta())));
0422 Q2 = (HvecElectron.Px()*HvecElectron.Px() + HvecElectron.Py()*HvecElectron.Py())/(1-y);
0423 break;
0424 }
0425 x = Q2/(s*y);
0426 switch(qComponentsMethod) {
0427 case qQuadratic: this->GetQWNu_quadratic(); break;
0428 case qHadronic: this->GetQWNu_hadronic(); break;
0429 case qElectronic: this->GetQWNu_electronic(); break;
0430 }
0431 };
0432
0433
0434
0435 void Kinematics::CalculateDISbyeSigma(){
0436 this->CalculateDISbySigma();
0437 Double_t xsigma = x;
0438 Double_t ysigma = y;
0439 vecQ = vecEleBeam - vecElectron;
0440 Q2 = -1 * vecQ.M2();
0441 y = Q2/(s*xsigma);
0442 x = xsigma;
0443 switch(qComponentsMethod) {
0444 case qQuadratic: this->GetQWNu_quadratic(); break;
0445 case qHadronic: this->GetQWNu_hadronic(); break;
0446 case qElectronic: this->GetQWNu_electronic(); break;
0447 }
0448 };
0449
0450 void Kinematics::CalculateDISbyML() {
0451 this->GetQWNu_ML();
0452 Q2 = -1 * vecQ.M2();
0453 x = Q2 / ( 2 * vecQ.Dot(vecIonBeam) );
0454 y = vecIonBeam.Dot(vecQ) / vecIonBeam.Dot(vecEleBeam);
0455 };
0456
0457
0458
0459
0460
0461 void Kinematics::CalculateHadronKinematics() {
0462
0463 pLab = vecHadron.P();
0464 pTlab = vecHadron.Pt();
0465 phiLab = vecHadron.Phi();
0466 etaLab = vecHadron.Eta();
0467
0468 z = vecIonBeam.Dot(vecHadron) / vecIonBeam.Dot(vecQ);
0469
0470 mX = (vecW-vecHadron).M();
0471
0472 this->BoostToComFrame(vecHadron,CvecHadron);
0473 this->BoostToComFrame(vecQ,CvecQ);
0474 this->BoostToIonFrame(vecHadron,IvecHadron);
0475 this->BoostToIonFrame(vecQ,IvecQ);
0476 this->BoostToIonFrame(vecElectron,IvecElectron);
0477
0478 xF = 2 * CvecHadron.Vect().Dot(CvecQ.Vect()) /
0479 (W * CvecQ.Vect().Mag());
0480
0481 phiH = AdjAngle(PlaneAngle(
0482 IvecQ.Vect(), IvecElectron.Vect(),
0483 IvecQ.Vect(), IvecHadron.Vect()
0484 ));
0485
0486 tSpin = RNG->Uniform() < 0.5 ? 1 : -1;
0487 lSpin = RNG->Uniform() < 0.5 ? 1 : -1;
0488 vecSpin.SetXYZT(0,1,0,0);
0489 this->BoostToIonFrame(vecSpin,IvecSpin);
0490 phiS = AdjAngle(PlaneAngle(
0491 IvecQ.Vect(), IvecElectron.Vect(),
0492 IvecQ.Vect(), IvecSpin.Vect()
0493 ));
0494
0495 pT = Reject(
0496 IvecHadron.Vect(),
0497 IvecQ.Vect()
0498 ).Mag();
0499
0500 qT = pT / z;
0501 };
0502
0503
0504 void Kinematics::ValidateHeadOnFrame() {
0505 this->BoostToIonFrame(vecEleBeam,IvecEleBeam);
0506 this->BoostToIonFrame(vecIonBeam,IvecIonBeam);
0507 this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0508 this->TransformToHeadOnFrame(vecEleBeam,HvecEleBeam);
0509 this->TransformToHeadOnFrame(vecIonBeam,HvecIonBeam);
0510 this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0511 this->TransformToHeadOnFrame(vecHadron,HvecHadron);
0512 printf("\nVALIDATION:\n");
0513 printf("lab E: "); vecEleBeam.Print();
0514 printf("lab I: "); vecIonBeam.Print();
0515 printf("ion RF E: "); IvecEleBeam.Print();
0516 printf("ion RF I: "); IvecIonBeam.Print();
0517 printf("head-on E: "); HvecEleBeam.Print();
0518 printf("head-on I: "); HvecIonBeam.Print();
0519 printf("---\n");
0520 printf("lab electron: "); vecElectron.Print();
0521 printf("head-on electron: "); HvecElectron.Print();
0522 printf("difference: "); (vecElectron-HvecElectron).Print();
0523 printf("---\n");
0524 printf("lab hadron: "); vecHadron.Print();
0525 printf("head-on hadron: "); HvecHadron.Print();
0526 printf("difference: "); (vecHadron-HvecHadron).Print();
0527 };
0528
0529
0530
0531 void Kinematics::AddToHFS(TLorentzVector p4_) {
0532 TLorentzVector p4 = p4_;
0533 if(mainFrame==fHeadOn) this->TransformToHeadOnFrame(p4,p4);
0534 sigmah += (p4.E() - p4.Pz());
0535 Pxh += p4.Px();
0536 Pyh += p4.Py();
0537 hadronSumVec += p4;
0538 countHadrons++;
0539 };
0540
0541
0542
0543 void Kinematics::AddToHFSTree(TLorentzVector p4, int pid) {
0544 hfspx.push_back(p4.Px());
0545 hfspy.push_back(p4.Py());
0546 hfspz.push_back(p4.Pz());
0547 hfsE.push_back(p4.E());
0548 hfspid.push_back(pid);
0549
0550 nHFS++;
0551 };
0552
0553
0554
0555 void Kinematics::AddTrackToHFSTree(TLorentzVector p4, int pid) {
0556 selectedHadPx.push_back(p4.Px());
0557 selectedHadPy.push_back(p4.Py());
0558 selectedHadPz.push_back(p4.Pz());
0559 selectedHadE.push_back(p4.E());
0560 selectedHadPID.push_back(pid);
0561 };
0562
0563
0564 void Kinematics::SubtractElectronFromHFS() {
0565 if(!isnan(vecElectron.E())){
0566 switch(mainFrame) {
0567 case fLab:
0568 sigmah -= (vecElectron.E() - vecElectron.Pz());
0569 Pxh -= vecElectron.Px();
0570 Pyh -= vecElectron.Py();
0571 hadronSumVec -= vecElectron;
0572 break;
0573 case fHeadOn:
0574 this->TransformToHeadOnFrame(vecElectron,HvecElectron);
0575 sigmah -= (HvecElectron.E() - HvecElectron.Pz());
0576 Pxh -= HvecElectron.Px();
0577 Pyh -= HvecElectron.Py();
0578 hadronSumVec -= HvecElectron;
0579 break;
0580 }
0581 countHadrons--;
0582 } else {
0583 cerr << "ERROR: electron energy is NaN" << endl;
0584
0585 }
0586 };
0587
0588
0589
0590 void Kinematics::ResetHFS() {
0591 sigmah = Pxh = Pyh = 0;
0592 hadronSumVec.SetPxPyPzE(0,0,0,0);
0593 countHadrons = 0;
0594 nHFS = 0;
0595
0596 hfspx.clear();
0597 hfspy.clear();
0598 hfspz.clear();
0599 hfsE.clear();
0600 hfspid.clear();
0601
0602 selectedHadPx.clear();
0603 selectedHadPy.clear();
0604 selectedHadPz.clear();
0605 selectedHadE.clear();
0606 selectedHadPID.clear();
0607 };
0608
0609
0610
0611 #ifndef EXCLUDE_DELPHES
0612
0613
0614
0615 void Kinematics::GetHFS(
0616 TObjArrayIter itTrack,
0617 TObjArrayIter itEFlowTrack,
0618 TObjArrayIter itEFlowPhoton,
0619 TObjArrayIter itEFlowNeutralHadron,
0620 TObjArrayIter itpfRICHTrack,
0621 TObjArrayIter itDIRCepidTrack, TObjArrayIter itDIRChpidTrack,
0622 TObjArrayIter itBTOFepidTrack, TObjArrayIter itBTOFhpidTrack,
0623 TObjArrayIter itdualRICHagTrack, TObjArrayIter itdualRICHcfTrack
0624 ) {
0625
0626
0627 this->ResetHFS();
0628 itTrack.Reset();
0629 itEFlowTrack.Reset();
0630 itEFlowPhoton.Reset();
0631 itEFlowNeutralHadron.Reset();
0632
0633
0634 while(Track *track = (Track*)itTrack() ){
0635 TLorentzVector trackp4 = track->P4();
0636 if(!isnan(trackp4.E())){
0637 if( std::abs(track->Eta) < 4.0 ){
0638
0639 int pid = GetTrackPID(
0640 track,
0641 itpfRICHTrack,
0642 itDIRCepidTrack, itDIRChpidTrack,
0643 itBTOFepidTrack, itBTOFhpidTrack,
0644 itdualRICHagTrack, itdualRICHcfTrack
0645 );
0646
0647 if(pid != -1){
0648 trackp4.SetPtEtaPhiM(trackp4.Pt(),trackp4.Eta(),trackp4.Phi(),correctMass(pid));
0649 countPIDsmeared++;
0650 }
0651 else {
0652 countPIDtrue++;
0653
0654 }
0655
0656 this->AddToHFS(trackp4);
0657 this->AddToHFSTree(trackp4,pid);
0658 }
0659 }
0660 }
0661
0662
0663 while(Track *eflowTrack = (Track*)itEFlowTrack() ){
0664 TLorentzVector eflowTrackp4 = eflowTrack->P4();
0665 int pid = GetTrackPID(
0666 eflowTrack,
0667 itpfRICHTrack,
0668 itDIRCepidTrack, itDIRChpidTrack,
0669 itBTOFepidTrack, itBTOFhpidTrack,
0670 itdualRICHagTrack, itdualRICHcfTrack
0671 );
0672
0673 if(!isnan(eflowTrackp4.E())){
0674 if(std::abs(eflowTrack->Eta) >= 4.0){
0675 this->AddToHFS(eflowTrackp4);
0676 this->AddToHFSTree(eflowTrackp4, pid);
0677 }
0678 }
0679 }
0680
0681
0682 while(Tower* towerPhoton = (Tower*)itEFlowPhoton() ){
0683 TLorentzVector towerPhotonp4 = towerPhoton->P4();
0684 if(!isnan(towerPhotonp4.E())){
0685 if( std::abs(towerPhoton->Eta) < 4.0 ){
0686 this->AddToHFS(towerPhotonp4);
0687 this->AddToHFSTree(towerPhotonp4,22);
0688 }
0689 }
0690 }
0691
0692
0693 while(Tower* towerNeutralHadron = (Tower*)itEFlowNeutralHadron() ){
0694 TLorentzVector towerNeutralHadronp4 = towerNeutralHadron->P4();
0695 if(!isnan(towerNeutralHadronp4.E())){
0696 if( std::abs(towerNeutralHadron->Eta) < 4.0 ){
0697 this->AddToHFS(towerNeutralHadronp4);
0698 this->AddToHFSTree(towerNeutralHadronp4,-1);
0699 }
0700 }
0701 }
0702
0703
0704 this->SubtractElectronFromHFS();
0705 };
0706
0707
0708
0709 void Kinematics::GetTrueHFS(TObjArrayIter itParticle){
0710
0711
0712 this->ResetHFS();
0713 itParticle.Reset();
0714
0715
0716 while(GenParticle *partTrue = (GenParticle*)itParticle() ) {
0717 if(partTrue->Status == 1) this->AddToHFS(partTrue->P4());
0718 }
0719
0720
0721 this->SubtractElectronFromHFS();
0722 };
0723
0724
0725
0726 int Kinematics::GetTrackPID(
0727 Track *track,
0728 TObjArrayIter itpfRICHTrack,
0729 TObjArrayIter itDIRCepidTrack, TObjArrayIter itDIRChpidTrack,
0730 TObjArrayIter itBTOFepidTrack, TObjArrayIter itBTOFhpidTrack,
0731 TObjArrayIter itdualRICHagTrack, TObjArrayIter itdualRICHcfTrack
0732 ) {
0733
0734 itpfRICHTrack.Reset();
0735 itDIRCepidTrack.Reset(); itDIRChpidTrack.Reset();
0736 itBTOFepidTrack.Reset(); itBTOFhpidTrack.Reset();
0737 itdualRICHagTrack.Reset(); itdualRICHcfTrack.Reset();
0738 GenParticle *trackParticle = (GenParticle*)track->Particle.GetObject();
0739 GenParticle *detectorParticle;
0740
0741
0742
0743 while(Track *detectorTrack = (Track*)itpfRICHTrack() ){
0744 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0745 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0746 }
0747
0748 while(Track *detectorTrack = (Track*)itDIRCepidTrack() ){
0749 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0750 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0751 }
0752 while(Track *detectorTrack = (Track*)itDIRChpidTrack() ){
0753 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0754 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0755 }
0756
0757 while(Track *detectorTrack = (Track*)itBTOFepidTrack() ){
0758 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0759 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0760 }
0761 while(Track *detectorTrack = (Track*)itBTOFhpidTrack() ){
0762 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0763 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0764 }
0765
0766 while(Track *detectorTrack = (Track*)itdualRICHagTrack() ){
0767 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0768 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0769 }
0770 while(Track *detectorTrack = (Track*)itdualRICHcfTrack() ){
0771 detectorParticle = (GenParticle*)detectorTrack->Particle.GetObject();
0772 if( detectorParticle == trackParticle ) return detectorTrack->PID;
0773 }
0774
0775 return -1;
0776 }
0777
0778 #endif
0779
0780
0781
0782
0783
0784
0785
0786 void Kinematics::BoostToComFrame(TLorentzVector Lvec, TLorentzVector &Cvec) {
0787 Cvec=Lvec;
0788 Cvec.Boost(Cboost);
0789 };
0790
0791
0792 void Kinematics::BoostToIonFrame(TLorentzVector Lvec, TLorentzVector &Ivec) {
0793 Ivec=Lvec;
0794 Ivec.Boost(Iboost);
0795 };
0796
0797
0798 void Kinematics::BoostToBeamComFrame(TLorentzVector Lvec, TLorentzVector &Bvec) {
0799 Bvec=Lvec;
0800 Bvec.Boost(Bboost);
0801 };
0802
0803
0804 void Kinematics::TransformToHeadOnFrame(TLorentzVector Lvec, TLorentzVector &Hvec) {
0805 this->BoostToBeamComFrame(Lvec,Hvec);
0806 Hvec.RotateY(rotAboutY);
0807 Hvec.RotateX(rotAboutX);
0808 Hvec.Boost(Oboost);
0809 };
0810
0811
0812 void Kinematics::TransformBackToLabFrame(TLorentzVector Hvec, TLorentzVector &Lvec) {
0813 Lvec=Hvec;
0814 Lvec.Boost(-1*Oboost);
0815 Lvec.RotateX(-rotAboutX);
0816 Lvec.RotateY(-rotAboutY);
0817 Lvec.Boost(-1*Bboost);
0818 };
0819
0820
0821
0822
0823 void Kinematics::InjectFakeAsymmetry() {
0824
0825 moduVal[0] = TMath::Sin(phiH-phiS);
0826 moduVal[1] = TMath::Sin(phiH+phiS) * depolP1;
0827
0828 ampVal[0] = 0.1;
0829 ampVal[1] = 0.1;
0830
0831 asymInject = 0;
0832 asymInject += ampVal[0]/0.2 * x * moduVal[0];
0833 asymInject += -ampVal[1]/0.2 * x * moduVal[1];
0834
0835
0836 asymInject *= polT;
0837
0838 RN = RNG->Uniform();
0839 tSpin = (RN<0.5*(1+asymInject)) ? 1 : -1;
0840 };
0841
0842
0843 Kinematics::~Kinematics() {
0844 };
0845