[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / UNICOR / AliUnicorAnalCorrel.cxx

/************************************************************************* 
* Copyright(c) 1998-2048, 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.                  * 
**************************************************************************/

// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005

//=============================================================================
// two-particle correlation analyzer
// Loop over pairs and fill pair histograms. The first particle is always 
// taken from ev0 and the second from ev1 so for ev0=ev1 one is getting true 
// pairs, otherwise mixed ones. 
//=============================================================================

#include <TROOT.h>
#include <TMath.h>
#include <TRandom2.h>
#include "AliUnicorEvent.h"
#include "AliUnicorHN.h"
#include "AliUnicorAnalCorrel.h"

ClassImp(AliUnicorAnalCorrel)
 
//=============================================================================
AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t ema, 
			 Int_t pid0, Int_t pid1, AnalysisFrame frame): 
  AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fZ0(0), fZ1(0), 
  fFrame(frame), fPa() {
  // constructor
  // emi and ema define the rapidity range for histogram

  TParticlePDG *part0 = AliUnicorAnal::fgPDG.GetParticle(fPid0);
  TParticlePDG *part1 = AliUnicorAnal::fgPDG.GetParticle(fPid1);
  fMass0 = part0? part0->Mass() : 0;
  fMass1 = part1? part1->Mass() : 0;
  fZ0 = part0? part0->Charge()/3.0 : 0;
  fZ1 = part1? part1->Charge()/3.0 : 0;
  double pi = TMath::Pi();

  // correlation function

  int nce = 6; double cebins[]={0,0.05,0.1,0.2,0.4,0.6,1.0};  // centrality bins

  //int npt = 7; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
  //int npt = 6; double ptbins[]={0,0.1,0.25,0.35,0.55,1.0,2.0}; // like Adam, except last bin split
  //int npt = 7; double ptbins[]={0,0.1,0.25,0.40,0.55,0.7,1.0,2.0}; // like Adam in Mar-2010, + first+last
  int npt = 10; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1.0,2.0};  // for Pb

  double qbins[200] = {0};
  //  for (int i=0;i<60;i++) qbins[i]=i*0.005;
  //  for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
  //  for (int i=0;i<30;i++) qbins[i]=i*0.010;
  for (int i=0;i<60;i++) qbins[i]=i*0.005;
  for (int i=0;i<35;i++) qbins[60+i]=0.3+i*0.02;
  for (int i=0;i<20;i++) qbins[95+i]=1.0+i*0.05;
  for (int i=0;i<11;i++) qbins[115+i]=2.0+i*0.20;

  TAxis *ax[8];
  ax[0] = new TAxis(4,-0.5,3.5);ax[0]->SetTitle("trumixrot");
  //  ax[1] = new TAxis(5,0,1.0);   ax[1]->SetTitle("centrality");
  ax[1] = new TAxis(nce,cebins);ax[1]->SetTitle("centrality");
  ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("y");          // pair y
  //  ax[3] = new TAxis(8,-pi,pi);  ax[3]->SetTitle("phi");      // wrt event plane
  ax[3] = new TAxis(1,-pi,pi);  ax[3]->SetTitle("phi");        // wrt event plane
  ax[4] = new TAxis(npt,ptbins);ax[4]->SetTitle("kt (GeV/c)"); // pair pt/2
  ax[5] = new TAxis(8,0,pi);    ax[5]->SetTitle("q-theta");
  ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
  ax[7] = new TAxis(125,qbins); ax[7]->SetTitle("q (GeV/c)");
  //  ax[7] = new TAxis(700,0,3.5); ax[7]->SetTitle("q (GeV/c)");
  AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
  fHistos.Add(pair);

  // correlation function bin monitor (needed to get <kt> etc.)

  TAxis *bx[3]={0};
  //bx[0] = new TAxis(*(pair->GetAxis(1))); // wait until root bug (516-00..527-06) fixed. For now, do: 
  bx[0] = new TAxis(ax[1]->GetNbins(), ax[1]->GetXmin(), ax[1]->GetXmax()); bx[0]->SetTitle("centrality"); 
  bx[1] = new TAxis(10*ax[2]->GetNbins(),emi,ema); bx[1]->SetTitle("y");         // pair y
  bx[2] = new TAxis(100,0,2); bx[2]->SetTitle("kt (GeV/c)");                     // pair pt/2
  AliUnicorHN *bimo = new AliUnicorHN("bimo",3,bx);
  for (int i=0; i<8; i++) delete ax[i];
  for (int i=0; i<3; i++) delete bx[i];
  fHistos.Add(bimo);

  // two-track resolution monitoring histogram

  ax[0] = new TAxis(3,-0.5,2.5);    ax[0]->SetTitle("trumixrot");
  ax[1] = new TAxis(2,-0.5,1.5);    ax[1]->SetTitle("cut applied");
  ax[2] = new TAxis(npt,ptbins);    ax[2]->SetTitle("(pair pt)/2 (GeV)");
  ax[3] = new TAxis(80,-0.08,0.08); ax[3]->SetTitle("dtheta");
  ax[4] = new TAxis(80,-0.20,0.20); ax[4]->SetTitle("dphi");
  AliUnicorHN *twot = new AliUnicorHN("twot",5,ax);
  for (int i=0; i<5; i++) delete ax[i];
  fHistos.Add(twot);

  gROOT->cd();
  printf("%s object named %s created\n",ClassName(),GetName());
}
//=============================================================================
void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot) 
{
  // process pairs from one or two (if mixing) events
  // tmr tells which histogram (bins) to fill: tru,mix,rot

  // Could be possibly accelerated by checking the "good particle" only once 
  // and caching the result. (Maybe the optimizer does it already.)

  static TRandom2 ran;
  AliUnicorHN *pair = (AliUnicorHN*) fHistos.At(0);
  AliUnicorHN *bimo = (AliUnicorHN*) fHistos.At(1);
  AliUnicorHN *twot = (AliUnicorHN*) fHistos.At(2);

  // mixing-and-rotating-proof centrality and reaction plane angle
  // (but not rotation-proof for rotation angles much different from 0 and 180)
  // true and rotated pairs are within the triangle (j<i), mixed - all
  // thus, proper rotation is either by 180, or by 170 AND 190, etc. 

  double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
  double q0x=0,q0y=0,q1x=0,q1y=0;
  ev0->RP(q0x,q0y);
  ev1->RP(q1x,q1y); 
  double rpphi = atan2(q0y+q1y,q0x+q1x);

  // loop over pairs 

  for (int i=0; i<ev0->NParticles(); i++) {
    if (!ev0->ParticleGood(i,fPid0)) continue;
    for (int j=0; j<ev1->NParticles(); j++) {
      if (ev0==ev1 && j<i && fPid0==fPid1 ) continue; 
      if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
      if (!ev1->ParticleGood(j,fPid1)) continue;
      fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i));
      fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
      if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
      twot->Fill((double) tmr, 0.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
      if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),fZ0,
      			 ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot,fZ1)) continue;
      twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
      fPa.CalcLAB(); // this could be organized better. AliUnicorPair should give k*?
      fPa.CalcPairCM();
      double qcm = fPa.QCM(); // momdif in pair cm - argument for Coulomb correction

      fPa.CalcLAB();
      if (fFrame == kPairFrame) fPa.CalcPairCM();
      if (fFrame == kLCMS)      fPa.CalcLcmsCM();
      if (fPa.QCM()==0) {printf("AliUnicorAnalCorrel: Q=0\n"); return;} // should not be too frequent
      double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
      double weigth = 1.0;
      /*
      static TH2D *coul = 0; 
      if (!coul) {
	TFile::Open("coulomb.root","read");
	coul = (TH2D*) gFile->Get("co");
	coul->SetDirectory(gROOT);
	gFile->Close();
      }
      if (tmr==0 && fPid0==fPid1) {
	double co = 0;
	if (qcm>0.999) co = 1;
	else if (qcm>0.001) co = coul->Interpolate(7,qcm);
	weigth = 1.0-0.5+0.5*co*(1+exp(-pow(fPa.QCM()*7/0.197,2))); 
      }
      */
      pair->Fill((double) tmr,           // 0 for tru, 1 for mix, 2 for rot
		 cent,                   // centrality
		 fPa.Rapidity(),         // pair rapidity
		 phi,                    // pair phi wrt reaction plane
		 fPa.Pt()/2.0,           // half of pair pt
		 fPa.QCMTheta(),         // polar angle of Q
		 fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
		 fPa.QCM(),              // |p2-p1| in c.m.s.
		 weigth);                // weigth
      if (tmr==0 && fPa.QCM()<0.2) bimo->Fill(cent, fPa.Rapidity(), fPa.Pt()/2.0, weigth);
      if (tmr==0) pair->Fill((double) 3,             // this is for Coulomb correction, maybe not necessary
			     cent,                   // centrality
			     fPa.Rapidity(),         // pair rapidity
			     phi,                    // pair phi wrt reaction plane
			     fPa.Pt()/2.0,           // half of pair pt
			     fPa.QCMTheta(),         // polar angle of Q
			     fPa.QCMPhiOut(),        // azimuthal angle of Q w.r.t. out
			     fPa.QCM(),              // |p2-p1| in c.m.s.
      			     weigth*qcm);            // weigth*Q
    }
  }
}
//=============================================================================
Commit	Line	Data
621688e4	1	/*************************************************************************
	2	* Copyright(c) 1998-2048, 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	// Author: Dariusz Miskowiec <mailto:d.miskowiec@gsi.de> 2005
	17
	18	//=============================================================================
	19	// two-particle correlation analyzer
	20	// Loop over pairs and fill pair histograms. The first particle is always
	21	// taken from ev0 and the second from ev1 so for ev0=ev1 one is getting true
	22	// pairs, otherwise mixed ones.
	23	//=============================================================================
	24
	25	#include <TROOT.h>
	26	#include <TMath.h>
	27	#include <TRandom2.h>
	28	#include "AliUnicorEvent.h"
	29	#include "AliUnicorHN.h"
	30	#include "AliUnicorAnalCorrel.h"
	31
	32	ClassImp(AliUnicorAnalCorrel)
	33
	34	//=============================================================================
61e4657c	35	AliUnicorAnalCorrel::AliUnicorAnalCorrel(const char *nam, Double_t emi, Double_t ema,
28eee19b	36	Int_t pid0, Int_t pid1, AnalysisFrame frame):
	37	AliUnicorAnal(nam), fPid0(pid0), fPid1(pid1), fMass0(0), fMass1(0), fZ0(0), fZ1(0),
	38	fFrame(frame), fPa() {
621688e4	39	// constructor
	40	// emi and ema define the rapidity range for histogram
	41
	42	TParticlePDG *part0 = AliUnicorAnal::fgPDG.GetParticle(fPid0);
	43	TParticlePDG *part1 = AliUnicorAnal::fgPDG.GetParticle(fPid1);
	44	fMass0 = part0? part0->Mass() : 0;
	45	fMass1 = part1? part1->Mass() : 0;
28eee19b	46	fZ0 = part0? part0->Charge()/3.0 : 0;
28eee19b	47	fZ1 = part1? part1->Charge()/3.0 : 0;
621688e4	48	double pi = TMath::Pi();
c6fc7f72	49
	50	// correlation function
	51
28eee19b	52	int nce = 6; double cebins[]={0,0.05,0.1,0.2,0.4,0.6,1.0}; // centrality bins
28eee19b	53
76d78859	54	//int npt = 7; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.7,1.0};
76d78859	55	//int npt = 6; double ptbins[]={0,0.1,0.25,0.35,0.55,1.0,2.0}; // like Adam, except last bin split
61e4657c	56	//int npt = 7; double ptbins[]={0,0.1,0.25,0.40,0.55,0.7,1.0,2.0}; // like Adam in Mar-2010, + first+last
28eee19b	57	int npt = 10; double ptbins[]={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1.0,2.0}; // for Pb
76d78859	58
61e4657c	59	double qbins[200] = {0};
28eee19b	60	// for (int i=0;i<60;i++) qbins[i]=i*0.005;
61e4657c	61	// for (int i=0;i<45;i++) qbins[20+i]=0.1+i*0.02;
28eee19b	62	// for (int i=0;i<30;i++) qbins[i]=i*0.010;
	63	for (int i=0;i<60;i++) qbins[i]=i*0.005;
	64	for (int i=0;i<35;i++) qbins[60+i]=0.3+i*0.02;
	65	for (int i=0;i<20;i++) qbins[95+i]=1.0+i*0.05;
	66	for (int i=0;i<11;i++) qbins[115+i]=2.0+i*0.20;
c6fc7f72	67
621688e4	68	TAxis *ax[8];
28eee19b	69	ax[0] = new TAxis(4,-0.5,3.5);ax[0]->SetTitle("trumixrot");
	70	// ax[1] = new TAxis(5,0,1.0); ax[1]->SetTitle("centrality");
	71	ax[1] = new TAxis(nce,cebins);ax[1]->SetTitle("centrality");
76d78859	72	ax[2] = new TAxis(3,emi,ema); ax[2]->SetTitle("y"); // pair y
	73	// ax[3] = new TAxis(8,-pi,pi); ax[3]->SetTitle("phi"); // wrt event plane
	74	ax[3] = new TAxis(1,-pi,pi); ax[3]->SetTitle("phi"); // wrt event plane
	75	ax[4] = new TAxis(npt,ptbins);ax[4]->SetTitle("kt (GeV/c)"); // pair pt/2
621688e4	76	ax[5] = new TAxis(8,0,pi); ax[5]->SetTitle("q-theta");
621688e4	77	ax[6] = new TAxis(16,-pi,pi); ax[6]->SetTitle("q-phi");
28eee19b	78	ax[7] = new TAxis(125,qbins); ax[7]->SetTitle("q (GeV/c)");
61e4657c	79	// ax[7] = new TAxis(700,0,3.5); ax[7]->SetTitle("q (GeV/c)");
621688e4	80	AliUnicorHN *pair = new AliUnicorHN("pair",8,ax);
621688e4	81	fHistos.Add(pair);
c6fc7f72	82
76d78859	83	// correlation function bin monitor (needed to get <kt> etc.)
76d78859	84
28eee19b	85	TAxis *bx[3]={0};
	86	//bx[0] = new TAxis(*(pair->GetAxis(1))); // wait until root bug (516-00..527-06) fixed. For now, do:
	87	bx[0] = new TAxis(ax[1]->GetNbins(), ax[1]->GetXmin(), ax[1]->GetXmax()); bx[0]->SetTitle("centrality");
	88	bx[1] = new TAxis(10*ax[2]->GetNbins(),emi,ema); bx[1]->SetTitle("y"); // pair y
	89	bx[2] = new TAxis(100,0,2); bx[2]->SetTitle("kt (GeV/c)"); // pair pt/2
	90	AliUnicorHN *bimo = new AliUnicorHN("bimo",3,bx);
	91	for (int i=0; i<8; i++) delete ax[i];
	92	for (int i=0; i<3; i++) delete bx[i];
76d78859	93	fHistos.Add(bimo);
76d78859	94
c6fc7f72	95	// two-track resolution monitoring histogram
	96
	97	ax[0] = new TAxis(3,-0.5,2.5); ax[0]->SetTitle("trumixrot");
	98	ax[1] = new TAxis(2,-0.5,1.5); ax[1]->SetTitle("cut applied");
76d78859	99	ax[2] = new TAxis(npt,ptbins); ax[2]->SetTitle("(pair pt)/2 (GeV)");
28eee19b	100	ax[3] = new TAxis(80,-0.08,0.08); ax[3]->SetTitle("dtheta");
28eee19b	101	ax[4] = new TAxis(80,-0.20,0.20); ax[4]->SetTitle("dphi");
c6fc7f72	102	AliUnicorHN *twot = new AliUnicorHN("twot",5,ax);
	103	for (int i=0; i<5; i++) delete ax[i];
	104	fHistos.Add(twot);
	105
621688e4	106	gROOT->cd();
	107	printf("%s object named %s created\n",ClassName(),GetName());
	108	}
	109	//=============================================================================
c6fc7f72	110	void AliUnicorAnalCorrel::Process(Int_t tmr, const AliUnicorEvent * const ev0, const AliUnicorEvent * const ev1, Double_t phirot)
621688e4	111	{
	112	// process pairs from one or two (if mixing) events
	113	// tmr tells which histogram (bins) to fill: tru,mix,rot
	114
	115	// Could be possibly accelerated by checking the "good particle" only once
	116	// and caching the result. (Maybe the optimizer does it already.)
	117
	118	static TRandom2 ran;
	119	AliUnicorHN pair = (AliUnicorHN) fHistos.At(0);
76d78859	120	AliUnicorHN bimo = (AliUnicorHN) fHistos.At(1);
76d78859	121	AliUnicorHN twot = (AliUnicorHN) fHistos.At(2);
621688e4	122
	123	// mixing-and-rotating-proof centrality and reaction plane angle
	124	// (but not rotation-proof for rotation angles much different from 0 and 180)
	125	// true and rotated pairs are within the triangle (j<i), mixed - all
	126	// thus, proper rotation is either by 180, or by 170 AND 190, etc.
	127
	128	double cent = (ev0->Centrality()+ev1->Centrality())/2.0;
61e4657c	129	double q0x=0,q0y=0,q1x=0,q1y=0;
621688e4	130	ev0->RP(q0x,q0y);
	131	ev1->RP(q1x,q1y);
	132	double rpphi = atan2(q0y+q1y,q0x+q1x);
	133
	134	// loop over pairs
	135
	136	for (int i=0; i<ev0->NParticles(); i++) {
	137	if (!ev0->ParticleGood(i,fPid0)) continue;
	138	for (int j=0; j<ev1->NParticles(); j++) {
	139	if (ev0==ev1 && j<i && fPid0==fPid1 ) continue;
	140	if (ev0==ev1 && j==i) continue; // beware, not even when rotated or non-identical
	141	if (!ev1->ParticleGood(j,fPid1)) continue;
621688e4	142	fPa.Set0(fMass0,ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i));
	143	fPa.Set1(fMass1,ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot);
	144	if (ev0==ev1 && fPid0==fPid1 && ran.Rndm()>=0.5) fPa.Swap();
c6fc7f72	145	twot->Fill((double) tmr, 0.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
28eee19b	146	if (!ev0->PairGood(ev0->ParticleP(i),ev0->ParticleTheta(i),ev0->ParticlePhi(i),fZ0,
28eee19b	147	ev1->ParticleP(j),ev1->ParticleTheta(j),ev1->ParticlePhi(j)+phirot,fZ1)) continue;
c6fc7f72	148	twot->Fill((double) tmr, 1.0, fPa.Pt()/2.0, fPa.DTheta(), fPa.DPhi(),1.0);
28eee19b	149	fPa.CalcLAB(); // this could be organized better. AliUnicorPair should give k*?
	150	fPa.CalcPairCM();
	151	double qcm = fPa.QCM(); // momdif in pair cm - argument for Coulomb correction
	152
621688e4	153	fPa.CalcLAB();
61e4657c	154	if (fFrame == kPairFrame) fPa.CalcPairCM();
	155	if (fFrame == kLCMS) fPa.CalcLcmsCM();
	156	if (fPa.QCM()==0) {printf("AliUnicorAnalCorrel: Q=0\n"); return;} // should not be too frequent
621688e4	157	double phi = TVector2::Phi_mpi_pi(fPa.Phi()-rpphi);
040da09c	158	double weigth = 1.0;
28eee19b	159	/*
	160	static TH2D *coul = 0;
	161	if (!coul) {
	162	TFile::Open("coulomb.root","read");
	163	coul = (TH2D*) gFile->Get("co");
	164	coul->SetDirectory(gROOT);
	165	gFile->Close();
	166	}
	167	if (tmr==0 && fPid0==fPid1) {
	168	double co = 0;
	169	if (qcm>0.999) co = 1;
	170	else if (qcm>0.001) co = coul->Interpolate(7,qcm);
	171	weigth = 1.0-0.5+0.5co(1+exp(-pow(fPa.QCM()*7/0.197,2)));
	172	}
	173	*/
c6fc7f72	174	pair->Fill((double) tmr, // 0 for tru, 1 for mix, 2 for rot
621688e4	175	cent, // centrality
	176	fPa.Rapidity(), // pair rapidity
	177	phi, // pair phi wrt reaction plane
	178	fPa.Pt()/2.0, // half of pair pt
	179	fPa.QCMTheta(), // polar angle of Q
	180	fPa.QCMPhiOut(), // azimuthal angle of Q w.r.t. out
	181	fPa.QCM(), // \|p2-p1\| in c.m.s.
040da09c	182	weigth); // weigth
28eee19b	183	if (tmr==0 && fPa.QCM()<0.2) bimo->Fill(cent, fPa.Rapidity(), fPa.Pt()/2.0, weigth);
	184	if (tmr==0) pair->Fill((double) 3, // this is for Coulomb correction, maybe not necessary
	185	cent, // centrality
	186	fPa.Rapidity(), // pair rapidity
	187	phi, // pair phi wrt reaction plane
	188	fPa.Pt()/2.0, // half of pair pt
	189	fPa.QCMTheta(), // polar angle of Q
	190	fPa.QCMPhiOut(), // azimuthal angle of Q w.r.t. out
	191	fPa.QCM(), // \|p2-p1\| in c.m.s.
	192	weigthqcm); // weigthQ
621688e4	193	}
	194	}
	195	}
	196	//=============================================================================