+++ /dev/null
-//////////////////////////////////////////////////////////////////////////
-// //
-// AliFast TrackMaker class. //
-// //
-// //
-//////////////////////////////////////////////////////////////////////////
-// ---------------------------------------------------------------------//
-// //
-// origin: "res.f" fortran by Karel Safarik which was used to //
-// calculate the track resolution for TP. //
-// Different detectors and material can be selected. //
-// The basic routines compute information and error matrices //
-// used for the calculation of momentum resolution. //
-// see references: ASK KAREL?? //
-// //
-// C++ in AliFast framework: Elzbieta Richter-Was and Yiota Foka //
-// following general structure od Makers in //
-// ATLFast by R. Brun. //
-// //
-// purpose: provide a Maker which by using general basic routines of //
-// "res.f" computes the necessary elements of covariance matrix//
-// for the calculation of Track Resolution. //
-// Those elements are the product of the TrackResolMaker and //
-// are hold in TrackResol class. They are expected to be used //
-// together with additional information for the calculation of //
-// the smeared momenta. //
-// Additional information necessary for this calculation //
-// will be provided via classes or functions specific to the //
-// specific study and/or detectors. //
-// One can select the detector and/or material for a specific //
-// study. //
-// //
-// starting point: res.f will be initialy partially contained in //
-// AliFTrackResolMaker and in AliFDet //
-// It will be reorganised further in the framework of //
-// AliFast according to the needs. //
-// Names of variables are kept as in fortran code. //
-// //
-//////////////////////////////////////////////////////////////////////////
-
-
-#ifdef WIN32
-// there is a bug in the Microsoft VisualC++ compiler
-// this class must be compiled with optimization off on Windows
-# pragma optimize( "", off )
-#endif
-
-#include <TFile.h>
-#include <TH1.h>
-#include <TMath.h>
-#include <TParticle.h>
-#include <TRandom.h>
-
-#include "AliFDet.h"
-#include "AliFTrack.h"
-#include "AliFTrackMaker.h"
-#include "AliFast.h"
-#include "AliMC.h"
-
-static const Double_t kPi = TMath::Pi();
-static const Double_t k2Pi = 2*kPi;
-static const Double_t kPiHalf = kPi/2.;
-extern AliFast * gAliFast;
-ClassImp(AliFTrackMaker)
-
-//_____________________________________________________________________________
-AliFTrackMaker::AliFTrackMaker()
-{
- //
- // Default constructor
- //
- fNTracks = 0;
- fResID1Test = 0;
- fResID2Test = 0;
- fResID3Test = 0;
- fResID4Test = 0;
- fResID5Test = 0;
-}
-
-//_____________________________________________________________________________
-AliFTrackMaker::AliFTrackMaker(const char *name, const char *title)
- :AliFMaker(name,title)
-{
- //
- // Standard Setters for tracks
- //
- fFruits = new TClonesArray("AliFTrack",100, kFALSE);
- fBranchName = "Tracks";
- fNTracks = 0;
-// Please, how to do this optionally ??!!!
- Save();
-}
-
-//_______________________________________________________________________
-AliFTrackMaker::AliFTrackMaker(const AliFTrackMaker& aftmk):
- AliFMaker(aftmk)
-{
- //
- // Copy constructor for AliRun
- //
- aftmk.Copy(*this);
-}
-
-//_____________________________________________________________________________
-AliFTrackMaker::~AliFTrackMaker()
-{
- //
- // Dummy constructor
- //
-}
-
-//_____________________________________________________________________________
-AliFTrack *AliFTrackMaker::AddTrack(Int_t code, Double_t charge,
- Double_t pT, Double_t eta,Double_t phi,
- Double_t v11, Double_t v22, Double_t v33,
- Double_t v12, Double_t v13, Double_t v23, Int_t iFlag)
-{
-// Add a new track to the list of tracks
-
- //Note the use of the "new with placement" to create a new track object.
- //This complex "new" works in the following way:
- // tracks[i] is the value of the pointer for track number i in the TClonesArray
- // if it is zero, then a new track must be generated. This typically
- // will happen only at the first events
- // If it is not zero, then the already existing object is overwritten
- // by the new track parameters.
- // This technique should save a huge amount of time otherwise spent
- // in the operators new and delete.
-
- TClonesArray &tracks = *(TClonesArray*)fFruits;
- return new(tracks[fNTracks++]) AliFTrack(code,charge,pT,eta,phi,
- v11, v22, v33, v12, v13, v23, iFlag);
-}
-
-//_____________________________________________________________________________
-void AliFTrackMaker::Clear(Option_t *option)
-{
- //Reset Track Maker
-
- fNTracks = 0;
- AliFMaker::Clear(option);
-}
-
-//_____________________________________________________________________________
-void AliFTrackMaker::Draw(Option_t *)
-{
-// Dummy Draw
-
-}
-
-//_____________________________________________________________________________
-void AliFTrackMaker::Init()
-{
- //Create control histograms
- if(gAliFast->TestTrack() == 0){
-
- fResID11 = new TH1D("ResID11","Elec: delta(1/pTotal)*pTotal",1000,-0.5,0.5);
- fResID12 = new TH1D("ResID12","Elec: delta(lambda)/lambda",1000,-0.01,0.01);
- fResID13 = new TH1D("ResID13","Elec: delta(phi)/phi",1000,-0.01,0.01);
-
- fResID21 = new TH1D("ResID21","Pion: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
- fResID22 = new TH1D("ResID22","Pion: delta(lambda)/lambda",1000,-1.0,1.0);
- fResID23 = new TH1D("ResID23","Pion: delta(phi)/phi",1000,-1.0,1.0);
-
- fResID31 = new TH1D("ResID31","Kaon: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
- fResID32 = new TH1D("ResID32","Kaon: delta(lambda)/lambda",1000,-1.0,1.0);
- fResID33 = new TH1D("ResID33","Kaon: delta(phi)/phi",1000,-1.0,1.0);
-
- fResID41 = new TH1D("ResID41","Proton: delta(1/pTotal)*pTotal",1000,-1.0,1.0);
- fResID42 = new TH1D("ResID42","Proton: delta(lambda)/lambda",1000,-1.0,1.0);
- fResID43 = new TH1D("ResID43","Proton: delta(phi)/phi",1000,-1.0,1.0);
-
- }
- //Create test histograms for TestJob only
- if(gAliFast->TestTrack() == 1){
- fResID1Test = new TH1D("ResID1Test","histogram21 from res.f",1000,0.075,10.075);
- fResID2Test = new TH1D("ResID2Test","histogram21 from res.f",1000,0.075,10.075);
- fResID3Test = new TH1D("ResID3Test","histogram21 from res.f",1000,0.075,10.075);
- fResID4Test = new TH1D("ResID4Test","histogram21 from res.f",1000,0.075,10.075);
- fResID5Test = new TH1D("ResID5Test","histogram21 from res.f",1000,0.075,10.075);
- }
-
- //Set particle masses
- SetPionMass();
- SetKaonMass();
- SetElectronMass();
- SetProtonMass();
-
- //Switch on/off tracks reconstruction
- SetRecTrack();
-
-}
-
-void AliFTrackMaker::Make()
-{
- //
- // Calculate track and its resolution
- //
- Double_t v11, v22, v33, v12, v13, v23;
- Int_t iFlag;
-
- fNTracks = 0;
-
- // Check if it is a TestJob
- if(gAliFast->TestTrack() == 1){
- // Run test job
- MakeTest(10);
- }else{
- // Run production job
- // Get pointers to Particles arrays and TClonesArray
-
- Int_t idPart, idTrack;
- Double_t charge, pT, eta, phi;
- TParticle *part;
- Int_t nparticles = gAlice->GetMCApp()->GetNtrack();
- printf("%10s%10d\n","nparticles",nparticles);
- for(Int_t ind=0;ind<nparticles;ind++) {
- part = gAlice->GetMCApp()->Particle(ind);
- idPart = part->GetPdgCode();
- charge = part->GetPDG()->Charge();
- pT = part->Pt();
- eta = part->Eta();
- phi = part->Phi();
- printf("%10s%10d%20.5e%20.5e%20.5e%20.5e\n","Particle",idPart,charge,pT,eta,phi);
- // Check convention for tracks reconstruction
- idTrack = 0;
- if(TMath::Abs(idPart) == 11) idTrack = 1;
- if(TMath::Abs(idPart) == 111 || TMath::Abs(idPart) == 211) idTrack = 2;
- if(TMath::Abs(idPart) == 311 || TMath::Abs(idPart) == 321) idTrack = 3;
- if(TMath::Abs(idPart) == 2212) idTrack = 4;
-
- if(idTrack > 0 && fRecTrack > 0){
- // Check if track should be reconstructed
- if((fRecTrack == 1 && idTrack == 1) ||
- (fRecTrack == 2 && idTrack == 2) ||
- (fRecTrack == 3 && idTrack == 3) ||
- (fRecTrack == 4 && idTrack == 4) ||
- fRecTrack == 100 ) {
- // Tracks are reconstructed
- ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
-
- // Calculate and smear track parameters
- Double_t lambda, cosLambda, pTotal,pInverse;
- Double_t pInverseSmea, lambdaSmea, phiSmea;
- Double_t a1, a2, a3, b2, b3, c3;
- Double_t rn1, rn2, rn3;
-
- lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
- cosLambda = TMath::Cos(lambda);
- pTotal = pT/cosLambda;
- pInverse = 1.0/pTotal;
-
- a1 = TMath::Sqrt(v11);
- if(a1 == 0.){
- a2 = 0;
- a3 = 0;
- }else{
- a2 = v12/a1;
- a3 = v13/a1;
- }
- b2 = TMath::Sqrt(v22-a2*a2);
- if(b2 == 0.){
- b3 = 0;
- }else{
- b3 = (v23 - a2*a3)/b2;
- }
- c3 = TMath::Sqrt(v33 - a3*a3 -b3*b3);
- rn1 = gRandom->Gaus(0,1);
- rn2 = gRandom->Gaus(0,1);
- rn3 = gRandom->Gaus(0,1);
-
- pInverseSmea = pInverse + a1*rn1;
- lambdaSmea = lambda + a2*rn1 + b2*rn2;
- phiSmea = phi + a3*rn1 + b3*rn2 + c3*rn3;
-
- // Fill control histograms
- if(idTrack == 1){
- fResID11->Fill((pInverseSmea-pInverse)/pInverse);
- fResID12->Fill((lambdaSmea-lambda)/lambda);
- fResID13->Fill((phiSmea-phi)/phi);
- }
- else if(idTrack == 2){
- fResID21->Fill((pInverseSmea-pInverse)/pInverse);
- fResID22->Fill((lambdaSmea-lambda)/lambda);
- fResID23->Fill((phiSmea-phi)/phi);
- }
- else if(idTrack == 3){
- fResID31->Fill((pInverseSmea-pInverse)/pInverse);
- fResID32->Fill((lambdaSmea-lambda)/lambda);
- fResID33->Fill((phiSmea-phi)/phi);
- }
- else if(idTrack == 4){
- fResID41->Fill((pInverseSmea-pInverse)/pInverse);
- fResID42->Fill((lambdaSmea-lambda)/lambda);
- fResID43->Fill((phiSmea-phi)/phi);
- }
- }else{
- // Tracks are not reconstructed
- v11=0.;
- v12=0.;
- v13=0.;
- v22=0.;
- v23=0.;
- v33=0.;
- iFlag=0;
- }
- // Store resolution variables to AliFTrack ClonesArray
- AddTrack(idTrack, charge, pT, eta, phi, v11, v22, v33, v12, v13, v23, iFlag);
- printf("%10s%10d%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%20.5e%10d\n",
- "Track",idTrack,charge,pT,eta,phi,v11, v22, v33, v12, v13, v23, iFlag);
- }
-
- }
- }
-}
-
-//_______________________________________________________________________
-void AliFTrackMaker::Copy(TObject &) const
-{
- Fatal("Copy","Not implemented!\n");
-}
-
-//_____________________________________________________________________________
-void AliFTrackMaker::Finish()
-{
- // For TestJob only
- if(gAliFast->TestTrack() == 1){
- /*
- // Draw test histograms
- TCanvas *c1 = new TCanvas("c1"," ",200,10,600,480);
- c1->Divide(2,3);
- c1->cd(1); fResID1Test->Draw();
- c1->cd(2); fResID2Test->Draw();
- c1->cd(3); fResID3Test->Draw();
- c1->cd(4); fResID4Test->Draw();
- c1->cd(5); fResID5Test->Draw();
- c1->Update();
- // Store TestRes.eps file
- c1->Print("TestRes.eps");
- */
- // Store histograms on file
- TFile f2("TestRes.root","RECREATE","Test Res.f");
- fResID1Test->Write();
- fResID2Test->Write();
- fResID3Test->Write();
- fResID4Test->Write();
- fResID5Test->Write();
- f2.Close();
- }
-}
-//_____________________________________________________________________________
-void AliFTrackMaker::ErrorMatrix(Int_t idTrack, Double_t pT, Double_t eta,
- Double_t &v11, Double_t &v22, Double_t &v33, Double_t &v12, Double_t &v13, Double_t &v23,
- Int_t &iFlag)
-{
- ///////////////////////////////////////////////
- //idTrack track type input //
- //pT transverse mom input //
- //lambda deep angle input //
- //v11,v22,v23 error matrix output //
- //v12,v13,v23 output //
- //iFlag output //
- ///////////////////////////////////////////////
-
- AliFDet *detector = gAliFast->Detector();
- Int_t nDet = detector->NDet();
- Int_t nDetActive = detector->NDetActive();
- Int_t nTwice = nDetActive + nDetActive;
-
- Double_t rTrack, rTrackInverse, pTotal, pInverse, diffPInverse;
- Double_t safety;
- Double_t cosLambda, tanLambda, diffLambda;
- Double_t rDet;
-
- Double_t hh0[kNMaxDet2][kNMaxDet2], hhi0[kNMaxDet2][kNMaxDet2];
- Double_t hh1[kNMaxDet2][kNMaxDet2], hhi1[kNMaxDet2][kNMaxDet2];
- Double_t dhhiOverPInverse[kNMaxDet2][kNMaxDet2];
- Double_t dhhiOverLambda[kNMaxDet2][kNMaxDet2];
- Double_t a1[kNMaxDet2][kNMaxDet2], a2[kNMaxDet2][kNMaxDet2];
- Double_t a0PInverse[kNMaxDet2];
- Double_t a0Lambda[kNMaxDet2];
- Double_t a0Phi[kNMaxDet2];
-
- Double_t vF11, vF12, vF13, vF22, vF23, vF33, d1, d2, d3, det;
- Int_t idet, icyl, im, in;
- Double_t phiHalf;
- Double_t lambda;
-
- lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
- rTrack = detector->ConstMag()*pT;
- safety = 10.0;
- if(2.0*rTrack < (detector->RDet(nDet) + safety)){
- iFlag = 0;
- v11 = 0;
- v22 = 0;
- v33 = 0;
- v12 = 0;
- v13 = 0;
- v23 = 0;
- return;
- }
- iFlag = 1;
- cosLambda = TMath::Cos(lambda);
- pTotal = pT/cosLambda;
- pInverse = 1.0/pTotal;
- diffPInverse = pInverse*1.0e-5;
- diffLambda = 1.0e-4;
-
- // Compute likelihood and derivatives
-
- LogLikelyhood(idTrack, pInverse, lambda);
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<nTwice+1; im++){
- hh0[icyl][im] = HH(icyl,im);
- hhi0[icyl][im] = HHI(icyl,im);
- }
- }
- LogLikelyhood(idTrack, pInverse+diffPInverse,lambda);
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<nTwice+1; im++){
- hh1[icyl][im] = HH(icyl,im);
- hhi1[icyl][im] = HHI(icyl,im);
- }
- }
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<icyl+1; im++){
- dhhiOverPInverse[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffPInverse;
- }
- }
- LogLikelyhood(idTrack, pInverse, lambda+diffLambda);
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<nTwice+1; im++){
- hh1[icyl][im] = HH(icyl,im);
- hhi1[icyl][im] = HHI(icyl,im);
- }
- }
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<icyl+1; im++){
- dhhiOverLambda[icyl][im] = (hhi1[icyl][im]-hhi0[icyl][im])/diffLambda;
- }
- }
-
- // Compute additional derivatives
- rTrackInverse = 1.0/rTrack;
- tanLambda = TMath::Tan(lambda);
- icyl = 0;
- for(idet=1; idet<nDet+1;idet++){
- if(detector->IFlagDet(idet) > 0){
- icyl = icyl + 1;
- rDet = detector->RDet(idet);
- phiHalf = TMath::ASin(0.5*rDet*rTrackInverse);
- Double_t rHelp = rDet /
- (2.0 * TMath::Sqrt(1.0-(0.5 *rDet*rTrackInverse)*
- (0.5 *rDet*rTrackInverse)));
- a0PInverse[icyl] = - rDet* rHelp
- /(detector->ConstMag()*cosLambda);
- a0Lambda[icyl] = - rDet* rHelp
- * tanLambda * rTrackInverse;
- a0Phi[icyl] = rDet;
- a0PInverse[nDetActive+icyl] = 2.0 * tanLambda
- *rTrack*(rHelp-rTrack*phiHalf)
- /(detector->ConstMag()*cosLambda);
- a0Lambda[nDetActive+icyl] = 2.0 * ( rHelp*tanLambda*tanLambda
- + rTrack*phiHalf);
- a0Phi[nDetActive+icyl] = 0.0 ;
- }
- }
-
- // Compute information matrix
-
- vF11=0.0;
- vF12=0.0;
- vF13=0.0;
- vF22=0.0;
- vF23=0.0;
- vF33=0.0;
- for(icyl=1; icyl<nTwice+1; icyl++){
- d1=0.0;
- d2=0.0;
- d3=0.0;
- for(im=1; im < icyl+1; im++){
- d1 = d1 + hhi0[icyl][im]*a0PInverse[im];
- d2 = d2 + hhi0[icyl][im]*a0Lambda[im];
- d3 = d3 + hhi0[icyl][im]*a0Phi[im];
- }
- vF11 =vF11 + d1*d1;
- vF12 =vF12 + d1*d2;
- vF13 =vF13 + d1*d3;
- vF22 =vF22 + d2*d2;
- vF23 =vF23 + d2*d3;
- vF33 =vF33 + d3*d3;
- }
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<icyl+1; im++){
- a1[icyl][im] = 0;
- a2[icyl][im] = 0;
- for(in=im; in<icyl+1;in++){
- a1[icyl][im]=a1[icyl][im]+dhhiOverPInverse[icyl][in]*hh0[im][in];
- a2[icyl][im]=a2[icyl][im]+dhhiOverLambda[icyl][in]*hh0[im][in];
- }
- vF11=vF11+a1[icyl][im]*a1[icyl][im];
- vF12=vF12+a1[icyl][im]*a2[icyl][im];
- vF22=vF22+a2[icyl][im]*a2[icyl][im];
- }
- vF11=vF11+a1[icyl][icyl]*a1[icyl][icyl];
- vF12=vF12+a1[icyl][icyl]*a2[icyl][icyl];
- vF22=vF22+a2[icyl][icyl]*a2[icyl][icyl];
- }
-
- // Invert information matrix
-
- det=( vF11*vF22 - vF12*vF12 ) *vF33 + (vF12*vF23 - vF13*vF22)*vF13
- + (vF12*vF13 - vF11*vF23)*vF23;
-
- v11 = (vF22*vF33 - vF23*vF23)/det;
- v22 = (vF11*vF33 - vF13*vF13)/det;
- v33 = (vF11*vF22 - vF12*vF12)/det;
- v12 = (vF13*vF23 - vF12*vF33)/det;
- v13 = (vF12*vF23 - vF13*vF22)/det;
- v23 = (vF12*vF13 - vF11*vF23)/det;
-
- }
-//_____________________________________________________________________________//
-void AliFTrackMaker::LogLikelyhood(Int_t idTrack, Double_t pInverse,Double_t lambda)
-{
- ///////////////////////////////////////////////
- //hh ?? output //
- //hhi ?? output //
- //idTrack track type input //
- //pInverse inverse momentum input //
- //lambda polar angle of track input //
- ///////////////////////////////////////////////
-
-
- AliFDet *detector = gAliFast->Detector();
- Int_t nDet = detector->NDet();
- Int_t nDetActive = detector->NDetActive();
- Int_t nTwice = nDetActive + nDetActive;
-
- Double_t rDet, rDetSQ;
- Int_t idet, icyl, im, imc;
- Double_t cosLambda, tanLambda, pTotal, pT, rTrack, rTrackSQ;
- Double_t beta, overPBeta, rTrackInv, thickCorr, temp1, temp2;
- Double_t partMassSQ;
- Double_t aShelp[kNMaxDet2], dShelp[kNMaxDet2];
- Double_t projXVXT[kNMaxDet2],projYVXT[kNMaxDet2], projZVXT[kNMaxDet2];
- Double_t proj[kNMaxDet2][kNMaxDet2];
- Double_t erroScatt[kNMaxDet2], variance[kNMaxDet2][kNMaxDet2];
- Double_t erroSQ[kNMaxDet2];
- Double_t hh[kNMaxDet2][kNMaxDet2];
- Double_t hhi[kNMaxDet2][kNMaxDet2];
- Double_t errorVX, errorVY, errorVZ;
-
- cosLambda = TMath::Cos(lambda);
- tanLambda = TMath::Tan(lambda);
- pTotal = 1.0/pInverse;
- pT = pTotal * cosLambda;
- rTrack = detector->ConstMag() * pTotal * cosLambda;
- rTrackSQ = rTrack * rTrack;
- partMassSQ= ParticleMass(idTrack)*ParticleMass(idTrack);
- beta = pTotal / TMath::Sqrt(partMassSQ+pTotal*pTotal);
- overPBeta = 1./(pTotal*beta);
- rTrackInv = 1./rTrack;
- errorVX = detector->ErrorVertexX();
- errorVY = detector->ErrorVertexY();
- errorVZ = detector->ErrorVertexZ();
-
-
- erroScatt[0]=0.0;
- erroScatt[1]=0.0;
- for(idet=1; idet < nDet; idet++){
- thickCorr = detector->ThickDet(idet)/TMath::Sqrt(cosLambda*
- TMath::Sqrt(1.0-0.25*(detector->RDetSQ(idet)/rTrackSQ)));
- if(detector->IFlagGas(idet) == 0){
- thickCorr = thickCorr * (1.3266 + 0.076 * TMath::Log(thickCorr));}
- thickCorr = overPBeta * thickCorr;
- erroScatt[idet+1]=thickCorr*thickCorr;
- }
-
-
- icyl = 0;
- for(idet=1; idet<nDet+1; idet++){
- rDet = detector->RDet(idet);
- rDetSQ = rDet*rDet;
- dShelp[idet] = TMath::Sqrt(4.0*rTrackSQ-rDetSQ);
- aShelp[idet] = TMath::ASin(rDet/(rTrack+rTrack));
- if(detector->IFlagDet(idet) > 0) {
- icyl = icyl + 1;
- projXVXT[icyl] = rDet * rTrackInv;
- projXVXT[nDetActive+icyl] = -tanLambda;
- temp1 = (rTrackSQ + rTrackSQ - rDetSQ)/dShelp[idet];
- temp2 = rDet/dShelp[idet];
- projYVXT[icyl] = temp1*rTrackInv;
- projYVXT[nDetActive+icyl] = tanLambda * temp2;
- projZVXT[icyl] = 0.0;
- projZVXT[nDetActive+icyl] = 1.0;
- proj[icyl][1] = 0.0;
- proj[nDetActive+icyl][0] = 0.0;
- proj[nDetActive+icyl][nDet] = 0.0;
- proj[icyl][nDet] = 0.0;
- for(im=2; im<idet+1; im++){
- proj[icyl][im]= (( rDet
- *(rTrackSQ+rTrackSQ-detector->RDetSQ(im-1))
- - detector->RDet(im-1)*temp1*dShelp[im-1])
- /((rTrackSQ + rTrackSQ)*cosLambda));
- proj[nDetActive+icyl][im]= 0.5 * detector->RDet(im-1)
- * rTrackInv
- * tanLambda * (detector->RDet(im-1)
- - dShelp[im-1]*temp2)/cosLambda;
- proj[nDetActive+icyl][nDet+im]= (rTrack+rTrack) * (aShelp[idet] - aShelp[im-1])
- + ( rDet*dShelp[im-1]-detector->RDet(im-1)*dShelp[idet])
- * dShelp[im-1] * tanLambda * tanLambda
- / (dShelp[idet] * (rTrack+rTrack));
- proj[icyl][nDet+im]= tanLambda
- * (rDet*detector->RDet(im-1)*dShelp[im-1]
- / (rTrackSQ+rTrackSQ)
- - (rDetSQ + detector->RDetSQ(im-1)
- - rDetSQ * detector->RDetSQ(im-1)
- / (rTrackSQ+rTrackSQ))
- / dShelp[idet]);
- }
- for(im=idet+1; im < nDet+1; im++){
- proj[icyl][im] = 0.0;
- proj[nDetActive+icyl][im] = 0.0;
- proj[nDetActive+icyl][nDet+im] = 0.0;
- proj[icyl][nDet+im] = 0.0;
- }
- if(detector->IFlagDet(idet) == 1){
- erroSQ[icyl] = detector->ErrorRPhi(idet);
- erroSQ[nDetActive+icyl] = detector->ErrorZ(idet);
- }else{
- TPCResolution(pT, rDet, lambda);
- erroSQ[icyl] = SigmaRPhiSQ();
- erroSQ[nDetActive+icyl] = SigmaZSQ();
- }
- erroSQ[icyl] = erroSQ[icyl] + detector->ErrorR(idet)*temp2*temp2;
- erroSQ[nDetActive+icyl] = erroSQ[nDetActive+icyl]
- + detector->ErrorR(idet)*tanLambda*tanLambda;
- }
- }
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<icyl+1; im++){
- variance[icyl][im]=
- projXVXT[icyl]*projXVXT[im]*errorVX
- +projYVXT[icyl]*projYVXT[im]*errorVY
- +projZVXT[icyl]*projZVXT[im]*errorVZ;
- for(imc=1; imc<nDet+1; imc++){
- variance[icyl][im]=variance[icyl][im]
- +(proj[icyl][imc]*proj[im][imc]
- + proj[icyl][nDet+imc]*proj[im][nDet+imc])
- * erroScatt[imc];
- }
- }
- variance[icyl][icyl] = variance[icyl][icyl]+erroSQ[icyl];
- }
-
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=icyl; im<nTwice+1; im++){
- hh[icyl][im]=variance[im][icyl];
- for(imc=1; imc<icyl;imc++){
- hh[icyl][im]=hh[icyl][im]-hh[imc][icyl]*hh[imc][im];
- }
- if(im == icyl){
- hh[icyl][im] = TMath::Sqrt(hh[icyl][im]);
- } else {
- hh[icyl][im] = hh[icyl][im]/hh[icyl][icyl];
- }
- }
- }
-
- for(icyl=1; icyl<nTwice+1; icyl++){
- hhi[icyl][icyl] = 1.0 / hh[icyl][icyl];
- for(im=1; im<icyl; im++){
- hhi[icyl][im] = 0.0;
- for(imc=im; imc<icyl; imc++){
- hhi[icyl][im] = hhi[icyl][im]-hh[imc][icyl]*hhi[imc][im];
- }
- hhi[icyl][im] = hhi[icyl][im]*hhi[icyl][icyl];
- }
- }
-
- for(icyl=1; icyl<nTwice+1; icyl++){
- for(im=1; im<nTwice+1; im++){
- SetHH(icyl,im,hh[icyl][im]);
- SetHHI(icyl,im,hhi[icyl][im]);
- }
- }
-
-
-}
-//_____________________________________________________________________________
-// translation of routine tpc_resolution of res.f
-//_____________________________________________________________________________
-void AliFTrackMaker::TPCResolution(Double_t pTransv, Double_t radiPad, Double_t lambda)
-{
- ///////////////////////////////////////////////
- //sigmaRPhiSQ resolution in r-phi output //
- //sigmaZSQ resolution in z output //
- //pTransv transverse momentum input //
- //radiPad radius of pad row input //
- //lambda polar angle of track input //
- //
- //units: cm, GeV/c, radian //
- //parametrisation of TPC resolution //
- //version of 03.07.1995 //
- //source: Marek Kowalski, Karel Safarik //
- ///////////////////////////////////////////////
-
- Double_t aRCoeff=0.41818e-2;
- Double_t bRCoeff=0.17460e-4;
- Double_t cRCoeff=0.30993e-8;
- Double_t dRCoeff=0.41061e-6;
- Double_t aZCoeff=0.39610e-2;
- Double_t bZCoeff=0.22443e-4;
- Double_t cZCoeff=0.51504e-1;
-
- Double_t sigmaRPhiSQ;
- Double_t sigmaZSQ;
-
- sigmaRPhiSQ = aRCoeff - bRCoeff * radiPad * TMath::Tan(lambda)+
- (cRCoeff * (radiPad/pTransv) + dRCoeff) * radiPad/pTransv;
-
- sigmaZSQ = aZCoeff - bZCoeff * radiPad * TMath::Tan(lambda)+
- cZCoeff * TMath::Tan(lambda)*TMath::Tan(lambda);
-
- if(sigmaRPhiSQ < 1.0e-6 ) sigmaRPhiSQ = 1.0e-6;
- if(sigmaZSQ < 1.0e-6 ) sigmaZSQ = 1.0e-6;
-
- sigmaRPhiSQ = (TMath::Sqrt(sigmaRPhiSQ) + 0.005)
- * (TMath::Sqrt(sigmaRPhiSQ) + 0.005);
- sigmaZSQ = (TMath::Sqrt(sigmaZSQ) + 0.005)
- * (TMath::Sqrt(sigmaZSQ) + 0.005);
-
- SetSigmaRPhiSQ(sigmaRPhiSQ);
- SetSigmaZSQ(sigmaZSQ);
-
-
-}
-
-//-----------------------------------------------------------------------------
-Double_t AliFTrackMaker::ParticleMass(Int_t idTrack) const
-{
- //
- // returns the mass given particle ID
- //
- Double_t mass = 0.0;
-
- if(idTrack == 2){ mass = fPionMass;}
- else if(idTrack == 3){ mass = fKaonMass;}
- else if(idTrack == 1) {mass = fElectronMass;}
- else if(idTrack == 4) {mass = fProtonMass;}
-
- return mass;
-
-}
-
-//____________________________________________________________________________
-Double_t AliFTrackMaker::Rapidity(Double_t pt, Double_t pz)
-{
- //
- // returns the rapidity given particle pT, pz
- // Compute rapidity
-
- Double_t etalog = TMath::Log((TMath::Sqrt(pt*pt + pz*pz) + TMath::Abs(pz))/pt);
- if (pz < 0 ) return -TMath::Abs(etalog);
- else return TMath::Abs(etalog);
-}
-
-//_____________________________________________________________________________
-// returns the phi angle given particle px, py
-//-----------------------------------------------------------------------------
-Double_t AliFTrackMaker::Angle(Double_t x, Double_t y)
-{
-// Compute phi angle of particle
-// ... this is a copy of function ULANGL
-// .. sign(a,b) = -abs(a) if b < 0
-// = abs(a) if b >= 0
-
- Double_t angle = 0;
- Double_t r = TMath::Sqrt(x*x + y*y);
- if (r < 1e-20) return angle;
- if (TMath::Abs(x)/r < 0.8) {
- angle = TMath::Sign((Double_t)TMath::Abs(TMath::ACos(x/r)), y);
- } else {
- angle = TMath::ASin(y/r);
- if (x < 0 ) {
- if(angle >= 0) angle = kPi - angle;
- else angle = -kPi - angle;
- }
- }
- return angle;
-}
-//_____________________________________________________________________________
-Int_t AliFTrackMaker::Charge(Int_t kf)
-{
-//...this is copy of function LUCHGE
-//...Purpose: to give three times the charge for a particle/parton.
-
- static Int_t kchg[500] = { -1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,-3,0,
- 0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,3,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 2,-1,2,-1,2,3,0,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,
- 0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,
- 0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,
- 0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,
- 0,3,0,0,0,0,0,0,0,0,-3,0,0,0,0,0,0,0,0,3,0,-3,0,3,-3,0,0,0,3,6,0,
- 3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,-3,0,3,6,-3,0,3,-3,0,-3,0,3,6,
- 0,3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
-
-// extern integer kfcomp_(integer *);
- Int_t ipower;
- Int_t ret = 0;
- Int_t kfa = TMath::Abs(kf);
- Int_t kc = Compress(kfa);
-
-//...Initial values. Simple case of direct readout.
- if (kc == 0) {
- } else if (kfa <= 100 || kc <= 80 || kc > 100) {
- ret = kchg[kc-1];
-
-// ...Construction from quark content for heavy meson, diquark, baryon.
- } else if (kfa/1000 % 10 == 0) {
- ipower = kfa/100 % 10;
- ret = (kchg[kfa / 100 % 10 - 1] - kchg[kfa/10 % 10 - 1])*Int_t(TMath::Power(-1, ipower));
- } else if (kfa / 10 % 10 == 0) {
- ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1];
- } else {
- ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1] + kchg[kfa/10 % 10 - 1];
- }
-
-// ...Add on correct sign.
- if (kf > 0) return ret;
- else return -ret;
-}
-//_____________________________________________________________________________
-Int_t AliFTrackMaker::Compress(Int_t kf)
-{
-//...this is copy of function LUCOMP
-//...Purpose: to compress the standard KF codes for use in mass and decay
-//...arrays; also to check whether a given code actually is defined.
-// from BLOCK LUDATA
- static Int_t kchg[500] = { 1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0,
- 0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,1,
- 1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,
- 1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
- 0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,
- 1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,
- 0,0,1,0,1,1,0,0,0,0,0,0,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,0,0,0,0,1,1,
- 1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,
- 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
- static Int_t kftab[25] = { 211,111,221,311,321,130,310,213,113,223,313,
- 323,2112,2212,210,2110,2210,110,220,330,440,30443,30553,0,0 };
- static Int_t kctab[25] = { 101,111,112,102,103,221,222,121,131,132,122,
- 123,332,333,281,282,283,284,285,286,287,231,235,0,0 };
-
- Int_t ret = 0;
- Int_t kfla, kflb, kflc, kflr, kfls, kfa, ikf;
-
- kfa = TMath::Abs(kf);
-//...Simple cases: direct translation or table.
- if (kfa == 0 || kfa >= 100000) {
- return ret;
- } else if (kfa <= 100) {
- ret = kfa;
- if (kf < 0 && kchg[kfa - 1] == 0) ret = 0;
- return ret;
- } else {
- for (ikf = 1; ikf <= 23; ++ikf) {
- if (kfa == kftab[ikf-1]) {
- ret = kctab[ikf-1];
- if (kf < 0 && kchg[ret-1] == 0) ret = 0;
- return ret;
- }
- }
- }
-// ...Subdivide KF code into constituent pieces.
- kfla = kfa / 1000%10;
- kflb = kfa / 100%10;
- kflc = kfa / 10%10;
- kfls = kfa%10;
- kflr = kfa / 10000%10;
-// ...Mesons.
- if (kfa - kflr*10000 < 1000) {
- if (kflb == 0 || kflb == 9 || kflc == 0 || kflc == 9) {
- } else if (kflb < kflc) {
- } else if (kf < 0 && kflb == kflc) {
- } else if (kflb == kflc) {
- if (kflr == 0 && kfls == 1) { ret = kflb + 110;
- } else if (kflr == 0 && kfls == 3) { ret = kflb + 130;
- } else if (kflr == 1 && kfls == 3) { ret = kflb + 150;
- } else if (kflr == 1 && kfls == 1) { ret = kflb + 170;
- } else if (kflr == 2 && kfls == 3) { ret = kflb + 190;
- } else if (kflr == 0 && kfls == 5) { ret = kflb + 210;
- }
- } else if (kflb <= 5) {
- if (kflr == 0 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 100 + kflc;
- } else if (kflr == 0 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 120 + kflc;
- } else if (kflr == 1 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 140 + kflc;
- } else if (kflr == 1 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 160 + kflc;
- } else if (kflr == 2 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 180 + kflc;
- } else if (kflr == 0 && kfls == 5) { ret = (kflb-1)*(kflb-2)/2 + 200 + kflc;
- }
- } else if (kfls == 1 && kflr <= 1 || kfls == 3 && kflr <= 2 || kfls == 5 && kflr == 0) {
- ret = kflb + 80;
- }
-// ...Diquarks.
- } else if ((kflr == 0 || kflr == 1) && kflc == 0) {
- if (kfls != 1 && kfls != 3) {
- } else if (kfla == 9 || kflb == 0 || kflb == 9) {
- } else if (kfla < kflb) {
- } else if (kfls == 1 && kfla == kflb) {
- } else { ret = 90;
- }
-// ...Spin 1/2 baryons.
- } else if (kflr == 0 && kfls == 2) {
- if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
- } else if (kfla <= kflc || kfla < kflb) {
- } else if (kfla >= 6 || kflb >= 4 || kflc >= 4) {
- ret = kfla + 80;
- } else if (kflb < kflc) {
- ret = (kfla+1)*kfla*(kfla-1)/6 + 300 + kflc*(kflc-1)/2 + kflb;
- } else {
- ret = (kfla+1)*kfla*(kfla-1)/6 + 330 + kflb*(kflb-1)/2 + kflc;
- }
-// ...Spin 3/2 baryons.
- } else if (kflr == 0 && kfls == 4) {
- if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
- } else if (kfla < kflb || kflb < kflc) {
- } else if (kfla >= 6 || kflb >= 4) {
- ret = kfla + 80;
- } else {
- ret = (kfla+1)*kfla*(kfla-1) / 6 + 360 + kflb*(kflb -1) / 2 + kflc;
- }
- }
- return ret;
-}
-
-//_____________________________________________________________________________
-void AliFTrackMaker::MakeTest(Int_t n)
-{
- //
- // TEST JOB: Calculate tracks resolution
- //
- Double_t v11, v22, v33, v12, v13, v23;
- Int_t iFlag;
- Int_t idTrack;
- Double_t pTStart, pT, eta;
-
- Double_t sumDPop,sumDDip,sumDPhi;
- Double_t isum,fm;
- Double_t pTotal,partMassSQ,beta,lambda;
- Double_t dPop,dLop,dDip,dPhi,rho12,rho13,rho23;
- Double_t dPPStrag,dPPTot=0;
- // Double_t resol1[1001][11],resol2[10001][11],resol3[1001][11],
- // resol4[1001][11],resol5[10001][11]
- Double_t store1[1001],store2[10001],store3[1001],
- store4[1001],store5[10001];
-
-
- idTrack = 2;
- pTStart = 0.07;
- for(Int_t istep=1; istep<n; istep++){
- if(istep < 100 && istep > 20) istep = istep -1 + 5;
- if(istep < 500 && istep > 100) istep = istep -1 + 25;
- if(istep <1000 && istep > 500) istep = istep -1 + 100;
- pT = pTStart + 0.01*istep;
- eta = - 0.044;
- sumDPop = 0;
- sumDDip = 0;
- sumDPhi = 0;
- isum = 0;
- for(Int_t in=1; in<11; in++){
- eta = eta + 0.088;
- lambda = kPiHalf -2.0*TMath::ATan(TMath::Exp(-eta));
- pTotal = pT / TMath::Cos(lambda);
- if(idTrack == 1){
- dPPStrag = 0.055 /pT;}
- else{
- partMassSQ = ParticleMass(idTrack)*ParticleMass(idTrack);
- beta = pTotal/ TMath::Sqrt(pTotal*pTotal + partMassSQ);
- dPPStrag = 0.04/(pT*TMath::Power(beta,2.6666666667));
- }
- ErrorMatrix(idTrack,pT,eta, v11, v22, v33, v12, v13, v23, iFlag);
- if(iFlag == 1){
- dLop = TMath::Sqrt(v11);
- dDip = TMath::Sqrt(v22);
- dPhi = TMath::Sqrt(v33);
- rho12 = v12/(dLop*dDip);
- rho13 = v13/(dLop*dPhi);
- rho23 = v23/(dDip*dPhi);
- dPop = 100. *dLop * pTotal;
- dDip = 1000. * dDip;
- dPhi = 1000. * dPhi;
- dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);
- // resol1[istep][in] = dPop;
- // resol2[istep][in] = dDip;
- // resol3[istep][in] = dPhi;
- // resol4[istep][in] = dPPTot;
- // resol5[istep][in] = dPPStrag;
- sumDPop = sumDPop + dPop;
- sumDDip = sumDDip + dDip;
- sumDPhi = sumDPhi + dPhi;
- isum = isum + 1;}
- else{
- printf("%20s %10.5f %10.5f %20s\n","pT,eta",pT,eta,"cannot smear");
- }
- }
- if(isum > 0){
- dPop = sumDPop/isum;
- dDip = sumDDip/isum;
- dPhi = sumDPhi/isum;
- dPPTot = TMath::Sqrt(dPop*dPop + dPPStrag*dPPStrag);}
- else{
- dPop = 0;
- dDip = 0;
- dPhi = 0;
- }
- store1[istep] = dPop;
- store2[istep] = dDip;
- store3[istep] = dPhi;
- store4[istep] = dPPTot;
- store5[istep] = dPPStrag;
- if(istep > 20 ){
- Int_t im = 5;
- if(istep > 100) {im = 25;}
- if(istep > 500) {im = 100;}
- fm = 1./(1.*im);
- for(Int_t ist=1; ist<im; ist++){
- // for(Int_t in=1; in < 11; in++){
- // resol1[istep-im+ist][in] = resol1[istep-im][in]+
- // ist*fm*(resol1[istep][in]-resol1[istep-im][in]);
- // resol2[istep-im+ist][in] = resol2[istep-im][in]+
- // ist*fm*(resol2[istep][in]-resol2[istep-im][in]);
- // resol3[istep-im+ist][in] = resol3[istep-im][in]+
- // ist*fm*(resol3[istep][in]-resol3[istep-im][in]);
- // resol4[istep-im+ist][in] = resol4[istep-im][in]+
- // ist*fm*(resol4[istep][in]-resol4[istep-im][in]);
- // resol5[istep-im+ist][in] = resol5[istep-im][in]+
- // ist*fm*(resol5[istep][in]-resol5[istep-im][in]);
- // }
- store1[istep-im+ist]=store1[istep-im]+
- ist*fm*(store1[istep]-store1[istep-im]);
- store2[istep-im+ist]=store2[istep-im]+
- ist*fm*(store2[istep]-store2[istep-im]);
- store3[istep-im+ist]=store3[istep-im]+
- ist*fm*(store3[istep]-store3[istep-im]);
- store4[istep-im+ist]=store4[istep-im]+
- ist*fm*(store4[istep]-store4[istep-im]);
- store5[istep-im+ist]=store5[istep-im]+
- ist*fm*(store5[istep]-store5[istep-im]);
- // Fill control histograms
- fResID1Test->Fill(pTStart + 0.01*(istep-im+ist),store1[istep-im+ist]);
- fResID2Test->Fill(pTStart + 0.01*(istep-im+ist),store2[istep-im+ist]);
- fResID3Test->Fill(pTStart + 0.01*(istep-im+ist),store3[istep-im+ist]);
- fResID4Test->Fill(pTStart + 0.01*(istep-im+ist),store4[istep-im+ist]);
- fResID5Test->Fill(pTStart + 0.01*(istep-im+ist),store5[istep-im+ist]);
- }
- printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n",
- "TestTrack:",istep,store1[istep],store2[istep],store3[istep],
- store4[istep],store5[istep]);
- } else {
- printf("%10s %10d %20.15f %20.15f %20.15f %20.15f %20.15f \n",
- "TestTrack:",istep,store1[istep],store2[istep],store3[istep],
- store4[istep],store5[istep]);
- fResID1Test->Fill(pT,store1[istep]);
- fResID2Test->Fill(pT,store2[istep]);
- fResID3Test->Fill(pT,store3[istep]);
- fResID4Test->Fill(pT,store4[istep]);
- fResID5Test->Fill(pT,store5[istep]);
- }
- }
-}
-//_____________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff