+ // create the chain for the tree of compared tracks
+ TChain ch1("cmptrktree");
+ TChain ch2("cmptrktree");
+ TChain ch3("cmptrktree");
+
+ for(Int_t j=evFirst; j<=evLast; j++) {
+ cerr<<"Processing event "<<j<<endl;
+
+ TString covName("CovMatrix.");
+ covName+=j;
+ covName.Append(".root");
+
+ if(gSystem->AccessPathName(covName.Data(),kFileExists)) continue;
+
+
+ if(j<=100) ch1.Add(covName.Data());
+ if(j>100 && j<=200) ch2.Add(covName.Data());
+ if(j>200) ch3.Add(covName.Data());
+
+ } // loop on events
+
+ // merge chain in one file
+ TFile *covOut=0;
+ covOut = new TFile("CovMatrix_AllEvts_1.root","recreate");
+ ch1.Merge(covOut,1000000000);
+ covOut->Close();
+ delete covOut;
+ covOut = new TFile("CovMatrix_AllEvts_2.root","recreate");
+ ch2.Merge(covOut,1000000000);
+ covOut->Close();
+ delete covOut;
+ covOut = new TFile("CovMatrix_AllEvts_3.root","recreate");
+ ch3.Merge(covOut,1000000000);
+ covOut->Close();
+ delete covOut;
+
+
+ // efficiencies
+ cerr<<" ***** EFFICIENCIES ******\n\n";
+
+ ReadEffs("TPCeff.1.root");
+
+ Int_t n = fEffPi.GetTotPoints();
+ Double_t *avEffPi = new Double_t[n];
+ Double_t *avEffKa = new Double_t[n];
+ Double_t *avEffPr = new Double_t[n];
+ Double_t *avEffEl = new Double_t[n];
+ Double_t *avEffMu = new Double_t[n];
+ Int_t *evtsPi = new Int_t[n];
+ Int_t *evtsKa = new Int_t[n];
+ Int_t *evtsPr = new Int_t[n];
+ Int_t *evtsEl = new Int_t[n];
+ Int_t *evtsMu = new Int_t[n];
+
+ for(Int_t j=0; j<n; j++) {
+ avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
+ evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
+ }
+
+ for(Int_t j=evFirst; j<=evLast; j++) {
+ cerr<<"Processing event "<<j<<endl;
+
+ TString effName("TPCeff.");
+ effName+=j;
+ effName.Append(".root");
+
+ if(gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
+
+ ReadEffs(effName.Data());
+
+ for(Int_t k=0; k<n; k++) {
+ if(fEffPi.GetParam(k)>=0.) {avEffPi[k]+=fEffPi.GetParam(k); evtsPi[k]++;}
+ if(fEffKa.GetParam(k)>=0.) {avEffKa[k]+=fEffKa.GetParam(k); evtsKa[k]++;}
+ if(fEffPr.GetParam(k)>=0.) {avEffPr[k]+=fEffPr.GetParam(k); evtsPr[k]++;}
+ if(fEffEl.GetParam(k)>=0.) {avEffEl[k]+=fEffEl.GetParam(k); evtsEl[k]++;}
+ if(fEffMu.GetParam(k)>=0.) {avEffMu[k]+=fEffMu.GetParam(k); evtsMu[k]++;}
+ }
+
+ } // loop on events
+
+ // compute average efficiencies
+ for(Int_t j=0; j<n; j++) {
+ if(evtsPi[j]==0) evtsPi[j]++;
+ fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
+ if(evtsKa[j]==0) evtsKa[j]++;
+ fEffKa.SetParam(j,(Double_t)avEffKa[j]/evtsKa[j]);
+ if(evtsPr[j]==0) evtsPr[j]++;
+ fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
+ if(evtsEl[j]==0) evtsEl[j]++;
+ fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
+ if(evtsMu[j]==0) evtsMu[j]++;
+ fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
+ }
+
+ // write efficiencies to a file
+ WriteEffs(outName);
+
+ delete [] avEffPi;
+ delete [] avEffKa;
+ delete [] avEffPr;
+ delete [] avEffEl;
+ delete [] avEffMu;
+ delete [] evtsPi;
+ delete [] evtsKa;
+ delete [] evtsPr;
+ delete [] evtsEl;
+ delete [] evtsMu;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadAllData(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads all parameters from the DB
+//-----------------------------------------------------------------------------
+
+ if(!ReadEffs(inName)) return 0;
+ if(!ReadPulls(inName)) return 0;
+ if(!ReadRegParams(inName,211)) return 0;
+ if(!ReadRegParams(inName,321)) return 0;
+ if(!ReadRegParams(inName,2212)) return 0;
+ if(!ReadRegParams(inName,11)) return 0;
+ if(!ReadRegParams(inName,13)) return 0;
+ if(!ReaddEdx(inName,211)) return 0;
+ if(!ReaddEdx(inName,321)) return 0;
+ if(!ReaddEdx(inName,2212)) return 0;
+ if(!ReaddEdx(inName,11)) return 0;
+ if(!ReaddEdx(inName,13)) return 0;
+ if(!ReadDBgrid(inName)) return 0;
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReaddEdx(const Char_t* inName,Int_t pdg) {
+//-----------------------------------------------------------------------------
+// This function reads the dEdx parameters from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReaddEdx: "<<inName<<" not found\n";
+ return 0; }
+
+ TFile *inFile = TFile::Open(inName);
+ switch (pdg) {
+ case 211:
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
+ fdEdxMeanPi.~AliTPCkineGrid();
+ new(&fdEdxMeanPi) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
+ fdEdxRMSPi.~AliTPCkineGrid();
+ new(&fdEdxRMSPi) AliTPCkineGrid(*fdEdxRMS);
+ break;
+ case 321:
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxMeanKa");
+ fdEdxMeanKa.~AliTPCkineGrid();
+ new(&fdEdxMeanKa) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Kaons/dEdxRMSKa");
+ fdEdxRMSKa.~AliTPCkineGrid();
+ new(&fdEdxRMSKa) AliTPCkineGrid(*fdEdxRMS);
+ break;
+ case 2212:
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxMeanPr");
+ fdEdxMeanPr.~AliTPCkineGrid();
+ new(&fdEdxMeanPr) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Protons/dEdxRMSPr");
+ fdEdxRMSPr.~AliTPCkineGrid();
+ new(&fdEdxRMSPr) AliTPCkineGrid(*fdEdxRMS);
+ break;
+ case 11:
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxMeanEl");
+ fdEdxMeanEl.~AliTPCkineGrid();
+ new(&fdEdxMeanEl) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Electrons/dEdxRMSEl");
+ fdEdxRMSEl.~AliTPCkineGrid();
+ new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
+ break;
+ case 13:
+ if(fColl==0) {
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
+ } else {
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxMeanMu");
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxRMSMu");
+ }
+ fdEdxMeanMu.~AliTPCkineGrid();
+ new(&fdEdxMeanMu) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMSMu.~AliTPCkineGrid();
+ new(&fdEdxRMSMu) AliTPCkineGrid(*fdEdxRMS);
+ break;
+ }
+ inFile->Close();
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadDBgrid(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the kine grid from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReadCovMatrices: "<<inName<<" not found\n";
+ return 0; }
+
+ TFile *inFile = TFile::Open(inName);
+
+ // first read the DB grid for the different particles
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
+ fDBgridPi.~AliTPCkineGrid();
+ new(&fDBgridPi) AliTPCkineGrid(*fDBgrid);
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
+ fDBgridKa.~AliTPCkineGrid();
+ new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
+ if(fColl==0) {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
+ } else {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Protons/DBgridPr");
+ }
+ fDBgridPr.~AliTPCkineGrid();
+ new(&fDBgridPr) AliTPCkineGrid(*fDBgrid);
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
+ fDBgridEl.~AliTPCkineGrid();
+ new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
+ if(fColl==0) {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
+ } else {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Muons/DBgridMu");
+ }
+ fDBgridMu.~AliTPCkineGrid();
+ new(&fDBgridMu) AliTPCkineGrid(*fDBgrid);
+
+ inFile->Close();
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadEffs(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the TPC efficiencies from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReadEffs: "<<inName<<" not found\n";
+ return 0; }
+
+ TFile *inFile = TFile::Open(inName);
+
+ fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Pions/EffPi");
+ fEffPi.~AliTPCkineGrid();
+ new(&fEffPi) AliTPCkineGrid(*fEff);
+ fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Kaons/EffKa");
+ fEffKa.~AliTPCkineGrid();
+ new(&fEffKa) AliTPCkineGrid(*fEff);
+ fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Protons/EffPr");
+ fEffPr.~AliTPCkineGrid();
+ new(&fEffPr) AliTPCkineGrid(*fEff);
+ fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Electrons/EffEl");
+ fEffEl.~AliTPCkineGrid();
+ new(&fEffEl) AliTPCkineGrid(*fEff);
+ fEff = (AliTPCkineGrid*)inFile->Get("/Efficiencies/Muons/EffMu");
+ fEffMu.~AliTPCkineGrid();
+ new(&fEffMu) AliTPCkineGrid(*fEff);
+
+ inFile->Close();
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadPulls(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the pulls from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReadPulls: "<<inName<<" not found\n";
+ return 0; }
+
+ TFile *inFile = TFile::Open(inName);
+
+ for(Int_t i=0; i<5; i++) {
+ TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
+ TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
+ TString pr("/Pulls/Protons/PullsPr_"); pr+=i;
+ TString el("/Pulls/Electrons/PullsEl_"); el+=i;
+ TString mu("/Pulls/Muons/PullsMu_"); mu+=i;
+
+ fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
+ fPullsPi[i].~AliTPCkineGrid();
+ new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
+
+ fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
+ fPullsKa[i].~AliTPCkineGrid();
+ new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
+
+ if(fColl==0) {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
+ } else {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pr.Data());
+ }
+ fPullsPr[i].~AliTPCkineGrid();
+ new(&fPullsPr[i]) AliTPCkineGrid(*fPulls);
+
+ fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
+ fPullsEl[i].~AliTPCkineGrid();
+ new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
+
+ if(fColl==0) {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
+ } else {
+ fPulls = (AliTPCkineGrid*)inFile->Get(mu.Data());
+ }
+ fPullsMu[i].~AliTPCkineGrid();
+ new(&fPullsMu[i]) AliTPCkineGrid(*fPulls);
+ }
+
+ inFile->Close();
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) {
+//-----------------------------------------------------------------------------
+// This function reads the regularization parameters from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
+ return 0; }
+
+ TFile *inFile = TFile::Open(inName);
+ switch (pdg) {
+ case 211:
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
+ new(&fRegParPi) TMatrixD(*fRegPar);
+ break;
+ case 321:
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons");
+ new(&fRegParKa) TMatrixD(*fRegPar);
+ break;
+ case 2212:
+ if(fColl==0) {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
+ } else {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
+ }
+ new(&fRegParPr) TMatrixD(*fRegPar);
+ break;
+ case 11:
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
+ new(&fRegParEl) TMatrixD(*fRegPar);
+ break;
+ case 13:
+ if(fColl==0) {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
+ } else {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
+ }
+ new(&fRegParMu) TMatrixD(*fRegPar);
+ break;
+ }
+ inFile->Close();
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::RegularizeCovMatrix(const Char_t *outName,Int_t pdg) {
+//-----------------------------------------------------------------------------
+// This function regularizes the elements of the covariance matrix
+// that show a momentum depence:
+// c00, c11, c22, c33, c44, c20, c24, c40, c31
+// The regularization is done separately for pions, kaons, electrons:
+// give "Pion","Kaon" or "Electron" as first argument.
+//-----------------------------------------------------------------------------
+
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(10001);
+
+ Int_t thisPdg=211;
+ Char_t *part="Pions - ";
+
+ InitializeKineGrid("DB");
+ SetParticle(pdg);
+ const Int_t fitbins = fDBgrid->GetBinsPt();
+ cerr<<" Fit bins: "<<fitbins<<endl;
+
+ switch (pdg) {
+ case 211: // pions
+ thisPdg=211;
+ part="Pions - ";
+ cerr<<" Processing pions ...\n";
+ break;
+ case 321: //kaons
+ thisPdg=321;
+ part="Kaons - ";
+ cerr<<" Processing kaons ...\n";
+ break;
+ case 2212: //protons
+ thisPdg=2212;
+ part="Protons - ";
+ cerr<<" Processing protons ...\n";
+ break;
+ case 11: // electrons
+ thisPdg= 11;
+ part="Electrons - ";
+ cerr<<" Processing electrons ...\n";
+ break;
+ case 13: // muons
+ thisPdg= 13;
+ part="Muons - ";
+ cerr<<" Processing muons ...\n";
+ break;
+ }
+
+
+ /*
+ // create a chain with compared tracks
+ TChain *cmptrkchain = new ("cmptrktree");
+ cmptrkchain.Add("CovMatrix_AllEvts.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ */
+
+
+ TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
+
+ COMPTRACK cmptrk;
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
+
+
+ Int_t pbin;
+ Int_t *n = new Int_t[fitbins];
+ Int_t *n00 = new Int_t[fitbins];
+ Int_t *n11 = new Int_t[fitbins];
+ Int_t *n20 = new Int_t[fitbins];
+ Int_t *n22 = new Int_t[fitbins];
+ Int_t *n31 = new Int_t[fitbins];
+ Int_t *n33 = new Int_t[fitbins];
+ Int_t *n40 = new Int_t[fitbins];
+ Int_t *n42 = new Int_t[fitbins];
+ Int_t *n44 = new Int_t[fitbins];
+ Double_t *p = new Double_t[fitbins];
+ Double_t *ep = new Double_t[fitbins];
+ Double_t *mean00 = new Double_t[fitbins];
+ Double_t *mean11 = new Double_t[fitbins];
+ Double_t *mean20 = new Double_t[fitbins];
+ Double_t *mean22 = new Double_t[fitbins];
+ Double_t *mean31 = new Double_t[fitbins];
+ Double_t *mean33 = new Double_t[fitbins];
+ Double_t *mean40 = new Double_t[fitbins];
+ Double_t *mean42 = new Double_t[fitbins];
+ Double_t *mean44 = new Double_t[fitbins];
+ Double_t *sigma00 = new Double_t[fitbins];
+ Double_t *sigma11 = new Double_t[fitbins];
+ Double_t *sigma20 = new Double_t[fitbins];
+ Double_t *sigma22 = new Double_t[fitbins];
+ Double_t *sigma31 = new Double_t[fitbins];
+ Double_t *sigma33 = new Double_t[fitbins];
+ Double_t *sigma40 = new Double_t[fitbins];
+ Double_t *sigma42 = new Double_t[fitbins];
+ Double_t *sigma44 = new Double_t[fitbins];
+ Double_t *rmean = new Double_t[fitbins];
+ Double_t *rsigma = new Double_t[fitbins];
+ Double_t fitpar[3];
+
+ for(Int_t l=0; l<fitbins; l++) {
+ n[l]=1;
+ n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
+ p[l ]=ep[l]=0.;
+ mean00[l]=mean11[l]=mean20[l]=mean22[l]=mean31[l]=mean33[l]=mean40[l]=mean42[l]=mean44[l]=0.;
+ sigma00[l]=sigma11[l]=sigma20[l]=sigma22[l]=sigma31[l]=sigma33[l]=sigma40[l]=sigma42[l]=sigma44[l]=0.;
+ }
+
+ // loop on chain entries for mean
+ for(Int_t l=0; l<entries; l++) {
+ cmptrktree->GetEvent(l);
+ if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
+ pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
+ n[pbin]++;
+ p[pbin]+=cmptrk.p;
+ mean00[pbin]+=cmptrk.c00;
+ mean11[pbin]+=cmptrk.c11;
+ mean20[pbin]+=cmptrk.c20;
+ mean22[pbin]+=cmptrk.c22;
+ mean31[pbin]+=cmptrk.c31;
+ mean33[pbin]+=cmptrk.c33;
+ mean40[pbin]+=cmptrk.c40;
+ mean42[pbin]+=cmptrk.c42;
+ mean44[pbin]+=cmptrk.c44;
+ } // loop on chain entries
+
+ for(Int_t l=0; l<fitbins; l++) {
+ p[l]/=n[l];
+ mean00[l]/=n[l];
+ mean11[l]/=n[l];
+ mean20[l]/=n[l];
+ mean22[l]/=n[l];
+ mean31[l]/=n[l];
+ mean33[l]/=n[l];
+ mean40[l]/=n[l];
+ mean42[l]/=n[l];
+ mean44[l]/=n[l];
+ }
+
+ // loop on chain entries for sigma
+ for(Int_t l=0; l<entries; l++) {
+ cmptrktree->GetEvent(l);
+ if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
+ pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
+ if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
+ sigma00[pbin]+=(cmptrk.c00-mean00[pbin])*(cmptrk.c00-mean00[pbin]); }
+ if(TMath::Abs(cmptrk.c11-mean11[pbin])<0.4*mean11[pbin]) { n11[pbin]++;
+ sigma11[pbin]+=(cmptrk.c11-mean11[pbin])*(cmptrk.c11-mean11[pbin]); }
+ if(TMath::Abs(cmptrk.c20-mean20[pbin])<0.4*mean20[pbin]) { n20[pbin]++;
+ sigma20[pbin]+=(cmptrk.c20-mean20[pbin])*(cmptrk.c20-mean20[pbin]); }
+ if(TMath::Abs(cmptrk.c22-mean22[pbin])<0.4*mean22[pbin]) { n22[pbin]++;
+ sigma22[pbin]+=(cmptrk.c22-mean22[pbin])*(cmptrk.c22-mean22[pbin]); }
+ if(TMath::Abs(cmptrk.c31-mean31[pbin])<-0.4*mean31[pbin]) { n31[pbin]++;
+ sigma31[pbin]+=(cmptrk.c31-mean31[pbin])*(cmptrk.c31-mean31[pbin]); }
+ if(TMath::Abs(cmptrk.c33-mean33[pbin])<0.4*mean33[pbin]) { n33[pbin]++;
+ sigma33[pbin]+=(cmptrk.c33-mean33[pbin])*(cmptrk.c33-mean33[pbin]); }
+ if(TMath::Abs(cmptrk.c40-mean40[pbin])<0.4*mean40[pbin]) { n40[pbin]++;
+ sigma40[pbin]+=(cmptrk.c40-mean40[pbin])*(cmptrk.c40-mean40[pbin]); }
+ if(TMath::Abs(cmptrk.c42-mean42[pbin])<0.4*mean42[pbin]) { n42[pbin]++;
+ sigma42[pbin]+=(cmptrk.c42-mean42[pbin])*(cmptrk.c42-mean42[pbin]); }
+ if(TMath::Abs(cmptrk.c44-mean44[pbin])<0.4*mean44[pbin]) { n44[pbin]++;
+ sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
+ } // loop on chain entries
+
+ for(Int_t l=0; l<fitbins; l++) {
+ sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
+ sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
+ sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
+ sigma22[l] = TMath::Sqrt(sigma22[l]/n22[l]);
+ sigma31[l] = TMath::Sqrt(sigma31[l]/n31[l]);
+ sigma33[l] = TMath::Sqrt(sigma33[l]/n33[l]);
+ sigma40[l] = TMath::Sqrt(sigma40[l]/n40[l]);
+ sigma42[l] = TMath::Sqrt(sigma42[l]/n42[l]);
+ sigma44[l] = TMath::Sqrt(sigma44[l]/n44[l]);
+ }
+
+
+ // Fit function
+ TF1 *func = new TF1("RegFunc",RegFunc,0.23,50.,3);
+ func->SetParNames("A_meas","A_scatt","B");
+
+ // line to draw on the plots
+ TLine *lin = new TLine(-1,1,1.69,1);
+ lin->SetLineStyle(2);
+ lin->SetLineWidth(2);
+
+ // matrix used to store fit results
+ TMatrixD fitRes(9,3);
+
+ // --- c00 ---
+
+ // create the canvas
+ TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
+ canv00->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
+ TString title00("C(y,y)"); title00.Prepend(part);
+ TH2F *frame00 = new TH2F("frame00",title00.Data(),2,0.1,50,2,0,5e-3);
+ frame00->SetXTitle("p [GeV/c]");
+ canv00->cd(1); gPad->SetLogx();
+ frame00->Draw();
+ gr00->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(1.6e-3,1.9e-4,1.5);
+ // Fit points in range defined by function
+ gr00->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean00[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma00[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B})");
+ regtitle00.Prepend(part);
+ TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
+ frame00reg->SetXTitle("p [GeV/c]");
+ canv00->cd(2); gPad->SetLogx();
+ frame00reg->Draw();
+ gr00reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c11 ---
+
+ // create the canvas
+ TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
+ canv11->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
+ TString title11("C(z,z)"); title11.Prepend(part);
+ TH2F *frame11 = new TH2F("frame11",title11.Data(),2,0.1,50,2,0,6e-3);
+ frame11->SetXTitle("p [GeV/c]");
+ canv11->cd(1); gPad->SetLogx();
+ frame11->Draw();
+ gr11->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(1.2e-3,8.1e-4,1.);
+ // Fit points in range defined by function
+ gr11->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean11[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma11[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B})");
+ regtitle11.Prepend(part);
+ TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
+ frame11reg->SetXTitle("p [GeV/c]");
+ canv11->cd(2); gPad->SetLogx();
+ frame11reg->Draw();
+ gr11reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c20 ---
+
+ // create the canvas
+ TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
+ canv20->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
+ TString title20("C(#eta, y)"); title20.Prepend(part);
+ TH2F *frame20 = new TH2F("frame20",title20.Data(),2,0.1,50,2,0,2.5e-4);
+ frame20->SetXTitle("p [GeV/c]");
+ canv20->cd(1); gPad->SetLogx();
+ frame20->Draw();
+ gr20->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(7.3e-5,1.2e-5,1.5);
+ // Fit points in range defined by function
+ gr20->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean20[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma20[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B})");
+ regtitle20.Prepend(part);
+ TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
+ frame20reg->SetXTitle("p [GeV/c]");
+ canv20->cd(2); gPad->SetLogx();
+ frame20reg->Draw();
+ gr20reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c22 ---
+
+ // create the canvas
+ TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
+ canv22->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
+ TString title22("C(#eta, #eta)"); title22.Prepend(part);
+ TH2F *frame22 = new TH2F("frame22",title22.Data(),2,0.1,50,2,0,3e-5);
+ frame22->SetXTitle("p [GeV/c]");
+ canv22->cd(1); gPad->SetLogx();
+ frame22->Draw();
+ gr22->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(5.2e-6,1.1e-6,2.);
+ // Fit points in range defined by function
+ gr22->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean22[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma22[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B})");
+ regtitle22.Prepend(part);
+ TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
+ frame22reg->SetXTitle("p [GeV/c]");
+ canv22->cd(2); gPad->SetLogx();
+ frame22reg->Draw();
+ gr22reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c31 ---
+
+ // create the canvas
+ TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
+ canv31->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
+ TString title31("C(tg #lambda,z)"); title31.Prepend(part);
+ TH2F *frame31 = new TH2F("frame31",title31.Data(),2,0.1,50,2,-2e-4,0);
+ frame31->SetXTitle("p [GeV/c]");
+ canv31->cd(1); gPad->SetLogx();
+ frame31->Draw();
+ gr31->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(-1.2e-5,-1.2e-5,1.5);
+ // Fit points in range defined by function
+ gr31->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean31[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = -sigma31[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B})");
+ regtitle31.Prepend(part);
+ TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
+ frame31reg->SetXTitle("p [GeV/c]");
+ canv31->cd(2); gPad->SetLogx();
+ frame31reg->Draw();
+ gr31reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c33 ---
+
+ // create the canvas
+ TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
+ canv33->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
+ TString title33("C(tg #lambda,tg #lambda)"); title33.Prepend(part);
+ TH2F *frame33 = new TH2F("frame33",title33.Data(),2,0.1,50,2,0,1e-5);
+ frame33->SetXTitle("p [GeV/c]");
+ canv33->cd(1); gPad->SetLogx();
+ frame33->Draw();
+ gr33->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(1.3e-7,4.6e-7,1.7);
+ // Fit points in range defined by function
+ gr33->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean33[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma33[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B})");
+ regtitle33.Prepend(part);
+ TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
+ frame33reg->SetXTitle("p [GeV/c]");
+ canv33->cd(2); gPad->SetLogx();
+ frame33reg->Draw();
+ gr33reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c40 ---
+
+ // create the canvas
+ TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
+ canv40->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
+ TString title40("C(C,y)"); title40.Prepend(part);
+ TH2F *frame40 = new TH2F("frame40",title40.Data(),2,0.1,50,2,0,1e-6);
+ frame40->SetXTitle("p [GeV/c]");
+ canv40->cd(1); gPad->SetLogx();
+ frame40->Draw();
+ gr40->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(4.e-7,4.4e-8,1.5);
+ // Fit points in range defined by function
+ gr40->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean40[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma40[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B})");
+ regtitle40.Prepend(part);
+ TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
+ frame40reg->SetXTitle("p [GeV/c]");
+ canv40->cd(2); gPad->SetLogx();
+ frame40reg->Draw();
+ gr40reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c42 ---
+
+ // create the canvas
+ TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
+ canv42->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
+ TString title42("C(C, #eta)"); title42.Prepend(part);
+ TH2F *frame42 = new TH2F("frame42",title42.Data(),2,0.1,50,2,0,2.2e-7);
+ frame42->SetXTitle("p [GeV/c]");
+ canv42->cd(1); gPad->SetLogx();
+ frame42->Draw();
+ gr42->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(3.e-8,8.2e-9,2.);
+ // Fit points in range defined by function
+ gr42->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean42[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma42[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B})");
+ regtitle42.Prepend(part);
+ TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
+ frame42reg->SetXTitle("p [GeV/c]");
+ canv42->cd(2); gPad->SetLogx();
+ frame42reg->Draw();
+ gr42reg->Draw("P");
+ lin->Draw("same");
+
+
+ // --- c44 ---
+
+ // create the canvas
+ TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
+ canv44->Divide(1,2);
+ // create the graph for cov matrix
+ TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
+ TString title44("C(C,C)"); title44.Prepend(part);
+ TH2F *frame44 = new TH2F("frame44",title44.Data(),2,0.1,50,2,0,2e-9);
+ frame44->SetXTitle("p [GeV/c]");
+ canv44->cd(1); gPad->SetLogx();
+ frame44->Draw();
+ gr44->Draw("P");
+ // Sets initial values for parameters
+ func->SetParameters(1.8e-10,5.8e-11,2.);
+ // Fit points in range defined by function
+ gr44->Fit("RegFunc","R,Q");
+ func->GetParameters(fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
+ for(Int_t l=0; l<fitbins; l++) {
+ rmean[l] = mean44[l]/RegFunc(&p[l],fitpar);
+ rsigma[l] = sigma44[l]/RegFunc(&p[l],fitpar);
+ }
+ // create the graph the regularized cov. matrix
+ TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B})");
+ regtitle44.Prepend(part);
+ TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
+ frame44reg->SetXTitle("p [GeV/c]");
+ canv44->cd(2); gPad->SetLogx();
+ frame44reg->Draw();
+ gr44reg->Draw("P");
+ lin->Draw("same");
+
+
+
+
+ switch (pdg) {
+ case 211:
+ new(&fRegParPi) TMatrixD(fitRes);
+ break;
+ case 321:
+ new(&fRegParKa) TMatrixD(fitRes);
+ break;
+ case 2212:
+ new(&fRegParPr) TMatrixD(fitRes);
+ break;
+ case 11:
+ new(&fRegParEl) TMatrixD(fitRes);
+ break;
+ case 13:
+ new(&fRegParMu) TMatrixD(fitRes);
+ break;
+ }
+
+ // write fit parameters to file
+ WriteRegParams(outName,pdg);
+
+ delete [] n;
+ delete [] n00;
+ delete [] n11;
+ delete [] n20;
+ delete [] n22;
+ delete [] n31;
+ delete [] n33;
+ delete [] n40;
+ delete [] n42;
+ delete [] n44;
+ delete [] p;
+ delete [] ep;
+ delete [] mean00;
+ delete [] mean11;
+ delete [] mean20;
+ delete [] mean22;
+ delete [] mean31;
+ delete [] mean33;
+ delete [] mean40;
+ delete [] mean42;
+ delete [] mean44;
+ delete [] sigma00;
+ delete [] sigma11;
+ delete [] sigma20;
+ delete [] sigma22;
+ delete [] sigma31;
+ delete [] sigma33;
+ delete [] sigma40;
+ delete [] sigma42;
+ delete [] sigma44;
+ delete [] rmean;
+ delete [] rsigma;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+Bool_t AliTPCtrackerParam::SelectedTrack(Double_t pt,Double_t eta) const {
+//-----------------------------------------------------------------------------
+// This function makes a selection according to TPC tracking efficiency
+//-----------------------------------------------------------------------------
+
+ Double_t eff=0.;
+
+ eff = fEff->GetValueAt(pt,eta);
+
+ if(gRandom->Rndm() < eff) return kTRUE;
+
+ return kFALSE;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::SetParticle(Int_t pdg) {
+//-----------------------------------------------------------------------------
+// This function set efficiencies, pulls, etc... for the particle
+// specie of the current particle
+//-----------------------------------------------------------------------------
+
+ switch (pdg) {
+ case 211:
+ fDBgrid = &fDBgridPi;
+ fEff = &fEffPi;
+ fPulls = fPullsPi;
+ fRegPar = &fRegParPi;
+ fdEdxMean = &fdEdxMeanPi;
+ fdEdxRMS = &fdEdxRMSPi;
+ break;
+ case 321:
+ fDBgrid = &fDBgridKa;
+ fEff = &fEffKa;
+ fPulls = fPullsKa;
+ fRegPar = &fRegParKa;
+ fdEdxMean = &fdEdxMeanKa;
+ fdEdxRMS = &fdEdxRMSKa;
+ break;
+ case 2212:
+ fDBgrid = &fDBgridPr;
+ fEff = &fEffPr;
+ fPulls = fPullsPr;
+ fRegPar = &fRegParPr;
+ fdEdxMean = &fdEdxMeanPr;
+ fdEdxRMS = &fdEdxRMSPr;
+ break;
+ case 11:
+ fDBgrid = &fDBgridEl;
+ fEff = &fEffEl;
+ fPulls = fPullsEl;
+ fRegPar = &fRegParEl;
+ fdEdxMean = &fdEdxMeanEl;
+ fdEdxRMS = &fdEdxRMSEl;
+ break;
+ case 13:
+ fDBgrid = &fDBgridMu;
+ fEff = &fEffMu;
+ fPulls = fPullsMu;
+ fRegPar = &fRegParMu;
+ fdEdxMean = &fdEdxMeanMu;
+ fdEdxRMS = &fdEdxRMSMu;
+ break;
+ }
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov)
+ const {
+//-----------------------------------------------------------------------------
+// This function smears track parameters according to streched cov. matrix
+//-----------------------------------------------------------------------------