#include <iostream>
using namespace std;
//#include <unistd.h>   // for sleep()
//#include <stdlib.h>   // for atoi()
//#include <math.h>

#include "T49Dedx.h"

//*-*-* ROOT includes

#include "TROOT.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TH1.h"
#include "TF1.h"
#include "TF2.h"
#include "TProfile.h"
#include "TKey.h"
#include "TTree.h"
#include "TDirectory.h"
#include "TObjectTable.h"
#include "TClassTable.h"
#include "TMinuit.h"

#include "Rtypes.h"

#include "PhysConst.h"


////////////////////////////////////////////////////////////////////    

Float_t  RelRise(Float_t p, Float_t mass);
Double_t FuncDedx(Double_t *x,Double_t *par); 
void     TestFunc(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag); 

////////////////////////////////////////////////////////////////////    

#define MINPOINTS       25
#define MINDEDXPOINTS   50

#define MAXPART  15000
Float_t   val[MAXPART][2]; 
Int_t     NumPart;

Double_t Resolution = 0.04;
Double_t MinIon     = 326.0;
Double_t Norm       = 1.0/sqrt(2.0*3.14159);


TF2 *F2Dedx;

////////////////////////////////////////////////////////////////////

ClassImp(T49Dedx)

/*
$Log: T49Dedx.C,v $
Revision 1.4  2005/05/20 11:26:47  cblume
Get rid of deprecated warnings

Revision 1.3  1999/12/01 09:39:25  cblume
Changed to root 2.23-09. Updated the makefiles. Changes in T49CutTrack, T49Dedx, T49Proj1, T49Proj2, TRootDS

 * Revision 1.2  1999/11/24  16:03:44  cblume
 * Fixes in gteval_nt_root.h and T49Mixer. Addrd T49Dedx again in T49ANA. Some cosmetics.
 *
 * Revision 1.1  1999/11/24  15:33:44  cblume
 * Add T49Dedx.C
 *
*/

 T49Dedx::T49Dedx()
{
    fRes0 = 0.23;
    fRes1 = 0.0003;
    fMinIonPos = 303.5;
    fMinIonNeg = 298.5;
    sprintf(fCalibHome,"/home/data1");

}

///////////////////////// init //////////////////////

 Int_t T49Dedx::Init() 
{
    cerr << "T49Dedx::init() has been called" << endl;

    F2Dedx   = new TF2("F2Dedx",FuncDedx,4);
    F2Dedx->SetName("F2Dedx");
    F2Dedx->SetParameters(0.1,0.1,0.1,100);

    cerr << "T49Dedx::init() has finished" << endl;

    return kTRUE;
}

 TF2 *T49Dedx::GetFunc()
{
    return F2Dedx;
}

#define DEDXNORM 0.6245

 Float_t T49Dedx::GetRelRise(Float_t mom, Float_t mass)
{
    return RelRise(mom,mass);
}

Double_t FuncDedx(Double_t *x,Double_t *par)
{
    Int_t i;
    Double_t pos[4];
    Double_t sig[4];
    Double_t res;
    Double_t f;
    Double_t PiConst;

#ifdef DEBUGFUNC
    printf("FuncDedx called : %g %g,  %g %g %g %g\n",x[0],x[1],par[0],par[1],par[2],par[3]);
#endif

    for(i=0;i<3;i++)
    {
       if(par[i] < 0)
       {  
           return 100000000.0;
       }
    }
    PiConst = 1.0 - (par[0]+par[1]+par[2]);
    if(PiConst < 0)
      return 10000000.0;
      
    res = Resolution + 0.001 * x[0];
    pos[0] = MinIon * RelRise((Float_t) x[0],(Float_t) kPMASS);
    pos[1] = MinIon * RelRise((Float_t) x[0],(Float_t) kKMASS);
    pos[2] = MinIon * RelRise((Float_t) x[0],(Float_t) kPIMASS);
    pos[3] = MinIon * RelRise((Float_t) x[0],(Float_t) kEMASS);
    sig[0] = pos[0] * res;
    sig[1] = pos[1] * res;
    sig[2] = pos[2] * res;
    sig[3] = pos[3] * res;
    f = par[3]*par[0] *exp(-(x[1]-pos[0])*(x[1]-pos[0])/(2.0*sig[0]*sig[0])) +
        par[3]*par[1] *exp(-(x[1]-pos[1])*(x[1]-pos[1])/(2.0*sig[1]*sig[1])) +
        par[3]*PiConst*exp(-(x[1]-pos[2])*(x[1]-pos[2])/(2.0*sig[2]*sig[2])) +
        par[3]*par[2] *exp(-(x[1]-pos[3])*(x[1]-pos[3])/(2.0*sig[3]*sig[3]));

#ifdef DEBUGFUNC
    printf("f = %g\n",f);
#endif

    return f;
}

 void T49Dedx::EventFit(T49EventRoot *event, Double_t *p2pi, Double_t *k2pi)
{
    //Int_t     j;
    //Int_t     i;
    Int_t     ErrFlag;
    //Float_t  *row;
    Double_t arglist[10];
    Double_t we, al, bl;
    Double_t pNum = 0;
    Double_t kNum = 0;
    Double_t par;
    const Double_t Pmin = 5;
    const Double_t Pmax = 50;
    T49Particle  *track;

    //Double_t kPar[5] = {1.0443E2,2.0724E1,-4.067,-5.349E-1,-1.3876E-2};
    //Double_t pPar[5] = {2.46E2,-9.706E-1,3.1679E-1,0,0};
    //Double_t piPar[5] ={-1.3874E3,1.68399E3,-1.73976E2,6.25518,-7.50058E-2};
    //Double_t pNorm  = 15695.7; 
    //Double_t kNorm  = 30905.2; 
    //Double_t piNorm = 105105.2; 
    
    static Int_t index = 0;
    static TMinuit *myMinuit = NULL;

    NumPart = 0;
    //T49ParticleIter *iter = event->GetParticleIter();
    //while((track = (T49Particle *) iter->Next()))
    //{
    TClonesArray *ParticleList = event->GetPrimaryParticles();
    for (Int_t iTrack = 0; iTrack < ParticleList->GetEntries(); iTrack++) {
        track = (T49Particle *) ParticleList->UncheckedAt(iTrack);
        if(track->GetP() > Pmin && 
           track->GetP() < Pmax &&
           track->GetTmeanCharge() > 300 && 
           track->GetTmeanCharge() < 500 && 
           track->GetNFitPoint() > MINDEDXPOINTS)
        {
            val[NumPart][0] = track->GetP();
            val[NumPart][1] = track->GetTmeanCharge();
            NumPart++;
        }
    }
    index += MAXPART;
    printf("Number of particles is %dn",NumPart);

////// Check if Minuit is initialized /////////////////////////////


   printf("Create myMinuitn");
   if(myMinuit == NULL)
   {
       myMinuit = new TMinuit();
       printf("Init Minuitn");
       myMinuit->mninit(5,6,7);
       printf("Call SetFCNn");
       myMinuit->SetFCN(TestFunc);
    }

   printf("Start event by event dedx fitn");
   myMinuit->mncler();
   myMinuit->mnparm(1,"rypr",0.1,0.01,0.0,1.0,ErrFlag);
   myMinuit->mnparm(2,"ryka",0.1,0.01,0.0,1.0,ErrFlag);
   arglist[0] = 0.5;
   myMinuit->mnexcm("SET ERR",arglist,1,ErrFlag);

   arglist[0] = 2;
   myMinuit->mnexcm("SET STR",arglist,1,ErrFlag);

   arglist[0] = 1;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   arglist[0] = 0;
   myMinuit->mnexcm("MINI",arglist,0,ErrFlag);

   arglist[0] = 3;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   arglist[0] = 6;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   TString str;
   myMinuit->mnpout(1,str,par,we,al,bl,ErrFlag);
   //myMinuit->mnpout(1,"rypr",par,we,al,bl,ErrFlag);
   *p2pi = par;
   myMinuit->mnpout(2,str,par,we,al,bl,ErrFlag);
   //myMinuit->mnpout(2,"ryka",par,we,al,bl,ErrFlag);
   *k2pi = par;

   myMinuit->mnexcm("CLEAR",arglist,0,ErrFlag);
 
   printf("Finished event by event dedx fit, found %g %gn",pNum,kNum);

}

 Double_t T49Dedx::GetMinIon()
{
   return MinIon;
}

 Double_t T49Dedx::GetResolution()
{
   return Resolution;
}

 void T49Dedx::SetMinIon(Double_t x)
{
   MinIon = x;
}

 void T49Dedx::SetResolution(Double_t x)
{
   Resolution = x;
}

void TestFunc(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
    Int_t i;
    Int_t j;
    Double_t Const[3];
    Double_t Mean[3];
    Double_t Sigma[3];
    Double_t DedxProb[3];
    Double_t MomProb[3];
    Double_t Mom[5];
    Double_t SumProb;
    Double_t prob;
    Double_t pPar[5] = {2.46E2,-9.706E-1,3.1679E-1,0,0};
    Double_t kPar[5] = {2.98385E2,-3.14546E1,6.72776E-1,3.66295E-1,-1.19806E-2};
    Double_t piPar[5] ={-1.3874E3,1.68399E3,-1.73976E2,6.25518,-7.50058E-2};
    Double_t MomNorm[3] = {15695.7,30905.2,105105.2}; 
    
    //Int_t index = 0;
    //static TMinuit *myMinuit = NULL;

    Const[0] = u[0];
    Const[1] = u[1];
    Const[2] = 1.0 - (Const[0]+Const[1]);
 
    SumProb = 0;
    for(i=0;i<NumPart;i++)
    {
        Mom[0] = 1;
        for(j=1;j<5;j++)
           Mom[j] = Mom[j-1] * val[i][0];
        Mean[0] = MinIon * RelRise((Float_t) Mom[1],(Float_t) kPMASS);
        Mean[1] = MinIon * RelRise((Float_t) Mom[1],(Float_t) kKMASS);
        Mean[2] = MinIon * RelRise((Float_t) Mom[1],(Float_t) kPIMASS);
        MomProb[0] = 0;
        for(j=0;j<5;j++)
            MomProb[0] += pPar[j]*Mom[j];
        MomProb[1] = 0;
        for(j=0;j<5;j++)
            MomProb[1] += kPar[j]*Mom[j];
        MomProb[2] = 0;
        for(j=0;j<5;j++)
            MomProb[2] += piPar[j]*Mom[j];
        prob = 0;
//        printf("p,dedx = %g,%gn", Mom[1],val[i][1]);
        for(j=0;j<3;j++)
        {
            Sigma[j] = Mean[j] * (Resolution + 0.001 * Mom[1]);
            DedxProb[j] = Norm/Sigma[j] * 
                        exp(-(val[i][1]-Mean[j])*(val[i][1]-Mean[j])/(2.0*Sigma[j]*Sigma[j]));
            
            MomProb[j] /= MomNorm[j];
//            printf("MomProb, DedxProb %d: %g - %g n",j, MomProb[j],DedxProb[j]);
            prob += MomProb[j] * DedxProb[j] * Const[j];
        }
        prob /= (Const[0]+Const[1]+Const[2]);
        if(prob > 0)
            SumProb  -= log(prob);  
    }
    f = SumProb;

}

 void T49Dedx::FitPRatio(Double_t Pmin, Double_t Pmax)
{
  //char CutStr[256];
  //char DrawStr[256];
  //char hname[256];
  //Double_t p;
  //Double_t step = 2;
  //Double_t ratios[5];

}

 Float_t T49Dedx::GetMultCorr(Int_t Mult, Float_t x)
{

    Float_t slope = -0.032 * (1 - 0.7 * fabs(x)/350);

    return 1 - ((slope * (Mult-654))/420);
}

 void T49Dedx::GetDedxPred(T49Particle *track, Float_t mass, Float_t *pos, Float_t *sigma)
{
    Float_t res;

    if(track->GetCharge() > 0)
    {
        *pos = GetRelRise(track->GetP(),mass)*fMinIonPos *(1 + 0.0008*track->GetP());
    }
    else
    {
        *pos = GetRelRise(track->GetP(),mass)*fMinIonNeg *(1 + 0.0005*track->GetP());
    }
    res    = pow(track->GetNDedxPoint(),0.43);
    *sigma = *pos * (fRes0/res + fRes1 * track->GetP());
}

 void T49Dedx::ReadCalib(Int_t RunNumber)
{
    
     Char_t FileName[512];

     sprintf(FileName,"%s/%d_dedx.root",fCalibHome,RunNumber);
     
     printf("Read calibration for run %d from %sn",RunNumber,FileName);

     TFile *f = new TFile(FileName);
     if(!f)
     {
          fprintf(stderr,"Can't open file %sn",FileName);
     }
     else
     {
          TProfile *h1 = (TProfile *) f->Get("dedx_vs_num_p");
          TProfile *h2 = (TProfile *) f->Get("dedx_vs_num_n");
      
          h1->Fit("pol1","0Q");
          TF1 *f1 = h1->GetFunction("pol1");

          h2->Fit("pol1","0Q");
          TF1 *f2 = h2->GetFunction("pol1");

          fCalibPos0 = f1->GetParameter(0);
          fCalibPos1 = f1->GetParameter(1);
          fCalibNeg0 = f2->GetParameter(0);
          fCalibNeg1 = f2->GetParameter(1);
          printf("Found Calibration Parameters %g,%g and %g,%gn",fCalibPos0,fCalibPos1,fCalibNeg0,fCalibNeg1);
     }
     f->Close();
}

 void T49Dedx::DoEventCalib(T49EventRoot *event, Bool_t DoMultCorr)
{
    T49Particle *track;
    Float_t x;
    Float_t y;
    //Float_t z;
    Float_t dx;
    Float_t dy;
    Float_t dz;

    Float_t dE;
    //Int_t   i;
    Int_t   mult;
    //Int_t   max;
    Float_t   kZCenter = 550;

    static Int_t OldRun = -1;

//  Primitve rejection of MC runs

    if(event->GetNRun() < 100)
    {
        return;
    }

    if(event->GetNRun() != OldRun)
    {
        ReadCalib(event->GetNRun());
        OldRun = event->GetNRun();
    }
    Float_t corr_pos = 400/(fCalibPos0 + fCalibPos1 * event->GetNEvent());
    Float_t corr_neg = 400/(fCalibNeg0 + fCalibNeg1 * event->GetNEvent());

    mult  = event->GetNParticles();
    //T49ParticleIter *iter = event->GetParticleIter();
    if(DoMultCorr)
    {
      //while((track = (T49Particle *) iter->Next()))
      //{
        TClonesArray *ParticleList = event->GetPrimaryParticles();
        for (Int_t iTrack = 0; iTrack < ParticleList->GetEntries(); iTrack++) {
            track = (T49Particle *) ParticleList->UncheckedAt(iTrack);
            if(track->GetCharge() > 0)
                track->SetTmeanCharge(track->GetTmeanCharge()*corr_pos,0);
            else
                track->SetTmeanCharge(track->GetTmeanCharge()*corr_neg,0);

            dx = track->GetXLast()-track->GetXFirst();
            dy = track->GetYLast()-track->GetYFirst();
            dz = track->GetZLast()-track->GetZFirst();
            x  = track->GetXFirst() + dx/dz*(kZCenter - track->GetZFirst());
            y  = track->GetYFirst() + dy/dz*(kZCenter - track->GetZFirst());
            dE = track->GetTmeanCharge() * GetMultCorr(mult,x);
            track->SetTmeanCharge(dE,0);
        }
    }
    else
    {
      //while((track = (T49Particle *) iter->Next()))
      //{
        TClonesArray *ParticleList = event->GetPrimaryParticles();
        for (Int_t iTrack = 0; iTrack < ParticleList->GetEntries(); iTrack++) {
            track = (T49Particle *) ParticleList->UncheckedAt(iTrack);
            if(track->GetCharge() > 0)
                track->SetTmeanCharge(track->GetTmeanCharge()*corr_pos,0);
            else
                track->SetTmeanCharge(track->GetTmeanCharge()*corr_neg,0);
        }
    }
}

 void T49Dedx::EventFit(Int_t Num)
{
    //Int_t     j;
    Int_t     i;
    Int_t     ErrFlag;
    //Float_t  *row;
    Float_t  row[7] = { 0 };
    Double_t arglist[10];
    Double_t we, al, bl;
    Double_t pNum = 0;
    Double_t kNum = 0;
    Double_t par;
    const Double_t Pmin = 5;
    const Double_t Pmax = 25;
    static Int_t index = 0;
    //T49Particle  *track;

    //Double_t kPar[5] = {1.0443E2,2.0724E1,-4.067,-5.349E-1,-1.3876E-2};
    //Double_t pPar[5] = {2.46E2,-9.706E-1,3.1679E-1,0,0};
    //Double_t piPar[5] ={-1.3874E3,1.68399E3,-1.73976E2,6.25518,-7.50058E-2};
    //Double_t pNorm  = 15695.7; 
    //Double_t kNorm  = 30905.2; 
    //Double_t piNorm = 105105.2; 
    Double_t k2pi;
    Double_t p2pi;
    
    static TMinuit *myMinuit = NULL;

    NumPart = 0;
    for(i=0;i<Num;i++) 
    {
        if(row[0] > Pmin && 
           row[0] < Pmax &&
           row[4] > 300 && 
           row[4] < 500 && 
           row[6] > MINDEDXPOINTS)
        {
            val[NumPart][0] = row[0];
            val[NumPart][1] = row[4];
            NumPart++;
        }
    }
    index += Num;    
    printf("Number of particles is %dn",NumPart);

////// Check if Minuit is initialized /////////////////////////////


   printf("Create myMinuitn");
   if(myMinuit == NULL)
   {
       myMinuit = new TMinuit();
       printf("Init Minuitn");
       myMinuit->mninit(5,6,7);
       printf("Call SetFCNn");
       myMinuit->SetFCN(TestFunc);
    }

   printf("Start event by event dedx fitn");
   myMinuit->mncler();
   myMinuit->mnparm(1,"rypr",0.1,0.01,0.0,1.0,ErrFlag);
   myMinuit->mnparm(2,"ryka",0.1,0.01,0.0,1.0,ErrFlag);
   arglist[0] = 0.5;

   myMinuit->mnexcm("SET ERR",arglist,1,ErrFlag);

   arglist[0] = 2;
   myMinuit->mnexcm("SET STR",arglist,1,ErrFlag);

   arglist[0] = 1;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   arglist[0] = 0;
   myMinuit->mnexcm("MINI",arglist,0,ErrFlag);

   arglist[0] = 3;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   arglist[0] = 6;
   myMinuit->mnexcm("CALL",arglist,1,ErrFlag);

   TString str;
   myMinuit->mnpout(1,str,par,we,al,bl,ErrFlag);
   //myMinuit->mnpout(1,"rypr",par,we,al,bl,ErrFlag);
   p2pi = par;
   myMinuit->mnpout(2,str,par,we,al,bl,ErrFlag);
   //myMinuit->mnpout(2,"ryka",par,we,al,bl,ErrFlag);
   k2pi = par;

   myMinuit->mnexcm("CLEAR",arglist,0,ErrFlag);
 
   printf("Finished event by event dedx fit, found %g %gn",pNum,kNum);

}

Float_t RelRise(Float_t mom, Float_t mass)
{
	Float_t epsilon; 		/**..dE/dx..**/ 
	Float_t zeta, K, Q;             /**..Target density,atomic shell, **/
        Float_t delta = 0;              /**.. charge,density corr...**/
	Float_t gamma, beta;
	Float_t p1, p2;		/**..Pressure in TPC in bar..**/
	Float_t a, m, X1, X0, XA;	/**..XA determines plateau height..**/
	Float_t epsmin, logbg, X;
        //Float_t betagamma, 

        Float_t corr;
        Int_t   sign;

        if(mom < 0)
            sign = -1;
        else
            sign = 1;

        mom = fabs(mom);

	m = 3.0;	/**..for density correction..**/
	Q = 1.0;
	p1 = 4.0;
	p2 = 1.0;
	epsmin = 30.0;

/***............Constants measured for Ar(88%)CH4(9.4%)iso-C4H10(2.6%) with JADE...........***/

	zeta = 0.42;
	K = 13.85;
	XA = 2.28;
	a = 0.26;


/***............Constants measured for Ar(88%)CH4(9.4%)iso-C4H10(2.6%) with OPAL...........

	zeta = 0.5;
	K = 11.3;
	XA = 2.28;
	a = 0.19;
		...***/

/***............Pressure dependencies...........***/

	zeta = zeta * (p2/p1);
	K = K + log(p2/p1);
	XA = XA - 0.5*log10(p2/p1);

	X0 = 1.771087313;
	X1 = 4.200915361;





/***............Defines unnormalized dE/dx as function of beta*gamma or mom/mass...........***/

		X = mom/mass;
		beta = sqrt((X*X)/(1+X*X));
		gamma = sqrt(1 + X*X);

		if(beta*gamma > 0.0)
		{
		logbg = log10(beta*gamma);
		if(logbg < X0)
		{
			delta = 0;
		}
		else if(X0 < logbg && logbg < X1)
		{
			delta = 2.0*log(10.0)*(logbg-XA) + a*(pow(X1-logbg,m));
		}
		else if(X1 < logbg)
		{
			delta = 2.0*log(10.0) * (logbg-XA);
		}

		epsilon = zeta * (Q*Q/(beta*beta))
			       * ( K + log(Q*Q)
			       + log(gamma*gamma) 
			       - beta*beta - delta);
		}
		else
			epsilon = 0.0;
		
                if(mom > 5 && mom < 40)
                {
                    if(sign < 0)
                    {
                       corr = 9.20540e-01 + 
                           2.19456e-02 * mom +
                          -2.31632e-03 * mom * mom +
                           1.12909e-04 * mom * mom * mom +
                          -2.53639e-06 * mom * mom * mom * mom + 
                           2.11768e-08 * mom * mom * mom * mom * mom ;
                    }
                    else
                    {
                       corr = 9.54256e-01  + 
                           7.81015e-03 * mom +
                          -3.57403e-04 * mom * mom +
                           4.44680e-06 * mom * mom * mom;
                    }
               }
               else
                   corr = 1;

                

		return(epsilon * DEDXNORM / corr);
}



ROOT page - Class index - Class Hierarchy - Top of the page

This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.