X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCtrackerParam.cxx;h=693e8da5ad9399e782e04ddac6a96d9618ef3ee3;hb=72600597e34e128f6b1e58cc3ebd220aeec02a5b;hp=fc666be492d689924801df4484b92e9536185e5f;hpb=e130146ca44b52012bc61405243e6a585759c0e3;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCtrackerParam.cxx b/TPC/AliTPCtrackerParam.cxx index fc666be492d..693e8da5ad9 100644 --- a/TPC/AliTPCtrackerParam.cxx +++ b/TPC/AliTPCtrackerParam.cxx @@ -11,261 +11,1064 @@ * 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. * - ************************************************************************** - ************************************************************************** + **************************************************************************/ + +/* $Id$ */ + +/************************************************************************** * * * This class builds AliTPCtrack objects from generated tracks to feed * * ITS tracking (V2). The AliTPCtrack is built from its first hit in * * the TPC. The track is assigned a Kalman-like covariance matrix * * depending on its pT and pseudorapidity and track parameters are * * smeared according to this covariance matrix. * - * Output file contains sorted tracks, ready for matching with ITS * + * Output file contains sorted tracks, ready for matching with ITS. * * * * For details: * - * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm * + * Alice Internal Note 2003-011 * * * * Test macro is: AliBarrelRec_TPCparam.C * * * - * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it * + * 2002/10/01: Introduction of the database for pp collisions (B=0.4 T) * + * - Everything done separately for pions, kaons, protons, electrons and * + * muons. * + * - Now (only for pp) the tracks are built from the AliTrackReferences * + * which contain the position and momentum of all tracks at R = 83 cm; * + * This eliminates the loss of tracks due the dead zone of the TPC * + * where the 1st hit is not created. * + * - In AliBarrelRec_TPCparam.C there many possible ways of determining * + * the z coordinate of the primary vertex in pp events (from pixels, * + * from ITS tracks, smearing according to resolution given by tracks. * + * * + * 2002/04/28: Major upgrade of the class * + * - Covariance matrices and pulls are now separeted for pions, kaons * + * and electrons. * + * - A parameterization of dE/dx in the TPC has been included and it is * + * used to assign a mass to each track according to a rough dE/dx PID. * + * - All the "numbers" have been moved to the file with the DataBase, * + * they are read as objects of the class AliTPCkineGrid, and assigned * + * to data memebers of the class AliTPCtrackerParam. * + * - All the code necessary to create a BataBase has been included in * + * class (see the macro AliTPCtrackingParamDB.C for the details). * + * * + * 2006/03/16: Adapted to ESD input/output * * * + * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it * + * adapted to ESD output by Marcello Lunardon, Padova * **************************************************************************/ -#include "AliTPCtrackerParam.h" -#include "alles.h" +// * +// This is a dummy comment +// +// +// * +//------- Root headers -------- +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +//------ AliRoot headers ------ +#include "AliGausCorr.h" +#include "AliTracker.h" +#include "AliMC.h" #include "AliMagF.h" +#include "AliRun.h" +#include "AliRunLoader.h" +#include "AliTPC.h" +#include "AliTPCParamSR.h" +#include "AliTPCkineGrid.h" #include "AliTPCtrack.h" -#include "TMatrixD.h" -#include "AliKalmanTrack.h" -#include "AliMagFCM.h" -#include "AliGausCorr.h" +#include "AliTPCtrackerParam.h" +#include "AliTrackReference.h" +//----------------------------- + +Double_t RegFunc(Double_t *x,Double_t *par) { +// This is the function used to regularize the covariance matrix + Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2]); + return value; +} +// structure for DB building +typedef struct { + Int_t pdg; + Int_t bin; + Double_t r; + Double_t p; + Double_t pt; + Double_t cosl; + Double_t eta; + Double_t dpt; + Double_t dP0,dP1,dP2,dP3,dP4; + Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44; + Double_t dEdx; +} COMPTRACK; +// cov matrix structure +typedef struct { + Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44; +} COVMATRIX; ClassImp(AliTPCtrackerParam) -//----------------------------------------------------------------- -AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz) +//----------------------------------------------------------------------------- +AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz, + const char* evfoldname):TObject(), + fEvFolderName(evfoldname), + fBz(kBz), + fColl(kcoll), + fSelAndSmear(kTRUE), + fDBfileName(""), + fTrack(), + fCovTree(0), + fDBgrid(0), + fDBgridPi(), + fDBgridKa(), + fDBgridPr(), + fDBgridEl(), + fDBgridMu(), + fEff(0), + fEffPi(), + fEffKa(), + fEffPr(), + fEffEl(), + fEffMu(), + fPulls(0), + fRegPar(0), + fRegParPi(), + fRegParKa(), + fRegParPr(), + fRegParEl(), + fRegParMu(), + fdEdxMean(0), + fdEdxMeanPi(), + fdEdxMeanKa(), + fdEdxMeanPr(), + fdEdxMeanEl(), + fdEdxMeanMu(), + fdEdxRMS(0), + fdEdxRMSPi(), + fdEdxRMSKa(), + fdEdxRMSPr(), + fdEdxRMSEl(), + fdEdxRMSMu() { -//----------------------------------------------------------------- +//----------------------------------------------------------------------------- // This is the class conctructor -//----------------------------------------------------------------- +//----------------------------------------------------------------------------- - fColl = coll; // collision code (0: PbPb6000) - fBz = Bz; // value of the z component of L3 field (Tesla) - -} -//----------------------------------------------------------------- -AliTPCtrackerParam::~AliTPCtrackerParam() -{} -//----------------------------------------------------------------- - -Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n) -{ -//----------------------------------------------------------------- -// This function creates the TPC parameterized tracks -//----------------------------------------------------------------- + // fBz = kBz; // value of the z component of L3 field (Tesla) + // fColl = kcoll; // collision code (0: PbPb6000; 1: pp) + // fSelAndSmear = kTRUE; // by default selection and smearing are done - if(fColl!=0) { - cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid collision!\n"; - cerr<<" Available: 0 -> PbPb6000"< PbPb6000\n 1 -> pp"); } - TFile *infile=(TFile*)inp; + fDBfileName = gSystem->Getenv("ALICE_ROOT"); + fDBfileName.Append("/TPC/CovMatrixDB_"); + //fDBfileName = "CovMatrixDB_"; + if(fColl==0) fDBfileName.Append("PbPb6000"); + if(fColl==1) fDBfileName.Append("pp"); + if(fBz==0.4) fDBfileName.Append("_B0.4T.root"); + // use same DB for 0.4 and 0.5 T; for 0.5 T, correction done in CookTrack() + if(fBz==0.5) fDBfileName.Append("_B0.4T.root"); +} +//----------------------------------------------------------------------------- +AliTPCtrackerParam::~AliTPCtrackerParam() {} +//____________________________________________________________________________ +AliTPCtrackerParam::AliTPCtrackerParam( const AliTPCtrackerParam& p) + :TObject(p), + fEvFolderName(""), + fBz(0.), + fColl(0), + fSelAndSmear(0), + fDBfileName(""), + fTrack(), + fCovTree(0), + fDBgrid(0), + fDBgridPi(), + fDBgridKa(), + fDBgridPr(), + fDBgridEl(), + fDBgridMu(), + fEff(0), + fEffPi(), + fEffKa(), + fEffPr(), + fEffEl(), + fEffMu(), + fPulls(0), + fRegPar(0), + fRegParPi(), + fRegParKa(), + fRegParPr(), + fRegParEl(), + fRegParMu(), + fdEdxMean(0), + fdEdxMeanPi(), + fdEdxMeanKa(), + fdEdxMeanPr(), + fdEdxMeanEl(), + fdEdxMeanMu(), + fdEdxRMS(0), + fdEdxRMSPi(), + fdEdxRMSKa(), + fdEdxRMSPr(), + fdEdxRMSEl(), + fdEdxRMSMu() +{ + // dummy copy constructor +} +//---------------------------------------------------------------------------- +AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant( + Double_t x,Double_t y,Double_t z, + Double_t px,Double_t py,Double_t pz, + Int_t lab) + :TObject(), + fXg(x), + fYg(y), + fZg(z), + fPx(px), + fPy(py), + fPz(pz), + fAlpha(0.), + fLabel(lab), + fSector(0) + +{ +//---------------------------------------------------------------------------- +// Constructor of the geant seeds +//---------------------------------------------------------------------------- + + Double_t a = TMath::ATan2(y,x)*180./TMath::Pi(); + if(a<0) a += 360.; + fSector = (Int_t)(a/20.); + fAlpha = 10.+20.*fSector; + fAlpha /= 180.; + fAlpha *= TMath::Pi(); +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::Init() { +//----------------------------------------------------------------------------- +// This function reads the DB from file +//----------------------------------------------------------------------------- + + if(fSelAndSmear) { + printf("+++\n+++ Reading DataBase from:%s\n+++\n+++\n",fDBfileName.Data()); + // Read paramters from DB file + if(!ReadAllData(fDBfileName.Data())) { + printf("AliTPCtrackerParam::BuildTPCtracks: \ + Could not read data from DB\n\n"); return 1; + } + + } else printf("\n ! Creating ALL TRUE tracks at TPC inner radius !\n\n"); - // Get gAlice object from file - if(!(gAlice=(AliRun*)infile->Get("gAlice"))) { - cerr<<"gAlice has not been found on galice.root !\n"; - return 1; - } - AliMagFCM *fiel = (AliMagFCM*)gAlice->Field(); - Double_t fieval=(Double_t)fiel->SolenoidField()/10.; + // Check if value in the galice file is equal to selected one (fBz) + AliMagF *fiel = (AliMagF*)gAlice->Field(); + Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.); printf("Magnetic field is %6.2f Tesla\n",fieval); if(fBz!=fieval) { - cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!"<Branch("tracks","AliTPCtrack",&tpctrack,20000,0); + AliRunLoader* rl = AliRunLoader::GetRunLoader(fEvFolderName); + if (rl == 0x0) { + Error("BuildTPCtracks","Can not get Run Loader from event folder named %s", + fEvFolderName.Data()); + return 2; + } + AliLoader* tpcloader = rl->GetLoader("TPCLoader"); + if (tpcloader == 0x0) { + Error("BuildTPCtracks","Can not get TPC Loader from Run Loader."); + return 3; + } + + // Get gAlice object from file + if(!gAlice) rl->LoadgAlice(); + //rl->LoadHeader(); + rl->LoadKinematics(); + tpcloader->LoadHits("read"); + + if(!(gAlice=rl->GetAliRun())) { + printf("Cannot get gAlice from Run Loader !\n"); + return 1; + } + + // Get TPC detector + AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC"); + + rl->CdGAFile(); + AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60"); + if(digp){ + delete digp; + digp = new AliTPCParamSR(); + } + else digp=(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60"); + + if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; } + tpc->SetParam(digp); + + TParticle *part=0; + AliTPCseedGeant *seed=0; + AliTPCtrack *tpctrack=0; + Double_t sPt,sEta; + Int_t bin,label,pdg,charge; + Int_t tracks; + Int_t nParticles,nSeeds,arrentr; + //Int_t nSel=0,nAcc=0; + + Int_t evt=event->GetEventNumber(); + + tracks=0; + + // array for TPC tracks + TObjArray tArray(20000); + + // array for TPC seeds with geant info + TObjArray sArray(20000); + + // get the particles stack + nParticles = (Int_t)gAlice->GetEvent(evt); + + Bool_t *done = new Bool_t[nParticles]; + Int_t *pdgCodes = new Int_t[nParticles]; + + // loop on particles and store pdg codes + for(Int_t l=0; lGetMCApp()->Particle(l); + pdgCodes[l] = part->GetPdgCode(); + done[l] = kFALSE; + } + + printf("+++ Number of particles: %d\n",nParticles); - // array for TPC tracks - TObjArray tarray(20000); + // Create the seeds for the TPC tracks at the inner radius of TPC + if(fColl==0) { + // Get TreeH with hits + TTree *th = tpcloader->TreeH(); + MakeSeedsFromHits(tpc,th,sArray); + } else { + // Get TreeTR with track references + rl->LoadTrackRefs(); + TTree *ttr = rl->TreeTR(); + if (!ttr) return -3; + MakeSeedsFromRefs(ttr,sArray); + } + + nSeeds = sArray.GetEntries(); + printf("+++ Number of seeds: %d\n",nSeeds); + + // loop over entries in sArray + for(Int_t l=0; lGetLabel(); + + // check if this track has already been processed + if(done[label]) continue; + + // PDG code & electric charge + pdg = pdgCodes[label]; + if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; } + else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; } + else continue; + pdg = TMath::Abs(pdg); + if(pdg>3000) pdg=211; + + if(fSelAndSmear) SetParticle(pdg); + + sPt = seed->GetPt(); + sEta = seed->GetEta(); + + // Apply selection according to TPC efficiency + //if(TMath::Abs(pdg)==211) nAcc++; + if(fSelAndSmear && !SelectedTrack(sPt,sEta)) continue; + //if(TMath::Abs(pdg)==211) nSel++; + + // create AliTPCtrack object + BuildTrack(seed,charge); + + if(fSelAndSmear) { + bin = fDBgrid->GetBin(sPt,sEta); + switch (pdg) { + case 211: + //fCovTree = &(fCovTreePi[bin]); + fCovTree = fCovTreePi[bin]; + break; + case 321: + //fCovTree = &(fCovTreeKa[bin]); + fCovTree = fCovTreeKa[bin]; + break; + case 2212: + //fCovTree = &(fCovTreePr[bin]); + fCovTree = fCovTreePr[bin]; + break; + case 11: + //fCovTree = &(fCovTreeEl[bin]); + fCovTree = fCovTreeEl[bin]; + break; + case 13: + //fCovTree = &(fCovTreeMu[bin]); + fCovTree = fCovTreeMu[bin]; + break; + } + // deal with covariance matrix and smearing of parameters + CookTrack(sPt,sEta); + // assign the track a dE/dx and make a rough PID + CookdEdx(sPt,sEta); + } + + // put track in array + AliTPCtrack *iotrack = new AliTPCtrack(fTrack); + iotrack->SetLabel(label); + tArray.AddLast(iotrack); + // Mark track as "done" and register the pdg code + done[label] = kTRUE; + tracks++; + + } // loop over entries in sArray + + // sort array with TPC tracks (decreasing pT) + tArray.Sort(); + + // convert to AliESDtrack and write to AliESD + arrentr = tArray.GetEntriesFast(); + Int_t idx; + Double_t wgts[5]; + for(Int_t l=0; lGetMass()-0.9)<0.1) { + idx = 4; // proton + } else if(TMath::Abs(tpctrack->GetMass()-0.5)<0.1) { + idx = 3; // kaon + } else { + idx = 2; // pion + } + wgts[idx] = 1.; + ioESDtrack.SetESDpid(wgts); + event->AddTrack(&ioESDtrack); + } + + + delete [] done; + delete [] pdgCodes; + + printf("+++ Number of TPC tracks: %d\n",tracks); + //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<SetOptStat(0); + gStyle->SetOptFit(10001); + + const char *part="PIONS"; + Double_t ymax=500.; + + /* + // 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(); + cerr<<" Number of entries: "<GetEvent(evt); + SetParticle(pdg); - Bool_t done[500000]; - for(Int_t l=0; l<500000; l++) { done[l]=kFALSE; } + const Int_t knTotBins = fDBgrid->GetTotBins(); + cerr<<" Fit bins: "<GetDetector("TPC"); - Int_t ver = TPC->IsVersion(); - cerr<<"+++ TPC version "<Get("75x40_100x60"); - if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; } - TPC->SetParam(digp); + Int_t bin=0; + Int_t *n = new Int_t[knTotBins]; + Double_t *p = new Double_t[knTotBins]; + Double_t *ep = new Double_t[knTotBins]; + Double_t *mean = new Double_t[knTotBins]; + Double_t *sigma = new Double_t[knTotBins]; - // Get TreeH with hits - TTree *TH=gAlice->TreeH(); - Int_t ntracks=(Int_t)TH->GetEntries(); - cerr<<"+++\n+++ Number of particles in event "<ResetHits(); - TH->GetEvent(i); - // Get FirstHit - tpcHit=(AliTPChit*)TPC->FirstHit(-1); - for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) { - if(tpcHit->fQ !=0.) continue; - // Get particle momentum at hit - hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z(); - hPt=TMath::Sqrt(hPx*hPx+hPy*hPy); - // reject hits with PtTrack(); - // check if this track has already been processed - if(done[label]) continue; - // electric charge - Part = gAlice->Particle(label); - pdg = Part->GetPdgCode(); - if(pdg>200 || pdg==-11 || pdg==-13) { charge=1; } - else if(pdg<-200 || pdg==11 || pdg==13) { charge=-1; } - else continue; - - - if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break; - if(tpcHit->fQ != 0.) continue; - // Get global coordinates of hit - xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z(); - if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue; - - // Get TPC sector, Alpha angle and local coordinates - // printf("Sector %d\n",tpcHit->fSector); - digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha); - alpha = TMath::ATan2(sinAlpha,cosAlpha); - xl = xg*cosAlpha + yg*sinAlpha; - yl =-xg*sinAlpha + yg*cosAlpha; - zl = zg; - //printf("Alpha %f xl %f yl %f zl %f\n",Alpha,xl,yl,zl); - - // reject tracks which are not in the TPC acceptance - if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue; - - // Get bin in pT,eta - bin = GetBin(hPt,Part->Eta()); - - // Apply selection according to TPC efficiency - //if(TMath::Abs(pdg)==211) nAcc++; - if(!SelectedTrack(pdg,hPt,Part->Eta())) continue; - //if(TMath::Abs(pdg)==211) nSel++; - - // Mark track as "done" - done[label]=kTRUE; tracks++; - - // create AliTPCtrack object - tpctrack = BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge,label); + for(Int_t l=0; lGetEvent(l); + if(TMath::Abs(cmptrk.pdg)!=pdg) continue; + bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta); + p[bin] += cmptrk.p; + mean[bin] += cmptrk.dEdx; + n[bin]++; + } // loop on chain entries + + for(Int_t l=0; lGetEvent(l); + if(TMath::Abs(cmptrk.pdg)!=pdg) continue; + bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta); + if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue; + n[bin]++; + sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]); + } // loop on chain entries + + for(Int_t l=0; lSetXTitle("p [GeV/c]"); + frame->SetYTitle("dE/dx [a.u.]"); + canv->SetLogx(); + frame->Draw(); + gr->Draw("P"); + + switch(pdg) { + case 211: + for(Int_t i=0; iGet("cmptrktree"); + + COMPTRACK cmptrk; + cmptrktree->SetBranchAddress("comptracks",&cmptrk); + Int_t entries = (Int_t)cmptrktree->GetEntries(); + cerr<<" Number of entries: "<Fill(); + for(Int_t i=0;i<5;i++) { + pulls[i].~AliTPCkineGrid(); + new(&pulls[i]) AliTPCkineGrid(*(fPulls+i)); + } + nTotBins = fDBgrid->GetTotBins(); + cerr<<"nTotBins = "<SetName(hname); hDum->SetTitle(htitle); + hPulls0[i] = *hDum; + sprintf(hname,"hPulls1%d",i); + sprintf(htitle,"P1 pulls for bin %d",i); + hDum->SetName(hname); hDum->SetTitle(htitle); + hPulls1[i] = *hDum; + sprintf(hname,"hPulls2%d",i); + sprintf(htitle,"P2 pulls for bin %d",i); + hDum->SetName(hname); hDum->SetTitle(htitle); + hPulls2[i] = *hDum; + sprintf(hname,"hPulls3%d",i); + sprintf(htitle,"P3 pulls for bin %d",i); + hDum->SetName(hname); hDum->SetTitle(htitle); + hPulls3[i] = *hDum; + sprintf(hname,"hPulls4%d",i); + sprintf(htitle,"P4 pulls for bin %d",i); + hDum->SetName(hname); hDum->SetTitle(htitle); + hPulls4[i] = *hDum; } - // write the tree with tracks in the output file - out->cd(); - tracktree->Write(); + // loop on chain entries + for(Int_t i=0; iGetEvent(i); + if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue; + // fill histograms with the pulls + bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta); + //cerr<<" pt "<10) { + g->SetRange(-3.*hPulls0[i].GetRMS(),3.*hPulls0[i].GetRMS()); + hPulls0[i].Fit("g","R,Q,N"); + pulls[0].SetParam(i,g->GetParameter(2)); + } else pulls[0].SetParam(i,-1.); + if(hPulls1[i].GetEntries()>10) { + g->SetRange(-3.*hPulls1[i].GetRMS(),3.*hPulls1[i].GetRMS()); + hPulls1[i].Fit("g","R,Q,N"); + pulls[1].SetParam(i,g->GetParameter(2)); + } else pulls[1].SetParam(i,-1.); + if(hPulls2[i].GetEntries()>10) { + g->SetRange(-3.*hPulls2[i].GetRMS(),3.*hPulls2[i].GetRMS()); + hPulls2[i].Fit("g","R,Q,N"); + pulls[2].SetParam(i,g->GetParameter(2)); + } else pulls[2].SetParam(i,-1.); + if(hPulls3[i].GetEntries()>10) { + g->SetRange(-3.*hPulls3[i].GetRMS(),3.*hPulls3[i].GetRMS()); + hPulls3[i].Fit("g","R,Q,N"); + pulls[3].SetParam(i,g->GetParameter(2)); + } else pulls[3].SetParam(i,-1.); + if(hPulls4[i].GetEntries()>10) { + g->SetRange(-3.*hPulls4[i].GetRMS(),3.*hPulls4[i].GetRMS()); + hPulls4[i].Fit("g","R,Q,N"); + pulls[4].SetParam(i,g->GetParameter(2)); + } else pulls[4].SetParam(i,-1.); + } // loop on bins + + + switch (part) { + case 0: // pions + for(Int_t i=0;i<5;i++) { + fPullsPi[i].~AliTPCkineGrid(); + new(&fPullsPi[i]) AliTPCkineGrid(pulls[i]); + } + break; + case 1: // kaons + for(Int_t i=0;i<5;i++) { + fPullsKa[i].~AliTPCkineGrid(); + new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]); + } + break; + case 2: // protons + for(Int_t i=0;i<5;i++) { + fPullsPr[i].~AliTPCkineGrid(); + new(&fPullsPr[i]) AliTPCkineGrid(pulls[i]); + } + break; + case 3: // electrons + for(Int_t i=0;i<5;i++) { + fPullsEl[i].~AliTPCkineGrid(); + new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]); + } + break; + case 4: // muons + for(Int_t i=0;i<5;i++) { + fPullsMu[i].~AliTPCkineGrid(); + new(&fPullsMu[i]) AliTPCkineGrid(pulls[i]); + //cerr<<" mu pulls "<Get("cmptrktree"); + + COMPTRACK cmptrk; + cmptrktree->SetBranchAddress("comptracks",&cmptrk); + Int_t entries = (Int_t)cmptrktree->GetEntries(); + cerr<<" Number of entries: "<GetPointsPt(); + cerr<<"knPtBins = "<GetArrayPt(pt); + + + TH1F *hDumP0 = new TH1F("nameP0","dy",100,-0.3,0.3); + TH1F *hDumP4 = new TH1F("nameP4","dC",100,-0.0005,0.0005); + TH1F *hDumPt = new TH1F("namePt","dp_{T}/p_{T}",100,-0.5,0.5); + + TF1 *g = new TF1("g","gaus"); + + // create histograms for the all the bins + TH1F *hP0=NULL; + TH1F *hP4=NULL; + TH1F *hPt=NULL; + + hP0 = new TH1F[knPtBins]; + hP4 = new TH1F[knPtBins]; + hPt = new TH1F[knPtBins]; + + for(Int_t i=0; iGetEvent(i); + if(TMath::Abs(cmptrk.pdg)!=pdg) continue; + // fill histograms with the residuals + bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta(); + //cerr<<" pt "<Divide(5,2); + TCanvas *cP4res = new TCanvas("cP4res","cP4res",0,0,1200,700); + cP4res->Divide(5,2); + TCanvas *cPtres = new TCanvas("cPtres","cPtres",0,0,1200,700); + cPtres->Divide(5,2); + + // Draw histograms + for(Int_t i=0; icd(i+1); hP0[i].Draw(); + cP4res->cd(i+1); hP4[i].Draw(); + cPtres->cd(i+1); hPt[i].Draw(); + } + + + // compute the sigma of the distributions + for(Int_t i=0; i10) { + g->SetRange(-3.*hP0[i].GetRMS(),3.*hP0[i].GetRMS()); + hP0[i].Fit("g","R,Q,N"); + dP0[i] = g->GetParameter(2); + } else dP0[i] = 0.; + if(hP4[i].GetEntries()>10) { + g->SetRange(-3.*hP4[i].GetRMS(),3.*hP4[i].GetRMS()); + hP4[i].Fit("g","R,Q,N"); + dP4[i] = g->GetParameter(2); + } else dP4[i] = 0.; + if(hPt[i].GetEntries()>10) { + g->SetRange(-3.*hPt[i].GetRMS(),3.*hPt[i].GetRMS()); + hPt[i].Fit("g","R,Q,N"); + dPtToPt[i] = 100.*g->GetParameter(2); + } else dPtToPt[i] = 0.; + } // loop on bins + + + TGraph *grdP0 = new TGraph(knPtBins,pt,dP0); + TGraph *grdP4 = new TGraph(knPtBins,pt,dP4); + TGraph *grdPtToPt = new TGraph(knPtBins,pt,dPtToPt); + + grdP0->SetMarkerStyle(20); grdP0->SetMarkerColor(2); grdP0->SetMarkerSize(1.5); + grdP4->SetMarkerStyle(21); grdP4->SetMarkerColor(3); grdP4->SetMarkerSize(1.5); + grdPtToPt->SetMarkerStyle(22); grdPtToPt->SetMarkerColor(4); grdPtToPt->SetMarkerSize(1.5); + + // Plot Results + gStyle->SetOptStat(0); + TCanvas *c1 = new TCanvas("c1","dP0",0,0,900,700); + c1->SetLogx(); + c1->SetGridx(); + c1->SetGridy(); + + TH2F *frame1 = new TH2F("frame1","y resolution VS p_{T} in TPC",2,0.1,30,2,0,0.1); + frame1->SetXTitle("p_{T} [GeV/c]"); + frame1->SetYTitle("#sigma(y) [cm]"); + frame1->Draw(); + grdP0->Draw("P"); + + + TCanvas *c2 = new TCanvas("c2","dP4",0,0,900,700); + c2->SetLogx(); + c2->SetGridx(); + c2->SetGridy(); + + TH2F *frame2 = new TH2F("frame2","C resolution VS p_{T} in TPC",2,0.1,30,2,0,0.0001); + frame2->SetXTitle("p_{T} [GeV/c]"); + frame2->SetYTitle("#sigma(C) [1/cm]"); + frame2->Draw(); + grdP4->Draw("P"); + + TCanvas *c3 = new TCanvas("c3","dPtToPt",0,0,900,700); + c3->SetLogx(); + c3->SetLogy(); + c3->SetGridx(); + c3->SetGridy(); + + TH2F *frame3 = new TH2F("frame3","Relative p_{T} resolution VS p_{T} in TPC",2,0.1,30,2,0.1,30.); + frame3->SetXTitle("p_{T} [GeV/c]"); + frame3->SetYTitle("dp_{T}/p_{T} (%)"); + frame3->Draw(); + grdPtToPt->Draw("P"); + + + delete [] dP0; + delete [] dP4; + delete [] dPtToPt; + delete [] pt; + + + delete [] hP0; + delete [] hP4; + delete [] hPt; + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) { +//----------------------------------------------------------------------------- // This function uses GEANT info to set true track parameters -//----------------------------------------------------------------- - Double_t xref = x; +//----------------------------------------------------------------------------- + Double_t xref = s->GetXL(); Double_t xx[5],cc[15]; cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.; cc[1]=cc[3]=cc[4]=cc[6]=cc[7]=cc[8]=cc[10]=cc[11]=cc[12]=cc[13]=0.; // Magnetic field - TVector3 bfield(0.,0.,fBz); + TVector3 bfield(0.,0.,-fBz); // radius [cm] of track projection in (x,y) - Double_t rho = pt*100./0.299792458/bfield.Z(); + Double_t rho = s->GetPt()*100./0.299792458/TMath::Abs(bfield.Z()); // center of track projection in local reference frame - TVector3 hmom,hpos; + TVector3 sMom,sPos; - // position (local) and momentum (local) at the hit + // position (local) and momentum (local) at the seed // in the bending plane (z=0) - hpos.SetXYZ(x,y,0.); - hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.); - TVector3 vrho = hmom.Cross(bfield); + sPos.SetXYZ(s->GetXL(),s->GetYL(),0.); + sMom.SetXYZ(s->GetPx()*TMath::Cos(s->GetAlpha())+s->GetPy()*TMath::Sin(s->GetAlpha()),-s->GetPx()*TMath::Sin(s->GetAlpha())+s->GetPy()*TMath::Cos(s->GetAlpha()),0.); + TVector3 vrho = sMom.Cross(bfield); vrho *= ch; vrho.SetMag(rho); - TVector3 vcenter = hpos+vrho; + TVector3 vcenter = sPos+vrho; Double_t x0 = vcenter.X(); @@ -273,135 +1076,618 @@ AliTPCtrack* AliTPCtrackerParam::BuildTrack(Double_t alpha,Double_t x, // fAlpha = Alpha Rotation angle the local (TPC sector) // fP0 = YL Y-coordinate of a track // fP1 = ZG Z-coordinate of a track - // fP2 = C*x0 x0 is center x in rotated frame + // fP2 = sin(phi) sine of the (local) azimuthal angle // fP3 = Tgl tangent of the track momentum dip angle // fP4 = C track curvature - xx[0] = y; - xx[1] = z; - xx[3] = pz/pt; - xx[4] = -ch/rho; - xx[2] = xx[4]*x0; + xx[0] = s->GetYL(); + xx[1] = s->GetZL(); + xx[2] = ch/rho*(xref-x0); + xx[3] = s->GetPz()/s->GetPt(); + xx[4] = ch/rho; // create the object AliTPCtrack - AliTPCtrack* track = new AliTPCtrack(0,xx,cc,xref,alpha); - // set the label - track->SetLabel(lab); + AliTPCtrack track(xref,s->GetAlpha(),xx,cc,0); + new(&fTrack) AliTPCtrack(track); - return track; + return; } +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::CompareTPCtracks( + Int_t nEvents, + const Char_t* galiceName, + const Char_t* trkGeaName, + const Char_t* trkKalName, + const Char_t* covmatName, + const Char_t* tpceffasciiName, + const Char_t* tpceffrootName) { +//----------------------------------------------------------------------------- +// This function compares tracks from TPC Kalman Filter V2 with +// geant tracks at TPC 1st hit. It gives: +// - a tree with Kalman cov. matrix elements, resolutions, dEdx +// - the efficiencies as a function of particle type, pT, eta +//----------------------------------------------------------------------------- + + TFile *kalFile = TFile::Open(trkKalName); + TFile *geaFile = TFile::Open(trkGeaName); + TFile *galiceFile = TFile::Open(galiceName); + + // get the AliRun object + AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice"); + + + // create the tree for comparison results + COMPTRACK cmptrk; + TTree *cmptrktree = new TTree("cmptrktree","results of track comparison"); + cmptrktree->Branch("comptracks",&cmptrk,"pdg/I:bin:r/D:p:pt:cosl:eta:dpt:dP0:dP1:dP2:dP3:dP4:c00:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44:dEdx"); + + InitializeKineGrid("eff"); + InitializeKineGrid("DB"); + Int_t effBins = fEffPi.GetTotPoints(); + Int_t effBinsPt = fEffPi.GetPointsPt(); + Double_t *pt = new Double_t[effBinsPt]; + fEffPi.GetArrayPt(pt); + + TParticle *part; + Double_t ptgener; + Bool_t usethis; + Int_t label; + Double_t dAlpha; + Int_t pi=0,ka=0,mu=0,el=0,pr=0; + Int_t *geaPi = new Int_t[effBins]; + Int_t *geaKa = new Int_t[effBins]; + Int_t *geaPr = new Int_t[effBins]; + Int_t *geaEl = new Int_t[effBins]; + Int_t *geaMu = new Int_t[effBins]; + Int_t *kalPi = new Int_t[effBins]; + Int_t *kalKa = new Int_t[effBins]; + Int_t *kalPr = new Int_t[effBins]; + Int_t *kalEl = new Int_t[effBins]; + Int_t *kalMu = new Int_t[effBins]; + Float_t *effPi = new Float_t[effBins]; + Float_t *effKa = new Float_t[effBins]; + Float_t *effPr = new Float_t[effBins]; + Float_t *effEl = new Float_t[effBins]; + Float_t *effMu = new Float_t[effBins]; + Int_t bin; + for(Int_t j=0; jGetEvent(evt); + + Int_t *kalLab = new Int_t[knparticles]; + for(Int_t i=0; iGet(tname); + if(!kaltree) continue; + AliTPCtrack *kaltrack=new AliTPCtrack; + kaltree->SetBranchAddress("tracks",&kaltrack); + Int_t kalEntries = (Int_t)kaltree->GetEntries(); + + // tracks from 1st hit + TTree *geatree=(TTree*)geaFile->Get(tname); + if(!geatree) continue; + AliTPCtrack *geatrack=new AliTPCtrack; + geatree->SetBranchAddress("tracks",&geatrack); + Int_t geaEntries = (Int_t)geatree->GetEntries(); + + cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<GetEvent(j); + if(kaltrack->GetLabel()>=0) { + if(kalLab[kaltrack->GetLabel()]!=-1) mult++; + kalLab[kaltrack->GetLabel()] = j; + } + else fake++; + } + + cerr<<"+++ Number of fake tracks in TPC Kalman: "<0) hasSeed[sLabel]=kTRUE; + } + fclose(seedFile); + */ + + cerr<<"Doing track comparison...\n"; + // loop on tracks at 1st hit + for(Int_t j=0; jGetEvent(j); + + label = geatrack->GetLabel(); + part = (TParticle*)gAlice->GetMCApp()->Particle(label); + + // use only injected tracks with fixed values of pT + ptgener = part->Pt(); + usethis = kFALSE; + for(Int_t l=0; lGetX()<GetY() + ( + TMath::Sqrt(1-(geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())* + (geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())) + -TMath::Sqrt(1-(geatrack->GetC()*x-geatrack->GetEta())* + (geatrack->GetC()*x-geatrack->GetEta())) + )/geatrack->GetC(); + y = TMath::Abs(y); + //cerr<<"y "<GetY()< 0.8*x*TMath::Tan(TMath::Pi()/18)) { out = kTRUE; break; } + } + if(out) continue; + */ + + cmptrk.pdg = part->GetPdgCode(); + cmptrk.eta = part->Eta(); + cmptrk.r = TMath::Sqrt(part->Vx()*part->Vx()+part->Vy()*part->Vy()); + + cmptrk.pt = 1/TMath::Abs(geatrack->Get1Pt()); + cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl())); + cmptrk.p = cmptrk.pt/cmptrk.cosl; + + + bin = fDBgridPi.GetBin(cmptrk.pt,cmptrk.eta); + cmptrk.bin = bin; + if(abs(cmptrk.pdg)==211) geaPi[bin]++; + if(abs(cmptrk.pdg)==321) geaKa[bin]++; + if(abs(cmptrk.pdg)==2212) geaPr[bin]++; + if(abs(cmptrk.pdg)==11) geaEl[bin]++; + if(abs(cmptrk.pdg)==13) geaMu[bin]++; + + + // check if there is the corresponding track in the TPC kalman and get it + if(kalLab[label]==-1) continue; + kaltree->GetEvent(kalLab[label]); + + // and go on only if it has xref = 84.57 cm (inner pad row) + if(kaltrack->GetX()>90.) continue; + + if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; } + if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; } + if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; } + if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; } + if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; } + + kaltrack->PropagateTo(geatrack->GetX(),1.e-9,0.); + + cmptrk.dEdx = kaltrack->GetdEdx(); + + // compute errors on parameters + dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha(); + if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; } + + cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY(); + cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ(); + cmptrk.dP2 = kaltrack->GetSnp()-geatrack->GetSnp(); + cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl(); + cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC(); + cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt(); + + // get covariance matrix + // beware: lines 3 and 4 in the matrix are inverted! + //kaltrack->GetCovariance(cc); + + cmptrk.c00 = kaltrack->GetSigmaY2(); + cmptrk.c10 = kaltrack->GetSigmaZY(); + cmptrk.c11 = kaltrack->GetSigmaZ2(); + cmptrk.c20 = kaltrack->GetSigmaSnpY(); + cmptrk.c21 = kaltrack->GetSigmaSnpY(); + cmptrk.c22 = kaltrack->GetSigmaSnp2(); + cmptrk.c30 = kaltrack->GetSigmaTglY(); + cmptrk.c31 = kaltrack->GetSigmaTglZ(); + cmptrk.c32 = kaltrack->GetSigmaTglSnp(); + cmptrk.c33 = kaltrack->GetSigmaTgl2(); + cmptrk.c40 = kaltrack->GetSigma1PtY(); + cmptrk.c41 = kaltrack->GetSigma1PtZ(); + cmptrk.c42 = kaltrack->GetSigma1PtSnp(); + cmptrk.c43 = kaltrack->GetSigma1PtTgl(); + cmptrk.c44 = kaltrack->GetSigma1Pt2(); + + // fill tree + cmptrktree->Fill(); + + } // loop on tracks at TPC 1st hit - if(gRandom->Rndm() < eff) return kTRUE; + delete [] kalLab; + //delete [] hasSeed; - return kFALSE; + } // end loop on events in file + + + cerr<<"+++\n+++ Number of compared tracks: "<Write(); + outfile->Close(); + + + // Write efficiencies to ascii file + FILE *effFile = fopen(tpceffasciiName,"w"); + //fprintf(effFile,"%d\n",kalEntries); + for(Int_t j=0; j=100) effPi[j]=(Float_t)kalPi[j]/geaPi[j]; + if(geaKa[j]>=100) effKa[j]=(Float_t)kalKa[j]/geaKa[j]; + if(geaPr[j]>=100) effPr[j]=(Float_t)kalPr[j]/geaPr[j]; + if(geaEl[j]>=500) effEl[j]=(Float_t)kalEl[j]/geaEl[j]; + if(geaMu[j]>=100) effMu[j]=(Float_t)kalMu[j]/geaMu[j]; + fprintf(effFile,"%f %f %f %f %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]); + } + + for(Int_t j=0; jClose(); + geaFile->Close(); + galiceFile->Close(); + + delete [] pt; + delete [] effPi; + delete [] effKa; + delete [] effPr; + delete [] effEl; + delete [] effMu; + delete [] geaPi; + delete [] geaKa; + delete [] geaPr; + delete [] geaEl; + delete [] geaMu; + delete [] kalPi; + delete [] kalKa; + delete [] kalPr; + delete [] kalEl; + delete [] kalMu; + + return; } +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) { +//----------------------------------------------------------------------------- +// This function assigns the track a dE/dx and makes a rough PID +//----------------------------------------------------------------------------- -//----------------------------------------------------------------- -Double_t AliTPCtrackerParam::LinearInterpolation(Int_t ptBins,Double_t *value, - Double_t trkPt,Double_t trkEta) const -{ -//----------------------------------------------------------------- -// This function makes a linear interpolation -//----------------------------------------------------------------- - Double_t intValue=0,intValue1=0,intValue2=0; - Int_t etaSide = (TMath::Abs(trkEta)<.45 ? 0 : 1); - Double_t eta[3]={0.15,0.45,0.75}; - Double_t pt[9]={0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.}; - if(ptBins==6) pt[5]=10.; - - for(Int_t i=0; iGetValueAt(pt,eta); + Double_t rms = fdEdxRMS->GetValueAt(pt,eta); + + Double_t dEdx = gRandom->Gaus(mean,rms); + + fTrack.SetdEdx(dEdx); + + AliTPCtrackParam t(fTrack); + + //Very rough PID + Double_t p = TMath::Sqrt(1.+t.GetTgl()*t.GetTgl())*pt; + + Double_t massPi = (Double_t)TDatabasePDG::Instance()->GetParticle(211)->Mass(); + Double_t massKa = (Double_t)TDatabasePDG::Instance()->GetParticle(321)->Mass(); + Double_t massPr = (Double_t)TDatabasePDG::Instance()->GetParticle(2212)->Mass(); + + if (p<0.6) { + if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; + } + if (dEdx < 39.+ 12./p/p) { + t.AssignMass(massKa); new(&fTrack) AliTPCtrack(t); return; } - if(i==ptBins-1) { - i=ptBins-2; - intValue1 = value[3*i-3+etaSide]+(value[3*i-3+etaSide]-value[3*i+etaSide])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]); - intValue2 = value[3*i-3+etaSide+1]+(value[3*i-3+etaSide+1]-value[3*i+etaSide+1])/(pt[i-1]-pt[i])*(trkPt-pt[i-1]); + t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return; + } - intValue = intValue1+(intValue1-intValue2)/(eta[etaSide]-eta[etaSide+1])*(trkEta-eta[etaSide]); - break; + if (p<1.2) { + if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; } + t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return; } - return intValue; -} - -//----------------------------------------------------------------- -Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const { -//----------------------------------------------------------------- -// This function tells bin number in a grid (pT,eta) -//----------------------------------------------------------------- - if(TMath::Abs(eta)<0.3) { - if(pt<0.3) return 0; - if(pt>=0.3 && pt<0.5) return 3; - if(pt>=0.5 && pt<1.) return 6; - if(pt>=1. && pt<1.5) return 9; - if(pt>=1.5 && pt<3.) return 12; - if(pt>=3. && pt<5.) return 15; - if(pt>=5. && pt<7.) return 18; - if(pt>=7. && pt<15.) return 21; - if(pt>=15.) return 24; - } - if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) { - if(pt<0.3) return 1; - if(pt>=0.3 && pt<0.5) return 4; - if(pt>=0.5 && pt<1.) return 7; - if(pt>=1. && pt<1.5) return 10; - if(pt>=1.5 && pt<3.) return 13; - if(pt>=3. && pt<5.) return 16; - if(pt>=5. && pt<7.) return 19; - if(pt>=7. && pt<15.) return 22; - if(pt>=15.) return 25; - } - if(TMath::Abs(eta)>=0.6) { - if(pt<0.3) return 2; - if(pt>=0.3 && pt<0.5) return 5; - if(pt>=0.5 && pt<1.) return 8; - if(pt>=1. && pt<1.5) return 11; - if(pt>=1.5 && pt<3.) return 14; - if(pt>=3. && pt<5.) return 17; - if(pt>=5. && pt<7.) return 20; - if(pt>=7. && pt<15.) return 23; - if(pt>=15.) return 26; - } - - return -1; - -} - -//----------------------------------------------------------------- -TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t eta) const -{ -//----------------------------------------------------------------- + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) { +//----------------------------------------------------------------------------- +// This function deals with covariance matrix and smearing +//----------------------------------------------------------------------------- + + COVMATRIX covmat; + Double_t p,cosl; + Double_t trkKine[1],trkRegPar[3]; + Double_t xref,alpha,xx[5],xxsm[5],cc[15]; + Int_t treeEntries; + + fCovTree->SetBranchAddress("matrix",&covmat); + + // get random entry from the tree + treeEntries = (Int_t)fCovTree->GetEntries(); + fCovTree->GetEvent(gRandom->Integer(treeEntries)); + + // get P and Cosl from track + cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl())); + p = 1./TMath::Abs(fTrack.Get1Pt())/cosl; + + trkKine[0] = p; + + // get covariance matrix from regularized matrix + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(0,l); + cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar); + cc[1] = covmat.c10; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(1,l); + cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar); + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(2,l); + cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar); + cc[4] = covmat.c21; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(3,l); + cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar); + cc[6] = covmat.c30; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(4,l); + cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar); + cc[8] = covmat.c32; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(5,l); + cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar); + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(6,l); + cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar); + cc[11]= covmat.c41; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(7,l); + cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar); + cc[13]= covmat.c43; + for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l); + cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar); + + + TMatrixD covMatSmear(5,5); + covMatSmear = GetSmearingMatrix(cc,pt,eta); + + // get track original parameters + xref =fTrack.GetX(); + alpha=fTrack.GetAlpha(); + xx[0]=fTrack.GetY(); + xx[1]=fTrack.GetZ(); + xx[2]=fTrack.GetSnp(); + xx[3]=fTrack.GetTgl(); + xx[4]=fTrack.GetC(); + + // use smearing matrix to smear the original parameters + xxsm[0]=xref; + SmearTrack(xx,xxsm,covMatSmear); + + AliTPCtrack track(xref,alpha,xxsm,cc,0); + new(&fTrack) AliTPCtrack(track); + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::DrawEffs(const Char_t* inName,Int_t pdg) { +//----------------------------------------------------------------------------- +// This function draws the TPC efficiencies in the [pT,eta] bins +//----------------------------------------------------------------------------- + + ReadEffs(inName); + SetParticle(pdg); + + const Int_t kn = fEff->GetPointsPt(); + Double_t *effsA = new Double_t[kn]; + Double_t *effsB = new Double_t[kn]; + Double_t *effsC = new Double_t[kn]; + Double_t *pt = new Double_t[kn]; + + fEff->GetArrayPt(pt); + for(Int_t i=0;iGetParam(i,0); + effsB[i] = fEff->GetParam(i,1); + effsC[i] = fEff->GetParam(i,2); + } + + TGraph *grA = new TGraph(kn,pt,effsA); + TGraph *grB = new TGraph(kn,pt,effsB); + TGraph *grC = new TGraph(kn,pt,effsC); + + grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5); + grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5); + grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5); + + TString title("Distribution of the TPC efficiencies"); + switch (pdg) { + case 211: + title.Prepend("PIONS - "); + break; + case 321: + title.Prepend("KAONS - "); + break; + case 2212: + title.Prepend("PROTONS - "); + break; + case 11: + title.Prepend("ELECTRONS - "); + break; + case 13: + title.Prepend("MUONS - "); + break; + } + + + gStyle->SetOptStat(0); + TCanvas *c = new TCanvas("c","effs",0,0,900,700); + c->SetLogx(); + c->SetGridx(); + c->SetGridy(); + + TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,1); + frame->SetXTitle("p_{T} [GeV/c]"); + frame->Draw(); + grA->Draw("P"); + grB->Draw("P"); + grC->Draw("P"); + + TLegend *leg = new TLegend(0.2,0.2,0.4,0.4); + leg->AddEntry(grA,"|#eta|<0.3","p"); + leg->AddEntry(grB,"0.3<|#eta|<0.6","p"); + leg->AddEntry(grC,"0.6<|#eta|<0.9","p"); + leg->SetFillColor(0); + leg->Draw(); + + delete [] effsA; + delete [] effsB; + delete [] effsC; + delete [] pt; + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::DrawPulls(const Char_t* inName,Int_t pdg, + Int_t par) { +//----------------------------------------------------------------------------- +// This function draws the pulls in the [pT,eta] bins +//----------------------------------------------------------------------------- + + ReadPulls(inName); + SetParticle(pdg); + + const Int_t kn = (fPulls+par)->GetPointsPt(); + Double_t *pullsA = new Double_t[kn]; + Double_t *pullsB = new Double_t[kn]; + Double_t *pullsC = new Double_t[kn]; + Double_t *pt = new Double_t[kn]; + (fPulls+par)->GetArrayPt(pt); + for(Int_t i=0;iGetParam(i,0); + pullsB[i] = (fPulls+par)->GetParam(i,1); + pullsC[i] = (fPulls+par)->GetParam(i,2); + } + + TGraph *grA = new TGraph(kn,pt,pullsA); + TGraph *grB = new TGraph(kn,pt,pullsB); + TGraph *grC = new TGraph(kn,pt,pullsC); + + grA->SetMarkerStyle(20); grA->SetMarkerColor(2); grA->SetMarkerSize(1.5); + grB->SetMarkerStyle(21); grB->SetMarkerColor(3); grB->SetMarkerSize(1.5); + grC->SetMarkerStyle(22); grC->SetMarkerColor(4); grC->SetMarkerSize(1.5); + + TString title("Distribution of the pulls: "); + switch (pdg) { + case 211: + title.Prepend("PIONS - "); + break; + case 321: + title.Prepend("KAONS - "); + break; + case 2212: + title.Prepend("PROTONS - "); + break; + case 11: + title.Prepend("ELECTRONS - "); + break; + case 13: + title.Prepend("MUONS - "); + break; + } + switch (par) { + case 0: + title.Append("y"); + break; + case 1: + title.Append("z"); + break; + case 2: + title.Append(" #eta"); + break; + case 3: + title.Append("tg #lambda"); + break; + case 4: + title.Append("C"); + break; + } + + gStyle->SetOptStat(0); + TCanvas *c = new TCanvas("c","pulls",0,0,900,700); + c->SetLogx(); + c->SetGridx(); + c->SetGridy(); + + TH2F *frame = new TH2F("frame",title.Data(),2,0.1,30,2,0,2); + frame->SetXTitle("p_{T} [GeV/c]"); + frame->Draw(); + grA->Draw("P"); + grB->Draw("P"); + grC->Draw("P"); + + TLegend *leg = new TLegend(0.2,0.2,0.4,0.4); + leg->AddEntry(grA,"|#eta|<0.3","p"); + leg->AddEntry(grB,"0.3<|#eta|<0.6","p"); + leg->AddEntry(grC,"0.6<|#eta|<0.9","p"); + leg->SetFillColor(0); + leg->Draw(); + + delete [] pullsA; + delete [] pullsB; + delete [] pullsC; + delete [] pt; + + return; +} +//----------------------------------------------------------------------------- +TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt, + Double_t eta) const { +//----------------------------------------------------------------------------- // This function stretches the covariance matrix according to the pulls -//----------------------------------------------------------------- +//----------------------------------------------------------------------------- + TMatrixD covMat(5,5); covMat(0,0)=cc[0]; @@ -420,6 +1706,7 @@ TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3); covMat(4,4)=cc[14]; + TMatrixD stretchMat(5,5); for(Int_t k=0;k<5;k++) { for(Int_t l=0;l<5;l++) { @@ -427,139 +1714,1969 @@ TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,Double_t } } - - Double_t s00[27]={1.3553,1.2973,1.2446,1.3428,1.3031,1.2345,1.3244,1.2681,1.2046,1.3046,1.2430,1.1606,1.2462,1.2104,1.1082,1.2207,1.1189,1.0789,1.2079,1.1300,1.0426,1.1502,1.1122,1.0559,1.1419,1.1072,1.0208}; - Double_t s11[27]={1.0890,1.1463,1.2313,1.1149,1.1597,1.2280,1.1225,1.1584,1.2329,1.1149,1.1550,1.2369,1.1095,1.1353,1.2050,1.0649,1.0858,1.1546,1.0663,1.0672,1.1340,1.0416,1.0738,1.0945,1.0663,1.0654,1.0909}; - Double_t s22[27]={1.1709,1.1367,1.0932,1.2247,1.1832,1.1247,1.2470,1.2017,1.1441,1.2202,1.1653,1.1050,1.1819,1.1269,1.0583,1.1546,1.0621,0.9928,1.1305,1.0512,0.9576,1.0995,1.0445,0.9884,1.0968,1.0368,0.9574}; - Double_t s33[27]={0.9720,0.9869,1.0271,1.0030,1.0223,1.0479,1.0164,1.0305,1.0559,1.0339,1.0450,1.0686,1.0450,1.0507,1.0784,1.0334,1.0208,1.0863,1.0252,1.0114,1.0835,0.9854,1.0144,1.0507,1.0124,1.0159,1.0464}; - Double_t s44[27]={1.1104,1.0789,1.0479,1.1597,1.1234,1.0728,1.2087,1.1602,1.1041,1.1942,1.1428,1.0831,1.1572,1.1036,1.0431,1.1296,1.0498,0.9844,1.1145,1.0266,0.9489,1.0959,1.0450,0.9875,1.0775,1.0266,0.9406}; - - - stretchMat(0,0) = LinearInterpolation(9,s00,pt,eta); - stretchMat(1,1) = LinearInterpolation(9,s11,pt,eta); - stretchMat(2,2) = LinearInterpolation(9,s22,pt,eta); - stretchMat(3,3) = LinearInterpolation(9,s33,pt,eta); - stretchMat(4,4) = LinearInterpolation(9,s44,pt,eta); + for(Int_t i=0;i<5;i++) { + stretchMat(i,i) = (fPulls+i)->GetValueAt(pt,eta); + if(stretchMat(i,i)==0.) stretchMat(i,i) = 1.; + } TMatrixD mat(stretchMat,TMatrixD::kMult,covMat); TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat); - return covMatSmear; } +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::InitializeKineGrid(Option_t* which) { +//----------------------------------------------------------------------------- +// This function initializes ([pt,eta] points) the data members AliTPCkineGrid +// which = "DB" -> initialize fDBgrid... members +// "eff" -> initialize fEff... members +// "pulls" -> initialize fPulls... members +// "dEdx" -> initialize fdEdx... members +//----------------------------------------------------------------------------- + + const char *db = strstr(which,"DB"); + const char *eff = strstr(which,"eff"); + const char *pulls = strstr(which,"pulls"); + const char *dEdx = strstr(which,"dEdx"); + + + Int_t nEta=0, nPt=0; + + Double_t etaPoints[2] = {0.3,0.6}; + Double_t etaBins[3] = {0.15,0.45,0.75}; + + Double_t ptPoints[9] = {0.4,0.6,0.8,1.2,1.7,3.,5.,8.,15.}; + Double_t ptBins[10] = {0.3,0.5,0.7,1.,1.5,2.,4.,6.,10.,20.}; + + + Double_t *eta=0,*pt=0; + + if(db) { + nEta = 2; + nPt = 9; + eta = etaPoints; + pt = ptPoints; + } else { + nEta = 3; + nPt = 10; + eta = etaBins; + pt = ptBins; + } -//----------------------------------------------------------------- -void AliTPCtrackerParam::SmearTrack(Double_t* xx,Double_t* xxsm,TMatrixD cov) const { -//----------------------------------------------------------------- -// This function smears track parameters according to streched cov. matrix -//----------------------------------------------------------------- - AliGausCorr *corgen = new AliGausCorr(cov,5); - TArrayD corr(5); - corgen->GetGaussN(corr); - delete corgen; - corgen = 0; + AliTPCkineGrid *dummy=0; - for(Int_t l=0;l<5;l++) { - xxsm[l] = xx[l]+corr[l]; + if(db) { + dummy = new AliTPCkineGrid(nPt,nEta,pt,eta); + new(&fDBgridPi) AliTPCkineGrid(*dummy); + new(&fDBgridKa) AliTPCkineGrid(*dummy); + new(&fDBgridPr) AliTPCkineGrid(*dummy); + new(&fDBgridEl) AliTPCkineGrid(*dummy); + new(&fDBgridMu) AliTPCkineGrid(*dummy); + delete dummy; + } + if(eff) { + dummy = new AliTPCkineGrid(nPt,nEta,pt,eta); + new(&fEffPi) AliTPCkineGrid(*dummy); + new(&fEffKa) AliTPCkineGrid(*dummy); + new(&fEffPr) AliTPCkineGrid(*dummy); + new(&fEffEl) AliTPCkineGrid(*dummy); + new(&fEffMu) AliTPCkineGrid(*dummy); + delete dummy; + } + if(pulls) { + dummy = new AliTPCkineGrid(nPt,nEta,pt,eta); + for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy); + for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy); + for(Int_t i=0;i<5;i++) new(&fPullsPr[i]) AliTPCkineGrid(*dummy); + for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy); + for(Int_t i=0;i<5;i++) new(&fPullsMu[i]) AliTPCkineGrid(*dummy); + delete dummy; + } + if(dEdx) { + dummy = new AliTPCkineGrid(nPt,nEta,pt,eta); + new(&fdEdxMeanPi) AliTPCkineGrid(*dummy); + new(&fdEdxRMSPi) AliTPCkineGrid(*dummy); + new(&fdEdxMeanKa) AliTPCkineGrid(*dummy); + new(&fdEdxRMSKa) AliTPCkineGrid(*dummy); + new(&fdEdxMeanPr) AliTPCkineGrid(*dummy); + new(&fdEdxRMSPr) AliTPCkineGrid(*dummy); + new(&fdEdxMeanEl) AliTPCkineGrid(*dummy); + new(&fdEdxRMSEl) AliTPCkineGrid(*dummy); + new(&fdEdxMeanMu) AliTPCkineGrid(*dummy); + new(&fdEdxRMSMu) AliTPCkineGrid(*dummy); + delete dummy; } return; } +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::MakeDataBase() { +//----------------------------------------------------------------------------- +// This function creates the DB file and store in it: +// - TPC Efficiencies for pi,ka,pr,el,mu (AliTPCkineGrid class) +// - Pulls for pi,ka,el (AliTPCkineGrid class) +// - Regularization parameters for pi,ka,el (TMatrixD class) +// - dE/dx parameterization for pi,ka,pr,el (AliTPCkineGrid class) +// - Regularized cov. matrices for pi,ka,el (COVMATRIX structure) +//----------------------------------------------------------------------------- + + // define some file names + const char *effFile ="TPCeff.root"; + const char *pullsFile ="pulls.root"; + const char *regPiFile ="regPi.root"; + const char *regKaFile ="regKa.root"; + const char *regPrFile ="regPr.root"; + const char *regElFile ="regEl.root"; + const char *regMuFile ="regMu.root"; + const char *dEdxPiFile="dEdxPi.root"; + const char *dEdxKaFile="dEdxKa.root"; + const char *dEdxPrFile="dEdxPr.root"; + const char *dEdxElFile="dEdxEl.root"; + const char *dEdxMuFile="dEdxMu.root"; + const char *cmFile ="CovMatrix_AllEvts.root"; + /* + Char_t *cmFile1 ="CovMatrix_AllEvts_1.root"; + Char_t *cmFile2 ="CovMatrix_AllEvts_2.root"; + Char_t *cmFile3 ="CovMatrix_AllEvts_3.root"; + */ + + // store the effieciencies + ReadEffs(effFile); + WriteEffs(fDBfileName.Data()); + + // store the pulls + ReadPulls(pullsFile); + WritePulls(fDBfileName.Data()); + + //store the regularization parameters + ReadRegParams(regPiFile,211); + WriteRegParams(fDBfileName.Data(),211); + ReadRegParams(regKaFile,321); + WriteRegParams(fDBfileName.Data(),321); + ReadRegParams(regPrFile,2212); + WriteRegParams(fDBfileName.Data(),2212); + ReadRegParams(regElFile,11); + WriteRegParams(fDBfileName.Data(),11); + ReadRegParams(regMuFile,13); + WriteRegParams(fDBfileName.Data(),13); + + // store the dEdx parameters + ReaddEdx(dEdxPiFile,211); + WritedEdx(fDBfileName.Data(),211); + ReaddEdx(dEdxKaFile,321); + WritedEdx(fDBfileName.Data(),321); + ReaddEdx(dEdxPrFile,2212); + WritedEdx(fDBfileName.Data(),2212); + ReaddEdx(dEdxElFile,11); + WritedEdx(fDBfileName.Data(),11); + ReaddEdx(dEdxMuFile,13); + WritedEdx(fDBfileName.Data(),13); + + + // + // store the regularized covariance matrices + // + InitializeKineGrid("DB"); + + const Int_t knBinsPi = fDBgridPi.GetTotBins(); + const Int_t knBinsKa = fDBgridKa.GetTotBins(); + const Int_t knBinsPr = fDBgridPr.GetTotBins(); + const Int_t knBinsEl = fDBgridEl.GetTotBins(); + const Int_t knBinsMu = fDBgridMu.GetTotBins(); + + + // create the trees for cov. matrices + // trees for pions + TTree *covTreePi1 = NULL; + covTreePi1 = new TTree[knBinsPi]; + // trees for kaons + TTree *covTreeKa1 = NULL; + covTreeKa1 = new TTree[knBinsKa]; + // trees for protons + TTree *covTreePr1 = NULL; + covTreePr1 = new TTree[knBinsPr]; + // trees for electrons + TTree *covTreeEl1 = NULL; + covTreeEl1 = new TTree[knBinsEl]; + // trees for muons + TTree *covTreeMu1 = NULL; + covTreeMu1 = new TTree[knBinsMu]; + + Char_t hname[100], htitle[100]; + COVMATRIX covmat; -//----------------------------------------------------------------- -void AliTPCtrackerParam::CookTracks(TObjArray& tarray,TObjArray& newtarray) -const { -//----------------------------------------------------------------- -// This function deals with covariance matrix and smearing -//----------------------------------------------------------------- + + for(Int_t i=0; iGet("cmptrktree"); + + COMPTRACK cmptrk; + cmptrktree->SetBranchAddress("comptracks",&cmptrk); + Int_t entries = (Int_t)cmptrktree->GetEntries(); + cerr<<" Number of entries: "<GetEvent(l); + // get the pdg code + trkPdg = TMath::Abs(cmptrk.pdg); + // use only pions, kaons, protons, electrons, muons + if(trkPdg!=211 && trkPdg!=321 && trkPdg!=2212 && trkPdg!=11 && trkPdg!=13) continue; + SetParticle(trkPdg); + trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta); + //cerr<=5000) continue; + if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue; + if(trkPdg==2212 && nPerBinPr[trkBin]>=5000) continue; + if(trkPdg==11 && nPerBinEl[trkBin]>=5000) continue; + if(trkPdg==13 && nPerBinMu[trkBin]>=5000) continue; + + trkKine[0] = cmptrk.p; + + // get regularized covariance matrix + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(0,k); + covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar); + covmat.c10 = cmptrk.c10; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(1,k); + covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar); + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(2,k); + covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar); + covmat.c21 = cmptrk.c21; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(3,k); + covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar); + covmat.c30 = cmptrk.c30; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(4,k); + covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar); + covmat.c32 = cmptrk.c32; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(5,k); + covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar); + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(6,k); + covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar); + covmat.c41 = cmptrk.c41; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(7,k); + covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar); + covmat.c43 = cmptrk.c43; + for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(8,k); + covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar); + + // fill the tree + switch (trkPdg) { + case 211: // pions + covTreePi1[trkBin].Fill(); + nPerBinPi[trkBin]++; + break; + case 321: // kaons + covTreeKa1[trkBin].Fill(); + nPerBinKa[trkBin]++; + break; + case 2212: // protons + covTreePr1[trkBin].Fill(); + nPerBinPr[trkBin]++; + break; + case 11: // electrons + covTreeEl1[trkBin].Fill(); + nPerBinEl[trkBin]++; + break; + case 13: // muons + covTreeMu1[trkBin].Fill(); + nPerBinMu[trkBin]++; + break; + } + } // loop on chain entries + + // store all trees the DB file + TFile *dbfile = new TFile(fDBfileName.Data(),"update"); + dbfile->mkdir("CovMatrices"); + gDirectory->cd("/CovMatrices"); + gDirectory->mkdir("Pions"); + gDirectory->mkdir("Kaons"); + gDirectory->mkdir("Protons"); + gDirectory->mkdir("Electrons"); + gDirectory->mkdir("Muons"); + // store pions + gDirectory->cd("/CovMatrices/Pions"); + fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write(); + for(Int_t i=0;icd("/CovMatrices/Kaons"); + fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write(); + for(Int_t i=0;icd("/CovMatrices/Protons"); + fDBgridPr.SetName("DBgridPr"); fDBgridPr.Write(); + for(Int_t i=0;icd("/CovMatrices/Electrons"); + fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write(); + for(Int_t i=0;icd("/CovMatrices/Muons"); + fDBgridMu.SetName("DBgridMu"); fDBgridMu.Write(); + for(Int_t i=0;iClose(); + delete [] nPerBinPi; + delete [] nPerBinKa; + delete [] nPerBinPr; + delete [] nPerBinEl; + delete [] nPerBinMu; + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th, + TObjArray &seedArray) const { +//----------------------------------------------------------------------------- +// This function makes the seeds for tracks from the 1st hits in the TPC +//----------------------------------------------------------------------------- + + Double_t xg,yg,zg,px,py,pz,pt; + Int_t label; + Int_t nTracks=(Int_t)th->GetEntries(); + + cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<ResetHits(); + th->GetEvent(l); + // Get FirstHit + tpcHit=(AliTPChit*)tpc->FirstHit(-1); + for( ; tpcHit; tpcHit=(AliTPChit*)tpc->NextHit() ) { + if(tpcHit->fQ !=0.) continue; + // Get particle momentum at hit + px=tpcHit->X(); py=tpcHit->Y(); pz=tpcHit->Z(); + + pt=TMath::Sqrt(px*px+py*py); + // reject hits with PtTrack(); + + if((tpcHit=(AliTPChit*)tpc->NextHit())==0) break; + if(tpcHit->fQ != 0.) continue; + // Get global coordinates of hit + xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z(); + if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue; + + // create the seed + AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label); + + // reject tracks which are not in the TPC acceptance + if(!ioseed->InTPCAcceptance()) { delete ioseed; continue; } + + // put seed in array + seedArray.AddLast(ioseed); + } + + } // loop over entries in TreeH + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr, + TObjArray &seedArray) const { +//----------------------------------------------------------------------------- +// This function makes the seeds for tracks from the track references +//----------------------------------------------------------------------------- + + Double_t xg,yg,zg,px,py,pz,pt; + Int_t label; + Int_t nnn,k; + + TClonesArray *tkRefArray = new TClonesArray("AliTrackReference"); + + TBranch *b =(TBranch*)ttr->GetBranch("TPC"); + if(!b) {cerr<<"TPC branch of TreeTR not found"<SetAddress(&tkRefArray); + Int_t nTkRef = (Int_t)b->GetEntries(); + cerr<<"+++ Number of entries in TreeTR(TPC): "<GetEvent(l)) continue; + nnn = tkRefArray->GetEntriesFast(); + + if(nnn <= 0) continue; + for(k=0; kUncheckedAt(k); + + xg = tkref->X(); + yg = tkref->Y(); + zg = tkref->Z(); + px = tkref->Px(); + py = tkref->Py(); + pz = tkref->Pz(); + label = tkref->GetTrack(); + + pt=TMath::Sqrt(px*px+py*py); + // reject hits with PtGetXL()-83.8) > 0.2) { + //cerr<<"not at TPC inner part: XL = "<GetXL()<InTPCAcceptance()) { + delete ioseed; continue; + } + + // put seed in array + seedArray.AddLast(ioseed); - TString s("$ALICE_ROOT/TPC/CovMatrixDB_"); - if(fColl==0) s.Append("PbPb6000"); - if(fBz==0.4) s.Append("_B0.4T.root"); - // open file with matrixes DB - TFile* fileDB = TFile::Open(s.Data()); + } + + + } // loop on track references + + delete tkRefArray; + + + return; +} +//----------------------------------------------------------------------------- +void AliTPCtrackerParam::MergeEvents(Int_t evFirst,Int_t evLast) { +//----------------------------------------------------------------------------- +// This function: 1) merges the files from track comparison +// (beware: better no more than 100 events per file) +// 2) computes the average TPC efficiencies +//----------------------------------------------------------------------------- + + const char *outName="TPCeff.root"; + + // merge files with tracks + cerr<<" ******** MERGING FILES **********\n\n"; + + // 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 "<AccessPathName(covName.Data(),kFileExists)) continue; + - Double_t pt = 1/TMath::Abs(track->Get1Pt()); - Double_t eta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(track->GetTgl()))); + if(j<=100) ch1.Add(covName.Data()); + if(j>100 && j<=200) ch2.Add(covName.Data()); + if(j>200) ch3.Add(covName.Data()); - Int_t bin = GetBin(pt,eta); + } // 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; jAccessPathName(effName.Data(),kFileExists)) continue; + + ReadEffs(effName.Data()); - // string with tree name - TString str("Tree_bin"); - str+=bin; + for(Int_t k=0; 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; jGet(str.Data()); - TArrayF* matrix = new TArrayF; - covTree->SetBranchAddress("matrixes",&matrix); - - // get random entry from the tree - Int_t treeEntries = (Int_t)covTree->GetEntries(); - covTree->GetEvent(gRandom->Integer(treeEntries)); - - Double_t xref,alpha,xx[5],xxsm[5],cc[15]; - Int_t lab; - // get P and Cosl from track - Double_t cosl = TMath::Cos(TMath::ATan(track->GetTgl())); - Double_t p = 1/TMath::Abs(track->Get1Pt())/cosl; + // 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; + if(!ReadCovTrees(inName)) return 0; + + return 1; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::ReadCovTrees(const Char_t* inName) { +//----------------------------------------------------------------------------- +// This function reads the covariance matrix trees from the DB +//----------------------------------------------------------------------------- + + if(gSystem->AccessPathName(inName,kFileExists)) { + cerr<<"AliTPCtrackerParam::ReaddEdx: "<Get(str.Data()); + //fCovTree = (TTree*)inFile->Get(str.Data()); + //fCovTreePi[l] = new TTree(*fCovTree); + } + Int_t nBinsKa = fDBgridKa.GetTotBins(); + for(Int_t l=0; lGet(str.Data()); + //fCovTree = (TTree*)inFile->Get(str.Data()); + //fCovTreeKa[l] = new TTree(*fCovTree); + } + Int_t nBinsPr = fDBgridPr.GetTotBins(); + for(Int_t l=0; lGet(str.Data()); + //fCovTree = (TTree*)inFile->Get(str.Data()); + //fCovTreePr[l] = new TTree(*fCovTree); + } + Int_t nBinsEl = fDBgridEl.GetTotBins(); + for(Int_t l=0; lGet(str.Data()); + //fCovTree = (TTree*)inFile->Get(str.Data()); + //fCovTreeEl[l] = new TTree(*fCovTree); + } + Int_t nBinsMu = fDBgridMu.GetTotBins(); + for(Int_t l=0; lGet(str.Data()); + //fCovTree = (TTree*)inFile->Get(str.Data()); + //fCovTreeMu[l] = new TTree(*fCovTree); + } + + //inFile->Close(); + + 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: "<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: "<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: "<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: "<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: "<Get("/RegParams/Pions/RegPions"); + new(&fRegParPi) TMatrixD(*fRegPar); + //printf("before\n"); + //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2)); + dummyMatrix = fRegParPi; + fRegParPi(0,0) = dummyMatrix(0,0); + fRegParPi(0,1) = dummyMatrix(0,1); + fRegParPi(0,2) = dummyMatrix(0,2); + fRegParPi(1,0) = dummyMatrix(3,0); + fRegParPi(1,1) = dummyMatrix(1,1); + fRegParPi(1,2) = dummyMatrix(1,2); + fRegParPi(2,0) = dummyMatrix(6,0); + fRegParPi(2,1) = dummyMatrix(3,2); + fRegParPi(2,2) = dummyMatrix(2,2); + fRegParPi(3,0) = dummyMatrix(1,0); + fRegParPi(3,1) = dummyMatrix(4,0); + fRegParPi(3,2) = dummyMatrix(7,0); + fRegParPi(4,0) = dummyMatrix(3,1); + fRegParPi(4,1) = dummyMatrix(4,1); + fRegParPi(4,2) = dummyMatrix(7,1); + fRegParPi(5,0) = dummyMatrix(6,1); + fRegParPi(5,1) = dummyMatrix(4,2); + fRegParPi(5,2) = dummyMatrix(7,2); + fRegParPi(6,0) = dummyMatrix(2,0); + fRegParPi(6,1) = dummyMatrix(5,0); + fRegParPi(6,2) = dummyMatrix(8,0); + fRegParPi(7,0) = dummyMatrix(2,1); + fRegParPi(7,1) = dummyMatrix(5,1); + fRegParPi(7,2) = dummyMatrix(8,1); + fRegParPi(8,0) = dummyMatrix(6,2); + fRegParPi(8,1) = dummyMatrix(5,2); + fRegParPi(8,2) = dummyMatrix(8,2); + //printf("after\n"); + //for(Int_t jj=0;jj<9;jj++) printf("%13.10f %13.10f %13.10f\n",fRegParPi(jj,0),fRegParPi(jj,1),fRegParPi(jj,2)); + break; + case 321: + fRegPar = (TMatrixD*)inFile->Get("/RegParams/Kaons/RegKaons"); + new(&fRegParKa) TMatrixD(*fRegPar); + dummyMatrix = fRegParKa; + fRegParKa(0,0) = dummyMatrix(0,0); + fRegParKa(0,1) = dummyMatrix(0,1); + fRegParKa(0,2) = dummyMatrix(0,2); + fRegParKa(1,0) = dummyMatrix(3,0); + fRegParKa(1,1) = dummyMatrix(1,1); + fRegParKa(1,2) = dummyMatrix(1,2); + fRegParKa(2,0) = dummyMatrix(6,0); + fRegParKa(2,1) = dummyMatrix(3,2); + fRegParKa(2,2) = dummyMatrix(2,2); + fRegParKa(3,0) = dummyMatrix(1,0); + fRegParKa(3,1) = dummyMatrix(4,0); + fRegParKa(3,2) = dummyMatrix(7,0); + fRegParKa(4,0) = dummyMatrix(3,1); + fRegParKa(4,1) = dummyMatrix(4,1); + fRegParKa(4,2) = dummyMatrix(7,1); + fRegParKa(5,0) = dummyMatrix(6,1); + fRegParKa(5,1) = dummyMatrix(4,2); + fRegParKa(5,2) = dummyMatrix(7,2); + fRegParKa(6,0) = dummyMatrix(2,0); + fRegParKa(6,1) = dummyMatrix(5,0); + fRegParKa(6,2) = dummyMatrix(8,0); + fRegParKa(7,0) = dummyMatrix(2,1); + fRegParKa(7,1) = dummyMatrix(5,1); + fRegParKa(7,2) = dummyMatrix(8,1); + fRegParKa(8,0) = dummyMatrix(6,2); + fRegParKa(8,1) = dummyMatrix(5,2); + fRegParKa(8,2) = dummyMatrix(8,2); + 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); + dummyMatrix = fRegParPr; + fRegParPr(0,0) = dummyMatrix(0,0); + fRegParPr(0,1) = dummyMatrix(0,1); + fRegParPr(0,2) = dummyMatrix(0,2); + fRegParPr(1,0) = dummyMatrix(3,0); + fRegParPr(1,1) = dummyMatrix(1,1); + fRegParPr(1,2) = dummyMatrix(1,2); + fRegParPr(2,0) = dummyMatrix(6,0); + fRegParPr(2,1) = dummyMatrix(3,2); + fRegParPr(2,2) = dummyMatrix(2,2); + fRegParPr(3,0) = dummyMatrix(1,0); + fRegParPr(3,1) = dummyMatrix(4,0); + fRegParPr(3,2) = dummyMatrix(7,0); + fRegParPr(4,0) = dummyMatrix(3,1); + fRegParPr(4,1) = dummyMatrix(4,1); + fRegParPr(4,2) = dummyMatrix(7,1); + fRegParPr(5,0) = dummyMatrix(6,1); + fRegParPr(5,1) = dummyMatrix(4,2); + fRegParPr(5,2) = dummyMatrix(7,2); + fRegParPr(6,0) = dummyMatrix(2,0); + fRegParPr(6,1) = dummyMatrix(5,0); + fRegParPr(6,2) = dummyMatrix(8,0); + fRegParPr(7,0) = dummyMatrix(2,1); + fRegParPr(7,1) = dummyMatrix(5,1); + fRegParPr(7,2) = dummyMatrix(8,1); + fRegParPr(8,0) = dummyMatrix(6,2); + fRegParPr(8,1) = dummyMatrix(5,2); + fRegParPr(8,2) = dummyMatrix(8,2); + break; + case 11: + fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons"); + new(&fRegParEl) TMatrixD(*fRegPar); + dummyMatrix = fRegParEl; + fRegParEl(0,0) = dummyMatrix(0,0); + fRegParEl(0,1) = dummyMatrix(0,1); + fRegParEl(0,2) = dummyMatrix(0,2); + fRegParEl(1,0) = dummyMatrix(3,0); + fRegParEl(1,1) = dummyMatrix(1,1); + fRegParEl(1,2) = dummyMatrix(1,2); + fRegParEl(2,0) = dummyMatrix(6,0); + fRegParEl(2,1) = dummyMatrix(3,2); + fRegParEl(2,2) = dummyMatrix(2,2); + fRegParEl(3,0) = dummyMatrix(1,0); + fRegParEl(3,1) = dummyMatrix(4,0); + fRegParEl(3,2) = dummyMatrix(7,0); + fRegParEl(4,0) = dummyMatrix(3,1); + fRegParEl(4,1) = dummyMatrix(4,1); + fRegParEl(4,2) = dummyMatrix(7,1); + fRegParEl(5,0) = dummyMatrix(6,1); + fRegParEl(5,1) = dummyMatrix(4,2); + fRegParEl(5,2) = dummyMatrix(7,2); + fRegParEl(6,0) = dummyMatrix(2,0); + fRegParEl(6,1) = dummyMatrix(5,0); + fRegParEl(6,2) = dummyMatrix(8,0); + fRegParEl(7,0) = dummyMatrix(2,1); + fRegParEl(7,1) = dummyMatrix(5,1); + fRegParEl(7,2) = dummyMatrix(8,1); + fRegParEl(8,0) = dummyMatrix(6,2); + fRegParEl(8,1) = dummyMatrix(5,2); + fRegParEl(8,2) = dummyMatrix(8,2); + 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); + dummyMatrix = fRegParMu; + fRegParMu(0,0) = dummyMatrix(0,0); + fRegParMu(0,1) = dummyMatrix(0,1); + fRegParMu(0,2) = dummyMatrix(0,2); + fRegParMu(1,0) = dummyMatrix(3,0); + fRegParMu(1,1) = dummyMatrix(1,1); + fRegParMu(1,2) = dummyMatrix(1,2); + fRegParMu(2,0) = dummyMatrix(6,0); + fRegParMu(2,1) = dummyMatrix(3,2); + fRegParMu(2,2) = dummyMatrix(2,2); + fRegParMu(3,0) = dummyMatrix(1,0); + fRegParMu(3,1) = dummyMatrix(4,0); + fRegParMu(3,2) = dummyMatrix(7,0); + fRegParMu(4,0) = dummyMatrix(3,1); + fRegParMu(4,1) = dummyMatrix(4,1); + fRegParMu(4,2) = dummyMatrix(7,1); + fRegParMu(5,0) = dummyMatrix(6,1); + fRegParMu(5,1) = dummyMatrix(4,2); + fRegParMu(5,2) = dummyMatrix(7,2); + fRegParMu(6,0) = dummyMatrix(2,0); + fRegParMu(6,1) = dummyMatrix(5,0); + fRegParMu(6,2) = dummyMatrix(8,0); + fRegParMu(7,0) = dummyMatrix(2,1); + fRegParMu(7,1) = dummyMatrix(5,1); + fRegParMu(7,2) = dummyMatrix(8,1); + fRegParMu(8,0) = dummyMatrix(6,2); + fRegParMu(8,1) = dummyMatrix(5,2); + fRegParMu(8,2) = dummyMatrix(8,2); + 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; + const char *part="Pions - "; + + InitializeKineGrid("DB"); + SetParticle(pdg); + const Int_t kfitbins = fDBgrid->GetBinsPt(); + cerr<<" Fit bins: "<Get("cmptrktree"); + + COMPTRACK cmptrk; + cmptrktree->SetBranchAddress("comptracks",&cmptrk); + Int_t entries = (Int_t)cmptrktree->GetEntries(); + - // get covariance matrix from regularized matrix - cc[0]=matrix->At(0)*(1.588e-3+1.911e-4/TMath::Power(p,1.5)); - cc[1]=matrix->At(1); - cc[2]=matrix->At(2)*(1.195e-3+0.8102e-3/p); - cc[3]=matrix->At(3)*(7.275e-5+1.181e-5/TMath::Power(p,1.5)); - cc[4]=matrix->At(4); - cc[5]=matrix->At(5)*(5.163e-6+1.138e-6/TMath::Power(p,2.)/cosl); - cc[6]=matrix->At(6); - cc[7]=matrix->At(7)*(1.176e-5+1.175e-5/TMath::Power(p,1.5)/cosl/cosl/cosl); - cc[8]=matrix->At(8); - cc[9]=matrix->At(9)*(1.289e-7+4.636e-7/TMath::Power(p,1.7)/cosl/cosl/cosl/cosl); - cc[10]=matrix->At(10)*(4.056e-7+4.379e-8/TMath::Power(p,1.5)); - cc[11]=matrix->At(11); - cc[12]=matrix->At(12)*(3.049e-8+8.244e-9/TMath::Power(p,2.)/cosl); - cc[13]=matrix->At(13); - cc[14]=matrix->At(14)*(1.847e-10+5.822e-11/TMath::Power(p,2.)/cosl/cosl/cosl); + Int_t pbin; + Int_t *n = new Int_t[kfitbins]; + Int_t *n00 = new Int_t[kfitbins]; + Int_t *n11 = new Int_t[kfitbins]; + Int_t *n20 = new Int_t[kfitbins]; + Int_t *n22 = new Int_t[kfitbins]; + Int_t *n31 = new Int_t[kfitbins]; + Int_t *n33 = new Int_t[kfitbins]; + Int_t *n40 = new Int_t[kfitbins]; + Int_t *n42 = new Int_t[kfitbins]; + Int_t *n44 = new Int_t[kfitbins]; + Double_t *p = new Double_t[kfitbins]; + Double_t *ep = new Double_t[kfitbins]; + Double_t *mean00 = new Double_t[kfitbins]; + Double_t *mean11 = new Double_t[kfitbins]; + Double_t *mean20 = new Double_t[kfitbins]; + Double_t *mean22 = new Double_t[kfitbins]; + Double_t *mean31 = new Double_t[kfitbins]; + Double_t *mean33 = new Double_t[kfitbins]; + Double_t *mean40 = new Double_t[kfitbins]; + Double_t *mean42 = new Double_t[kfitbins]; + Double_t *mean44 = new Double_t[kfitbins]; + Double_t *sigma00 = new Double_t[kfitbins]; + Double_t *sigma11 = new Double_t[kfitbins]; + Double_t *sigma20 = new Double_t[kfitbins]; + Double_t *sigma22 = new Double_t[kfitbins]; + Double_t *sigma31 = new Double_t[kfitbins]; + Double_t *sigma33 = new Double_t[kfitbins]; + Double_t *sigma40 = new Double_t[kfitbins]; + Double_t *sigma42 = new Double_t[kfitbins]; + Double_t *sigma44 = new Double_t[kfitbins]; + Double_t *rmean = new Double_t[kfitbins]; + Double_t *rsigma = new Double_t[kfitbins]; + Double_t fitpar[3]; + + for(Int_t l=0; lGetEvent(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; lGetEvent(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; lGetX(); - alpha=track->GetAlpha(); - xx[0]=track->GetY(); - xx[1]=track->GetZ(); - xx[2]=track->GetX()*track->GetC()-track->GetSnp(); - xx[3]=track->GetTgl(); - xx[4]=track->GetC(); - lab=track->GetLabel(); - - // use smearing matrix to smear the original parameters - SmearTrack(xx,xxsm,covMatSmear); - - AliTPCtrack* tpctrack = new AliTPCtrack(0,xxsm,cc,xref,alpha); - tpctrack->SetLabel(lab); - - // fill the array - newtarray.AddLast(tpctrack); + + // 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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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(kfitbins,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; lSetXTitle("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.; - delete matrix; + 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 +//----------------------------------------------------------------------------- + Double_t xref=xxsm[0]; xxsm[0]=0; + AliGausCorr *corgen = new AliGausCorr(cov,5); + TArrayD corr(5); + corgen->GetGaussN(corr); + // check on fP2(ESD) = fX*fP4-fP2(TPC) + // with fX=xref (not smeared), fP4=xx[4]+corr[4] e fP2=xx[2]+corr[2]; + // if fP2(ESD)^2 > 1 -> resmear... + Double_t fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]); + while ( (fp2esd*fp2esd) > 1.0 ) { + Warning("SmearTrack()","Badly smeared track, retrying..."); + corgen->GetGaussN(corr); + fp2esd=xref*(xx[4]+corr[4])-(xx[2]+corr[2]); } - fileDB->Close(); - delete fileDB; + delete corgen; + corgen = 0; + for(Int_t l=0;l<5;l++) xxsm[l] = xx[l]+corr[l]; - return; + return; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::WritedEdx(const Char_t *outName,Int_t pdg) { +//----------------------------------------------------------------------------- +// This function writes the dEdx parameters to the DB +//----------------------------------------------------------------------------- + + Option_t *opt; + const char *dirName="Pions"; + const char *meanName="dEdxMeanPi"; + const char *rmsName="dEdxRMSPi"; + + SetParticle(pdg); + + if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate"; + } else { opt="update"; } + + switch (pdg) { + case 211: + dirName="Pions"; + meanName="dEdxMeanPi"; + rmsName="dEdxRMSPi"; + break; + case 321: + dirName="Kaons"; + meanName="dEdxMeanKa"; + rmsName="dEdxRMSKa"; + break; + case 2212: + dirName="Protons"; + meanName="dEdxMeanPr"; + rmsName="dEdxRMSPr"; + break; + case 11: + dirName="Electrons"; + meanName="dEdxMeanEl"; + rmsName="dEdxRMSEl"; + break; + case 13: + dirName="Muons"; + meanName="dEdxMeanMu"; + rmsName="dEdxRMSMu"; + break; + } + + TFile *outFile = new TFile(outName,opt); + if(!gDirectory->cd("/dEdx")) { + outFile->mkdir("dEdx"); + gDirectory->cd("/dEdx"); + } + TDirectory *dir2 = gDirectory->mkdir(dirName); + dir2->cd(); + fdEdxMean->SetName(meanName); fdEdxMean->Write(); + fdEdxRMS->SetName(rmsName); fdEdxRMS->Write(); + + outFile->Close(); + delete outFile; + + + return 1; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::WriteEffs(const Char_t *outName) { +//----------------------------------------------------------------------------- +// This function writes the TPC efficiencies to the DB +//----------------------------------------------------------------------------- + + + Option_t *opt; + + if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate"; + } else { opt="update"; } + + TFile *outFile = new TFile(outName,opt); + + outFile->mkdir("Efficiencies"); + gDirectory->cd("/Efficiencies"); + gDirectory->mkdir("Pions"); + gDirectory->mkdir("Kaons"); + gDirectory->mkdir("Protons"); + gDirectory->mkdir("Electrons"); + gDirectory->mkdir("Muons"); + + gDirectory->cd("/Efficiencies/Pions"); + fEffPi.SetName("EffPi"); + fEffPi.Write(); + gDirectory->cd("/Efficiencies/Kaons"); + fEffKa.SetName("EffKa"); + fEffKa.Write(); + gDirectory->cd("/Efficiencies/Protons"); + fEffPr.SetName("EffPr"); + fEffPr.Write(); + gDirectory->cd("/Efficiencies/Electrons"); + fEffEl.SetName("EffEl"); + fEffEl.Write(); + gDirectory->cd("/Efficiencies/Muons"); + fEffMu.SetName("EffMu"); + fEffMu.Write(); + + outFile->Close(); + + delete outFile; + + return 1; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::WritePulls(const Char_t *outName) { +//----------------------------------------------------------------------------- +// This function writes the pulls to the DB +//----------------------------------------------------------------------------- + + Option_t *opt; + + if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate"; + } else { opt="update"; } + + TFile *outFile = new TFile(outName,opt); + + outFile->mkdir("Pulls"); + gDirectory->cd("/Pulls"); + gDirectory->mkdir("Pions"); + gDirectory->mkdir("Kaons"); + gDirectory->mkdir("Protons"); + gDirectory->mkdir("Electrons"); + gDirectory->mkdir("Muons"); + + for(Int_t i=0;i<5;i++) { + TString pi("PullsPi_"); pi+=i; + TString ka("PullsKa_"); ka+=i; + TString pr("PullsPr_"); pr+=i; + TString el("PullsEl_"); el+=i; + TString mu("PullsMu_"); mu+=i; + fPullsPi[i].SetName(pi.Data()); + fPullsKa[i].SetName(ka.Data()); + fPullsPr[i].SetName(pr.Data()); + fPullsEl[i].SetName(el.Data()); + fPullsMu[i].SetName(mu.Data()); + gDirectory->cd("/Pulls/Pions"); + fPullsPi[i].Write(); + gDirectory->cd("/Pulls/Kaons"); + fPullsKa[i].Write(); + gDirectory->cd("/Pulls/Protons"); + fPullsPr[i].Write(); + gDirectory->cd("/Pulls/Electrons"); + fPullsEl[i].Write(); + gDirectory->cd("/Pulls/Muons"); + fPullsMu[i].Write(); + } + outFile->Close(); + delete outFile; + + return 1; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::WriteRegParams(const Char_t *outName,Int_t pdg) { +//----------------------------------------------------------------------------- +// This function writes the regularization parameters to the DB +//----------------------------------------------------------------------------- + + Option_t *opt; + const char *dirName="Pions"; + const char *keyName="RegPions"; + + SetParticle(pdg); + + if(gSystem->AccessPathName(outName,kFileExists)) { opt="recreate"; + } else { opt="update"; } + + switch (pdg) { + case 211: + dirName="Pions"; + keyName="RegPions"; + break; + case 321: + dirName="Kaons"; + keyName="RegKaons"; + break; + case 2212: + dirName="Protons"; + keyName="RegProtons"; + break; + case 11: + dirName="Electrons"; + keyName="RegElectrons"; + break; + case 13: + dirName="Muons"; + keyName="RegMuons"; + break; + } + + TFile *outFile = new TFile(outName,opt); + if(!gDirectory->cd("/RegParams")) { + outFile->mkdir("RegParams"); + gDirectory->cd("/RegParams"); + } + TDirectory *dir2 = gDirectory->mkdir(dirName); + dir2->cd(); + fRegPar->Write(keyName); + + outFile->Close(); + delete outFile; + + + return 1; } -//----------------------------------------------------------------- +//----------------------------------------------------------------------------- +//***************************************************************************** +//----------------------------------------------------------------------------- + + + + + + + + + + +