X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCtrackerParam.cxx;h=60983c3e8b8d95ebadda013e92694f8c1b65a1ff;hb=cc65e4f59e4fc079ebe8407b30c207571da26430;hp=c64d9533885a9b26df11ba44a79f187993051fd9;hpb=8e8eae84c82d2951bab0d875cdd742d409f5b196;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCtrackerParam.cxx b/TPC/AliTPCtrackerParam.cxx index c64d9533885..60983c3e8b8 100644 --- a/TPC/AliTPCtrackerParam.cxx +++ b/TPC/AliTPCtrackerParam.cxx @@ -51,8 +51,10 @@ * - All the code necessary to create a BataBase has been included in * * class (see the macro AliTPCtrackingParamDB.C for the details). * * * - * Origin: Andrea Dainese, Padova - e-mail: andrea.dainese@pd.infn.it * + * 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 * **************************************************************************/ // * // This is a dummy comment @@ -60,29 +62,36 @@ // // * //------- Root headers -------- +#include +#include #include #include #include #include +#include #include #include #include -#include +#include +#include +#include #include -#include +#include +#include //------ AliRoot headers ------ -#include "alles.h" #include "AliGausCorr.h" -#include "AliKalmanTrack.h" +#include "AliTracker.h" +#include "AliMC.h" #include "AliMagF.h" -#include "AliMagFCM.h" +#include "AliRun.h" #include "AliRunLoader.h" -#include "AliTPCLoader.h" +#include "AliTPC.h" +#include "AliTPCParamSR.h" #include "AliTPCkineGrid.h" #include "AliTPCtrack.h" -#include "AliTrackReference.h" #include "AliTPCtrackerParam.h" -#include "AliMC.h" +#include "AliTrackReference.h" +#include "AliESDtrack.h" //----------------------------- Double_t RegFunc(Double_t *x,Double_t *par) { @@ -114,25 +123,59 @@ ClassImp(AliTPCtrackerParam) //----------------------------------------------------------------------------- AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz, - Int_t kn, const char* evfoldname): - fEvFolderName(evfoldname) { + 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 //----------------------------------------------------------------------------- - fNevents = kn; // events to be processed - 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 + // 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(fBz!=0.4) { - cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n"; - cerr<<" Available: 0.4"< PbPb6000"< pp"< PbPb6000\n 1 -> pp"); } fDBfileName = gSystem->Getenv("ALICE_ROOT"); @@ -141,11 +184,52 @@ AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz, 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) +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 } @@ -153,17 +237,23 @@ AliTPCtrackerParam::AliTPCtrackerParam( const AliTPCtrackerParam& p):TObject(p) AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant( Double_t x,Double_t y,Double_t z, Double_t px,Double_t py,Double_t pz, - Int_t lab) { + 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 //---------------------------------------------------------------------------- - fXg = x; - fYg = y; - fZg = z; - fPx = px; - fPy = py; - fPz = pz; - fLabel = lab; + Double_t a = TMath::ATan2(y,x)*180./TMath::Pi(); if(a<0) a += 360.; fSector = (Int_t)(a/20.); @@ -172,307 +262,239 @@ AliTPCtrackerParam::AliTPCseedGeant::AliTPCseedGeant( fAlpha *= TMath::Pi(); } //----------------------------------------------------------------------------- -Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out) { +Int_t AliTPCtrackerParam::Init() { //----------------------------------------------------------------------------- -// This function creates the TPC parameterized tracks +// This function reads the DB from file //----------------------------------------------------------------------------- - Error("BuildTPCtracks","in and out parameters ignored. new io"); - -/********************************************/ - 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; - } - -/********************************************/ - - TFile *fileDB=0; - TTree *covTreePi[50]; - TTree *covTreeKa[50]; - TTree *covTreePr[50]; - TTree *covTreeEl[50]; - TTree *covTreeMu[50]; - if(fSelAndSmear) { - cerr<<"+++\n+++ Reading DataBase from:\n+++ "<< - fDBfileName.Data()<<"\n+++\n"; + printf("+++\n+++ Reading DataBase from:%s\n+++\n+++\n",fDBfileName.Data()); // Read paramters from DB file if(!ReadAllData(fDBfileName.Data())) { - cerr<<"AliTPCtrackerParam::BuildTPCtracks: \ - Could not read data from DB\n\n"; return 1; + 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"); - // Read the trees with regularized cov. matrices from DB - TString str; - fileDB = TFile::Open(fDBfileName.Data()); - Int_t nBinsPi = fDBgridPi.GetTotBins(); - for(Int_t l=0; lGet(str.Data()); - } - Int_t nBinsKa = fDBgridKa.GetTotBins(); - for(Int_t l=0; lGet(str.Data()); - } - Int_t nBinsPr = fDBgridPr.GetTotBins(); - for(Int_t l=0; lGet(str.Data()); - } - Int_t nBinsEl = fDBgridEl.GetTotBins(); - for(Int_t l=0; lGet(str.Data()); - } - Int_t nBinsMu = fDBgridMu.GetTotBins(); - for(Int_t l=0; lGet(str.Data()); - } + AliMagF *fiel = (AliMagF*)TGeoGlobalMagField::Instance()->GetField(); + Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.); + printf("Magnetic field is %6.2f Tesla\n",fieval); + if(fBz!=fieval) { + printf("AliTPCtrackerParam::BuildTPCtracks: Invalid field!"); + printf("Field selected is: %f T\n",fBz); + printf("Field found on file is: %f T\n",fieval); + return 1; + } - } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n"; + return 0; +} +//----------------------------------------------------------------------------- +Int_t AliTPCtrackerParam::BuildTPCtracks(AliESDEvent *event) { +//----------------------------------------------------------------------------- +// This function creates the TPC parameterized tracks and writes them +// as AliESDtrack objects in the ESD event +//----------------------------------------------------------------------------- - TFile *infile=(TFile*)inp; - infile->cd(); - // Get gAlice object from file + if(!event) return -1; + + 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; + } - rl->LoadgAlice(); - rl->LoadHeader(); + // Get gAlice object from file + if(!gAlice) rl->LoadgAlice(); + //rl->LoadHeader(); + rl->LoadKinematics(); tpcloader->LoadHits("read"); if(!(gAlice=rl->GetAliRun())) { - cerr<<"Can not get gAlice from Run Loader !\n"; + printf("Cannot get gAlice from Run Loader !\n"); return 1; } - // Check if value in the galice file is equal to selected one (fBz) - AliMagF *fiel = (AliMagF*)gAlice->Field(); - Double_t fieval=(Double_t)fiel->SolenoidField()/10.; - printf("Magnetic field is %6.2f Tesla\n",fieval); - if(fBz!=fieval) { - cerr<<"AliTPCtrackerParam::BuildTPCtracks: Invalid field!"<GetDetector("TPC"); - Int_t ver = tpc->IsVersion(); - cerr<<"+++ TPC version "<CdGAFile(); - AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60"); + AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60"); if(digp){ delete digp; digp = new AliTPCParamSR(); } - else digp=(AliTPCParam*)infile->Get("75x40_100x60_150x60"); + else digp=(AliTPCParam*)gDirectory->Get("75x40_100x60_150x60"); if(!digp) { cerr<<"TPC parameters have not been found !\n"; return 1; } tpc->SetParam(digp); - // Set the conversion constant between curvature and Pt - AliKalmanTrack::SetConvConst(100/0.299792458/fBz); - TParticle *part=0; AliTPCseedGeant *seed=0; AliTPCtrack *tpctrack=0; Double_t sPt,sEta; - Int_t cl=0,bin,label,pdg,charge; + Int_t bin,label,pdg,charge; Int_t tracks; Int_t nParticles,nSeeds,arrentr; - Char_t tname[100]; //Int_t nSel=0,nAcc=0; + Int_t evt=event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number. + + tracks=0; - // loop over first n events in file - for(Int_t evt=0; evtBranch("tracks","AliTPCtrack",&tpctrack,20000,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); + // 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]; - Double_t *ptkine = new Double_t[nParticles]; - Double_t *pzkine = new Double_t[nParticles]; - - // loop on particles and store pdg codes - for(Int_t l=0; lGetMCApp()->Particle(l); - pdgCodes[l] = part->GetPdgCode(); - ptkine[l] = part->Pt(); - pzkine[l] = part->Pz(); - done[l] = kFALSE; - } - cerr<<"+++\n+++ Number of particles in event "<TreeH(); - MakeSeedsFromHits(tpc,th,sArray); - } else { - // Get TreeTR with track references - TTree *ttr = rl->TreeTR(); - MakeSeedsFromRefs(ttr,sArray); - } - - - nSeeds = sArray.GetEntries(); - cerr<<"\n\n+++\n+++ Number of seeds: "<GetMCApp()->Particle(l); + pdgCodes[l] = part->GetPdgCode(); + done[l] = kFALSE; + } - // Get track label - label = seed->GetLabel(); - - // 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); + printf("+++ Number of particles: %d\n",nParticles); + // 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); + } - 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 = covTreePi[bin]; - break; - case 321: - fCovTree = covTreeKa[bin]; - break; - case 2212: - fCovTree = covTreePr[bin]; - break; - case 11: - fCovTree = covTreeEl[bin]; - break; - case 13: - fCovTree = covTreeMu[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); + 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); - // 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(); - - arrentr = tArray.GetEntriesFast(); - for(Int_t l=0; lFill(); + // assign the track a dE/dx and make a rough PID + CookdEdx(sPt,sEta); } - - - // write the tree with tracks in the output file - out->cd(); - tracktree->Write(); - delete tracktree; - delete [] done; - delete [] pdgCodes; - delete [] ptkine; - delete [] pzkine; - - printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks); - //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<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<cd(); - } // loop on events - - if(fileDB) fileDB->Close(); - return 0; } //----------------------------------------------------------------------------- @@ -1022,15 +1044,15 @@ void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) { //----------------------------------------------------------------------------- Double_t xref = s->GetXL(); Double_t xx[5],cc[15]; - cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.; + cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=0.; 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 = s->GetPt()*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 sMom,sPos; @@ -1051,56 +1073,24 @@ void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) { // 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] = s->GetYL(); xx[1] = s->GetZL(); + xx[2] = ch/rho*(xref-x0); xx[3] = s->GetPz()/s->GetPt(); - xx[4] = -ch/rho; - xx[2] = xx[4]*x0; + xx[4] = ch/rho; // create the object AliTPCtrack - AliTPCtrack track(0,xx,cc,xref,s->GetAlpha()); + AliTPCtrack track(xref,s->GetAlpha(),xx,cc,0); new(&fTrack) AliTPCtrack(track); return; } //----------------------------------------------------------------------------- -Int_t AliTPCtrackerParam::CheckLabel(AliTPCseedGeant *s,Int_t nPart, - Double_t *ptkine,Double_t *pzkine) const { -//----------------------------------------------------------------------------- -// This function checks if the label of the seed has been correctly -// assigned (to do only for pp charm production with AliRoot v3-08-02) -//----------------------------------------------------------------------------- - - Int_t sLabel = s->GetLabel(); - Double_t sPt = s->GetPt(); - Double_t sPz = s->GetPz(); - - // check if the label is correct (comparing momentum) - if(sLabel=nPart) return 1; - - Double_t diff=0,mindiff=1000.; - Int_t bestLabel=0; - - for(Int_t i=sLabel-30; i=nPart) continue; - diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]); - if(diff0.001) return 1; - s->SetLabel(bestLabel); - - return 0; -} -//----------------------------------------------------------------------------- void AliTPCtrackerParam::CompareTPCtracks( + Int_t nEvents, const Char_t* galiceName, const Char_t* trkGeaName, const Char_t* trkKalName, @@ -1119,7 +1109,7 @@ void AliTPCtrackerParam::CompareTPCtracks( TFile *galiceFile = TFile::Open(galiceName); // get the AliRun object - AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice"); + AliRun *lAlice = (AliRun*)galiceFile->Get("gAlice"); // create the tree for comparison results @@ -1138,7 +1128,7 @@ void AliTPCtrackerParam::CompareTPCtracks( Double_t ptgener; Bool_t usethis; Int_t label; - Double_t cc[15],dAlpha; + 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]; @@ -1165,12 +1155,12 @@ void AliTPCtrackerParam::CompareTPCtracks( Char_t tname[100]; // loop on events in file - for(Int_t evt=0; evtGetEvent(evt); + const Int_t knparticles = lAlice->GetEvent(evt); Int_t *kalLab = new Int_t[knparticles]; for(Int_t i=0; iGetEvent(j); label = geatrack->GetLabel(); - part = (TParticle*)gAlice->GetMCApp()->Particle(label); + part = (TParticle*)lAlice->GetMCApp()->Particle(label); // use only injected tracks with fixed values of pT ptgener = part->Pt(); @@ -1265,7 +1255,7 @@ void AliTPCtrackerParam::CompareTPCtracks( cmptrk.eta = part->Eta(); cmptrk.r = TMath::Sqrt(part->Vx()*part->Vx()+part->Vy()*part->Vy()); - cmptrk.pt = 1/TMath::Abs(geatrack->Get1Pt()); + cmptrk.pt = geatrack->Pt(); cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl())); cmptrk.p = cmptrk.pt/cmptrk.cosl; @@ -1302,30 +1292,30 @@ void AliTPCtrackerParam::CompareTPCtracks( cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY(); cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ(); - cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta(); + 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(); + cmptrk.dpt = 1/kaltrack->GetSigned1Pt()-1/geatrack->GetSigned1Pt(); // get covariance matrix // beware: lines 3 and 4 in the matrix are inverted! - kaltrack->GetCovariance(cc); - - cmptrk.c00 = cc[0]; - cmptrk.c10 = cc[1]; - cmptrk.c11 = cc[2]; - cmptrk.c20 = cc[3]; - cmptrk.c21 = cc[4]; - cmptrk.c22 = cc[5]; - cmptrk.c30 = cc[10]; - cmptrk.c31 = cc[11]; - cmptrk.c32 = cc[12]; - cmptrk.c33 = cc[14]; - cmptrk.c40 = cc[6]; - cmptrk.c41 = cc[7]; - cmptrk.c42 = cc[8]; - cmptrk.c43 = cc[13]; - cmptrk.c44 = cc[9]; + //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(); @@ -1379,7 +1369,7 @@ void AliTPCtrackerParam::CompareTPCtracks( WriteEffs(tpceffrootName); // delete AliRun object - delete gAlice; gAlice=0; + delete lAlice; lAlice=0; // close all input files kalFile->Close(); @@ -1423,24 +1413,28 @@ void AliTPCtrackerParam::CookdEdx(Double_t pt,Double_t eta) { //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(0.13957); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; } if (dEdx < 39.+ 12./p/p) { - t.AssignMass(0.49368); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massKa); new(&fTrack) AliTPCtrack(t); return; } - t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return; } if (p<1.2) { if (dEdx < 39.+ 12./(p+0.25)/(p+0.25)) { - t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; } - t.AssignMass(0.93827); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massPr); new(&fTrack) AliTPCtrack(t); return; } - t.AssignMass(0.13957); new(&fTrack) AliTPCtrack(t); return; + t.AssignMass(massPi); new(&fTrack) AliTPCtrack(t); return; } //----------------------------------------------------------------------------- void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) { @@ -1462,7 +1456,7 @@ void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) { // get P and Cosl from track cosl = TMath::Cos(TMath::ATan(fTrack.GetTgl())); - p = 1./TMath::Abs(fTrack.Get1Pt())/cosl; + p = fTrack.Pt()/cosl; trkKine[0] = p; @@ -1490,10 +1484,10 @@ void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) { 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); - + cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar); + + + TMatrixD covMatSmear(5,5); covMatSmear = GetSmearingMatrix(cc,pt,eta); // get track original parameters @@ -1501,14 +1495,15 @@ void AliTPCtrackerParam::CookTrack(Double_t pt,Double_t eta) { alpha=fTrack.GetAlpha(); xx[0]=fTrack.GetY(); xx[1]=fTrack.GetZ(); - xx[2]=fTrack.GetX()*fTrack.GetC()-fTrack.GetSnp(); + 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(0,xxsm,cc,xref,alpha); + AliTPCtrack track(xref,alpha,xxsm,cc,0); new(&fTrack) AliTPCtrack(track); return; @@ -1889,20 +1884,20 @@ void AliTPCtrackerParam::MakeDataBase() { // create the trees for cov. matrices // trees for pions - TTree *covTreePi_ = NULL; - covTreePi_ = new TTree[knBinsPi]; + TTree *covTreePi1 = NULL; + covTreePi1 = new TTree[knBinsPi]; // trees for kaons - TTree *covTreeKa_ = NULL; - covTreeKa_ = new TTree[knBinsKa]; + TTree *covTreeKa1 = NULL; + covTreeKa1 = new TTree[knBinsKa]; // trees for protons - TTree *covTreePr_ = NULL; - covTreePr_ = new TTree[knBinsPr]; + TTree *covTreePr1 = NULL; + covTreePr1 = new TTree[knBinsPr]; // trees for electrons - TTree *covTreeEl_ = NULL; - covTreeEl_ = new TTree[knBinsEl]; + TTree *covTreeEl1 = NULL; + covTreeEl1 = new TTree[knBinsEl]; // trees for muons - TTree *covTreeMu_ = NULL; - covTreeMu_ = new TTree[knBinsMu]; + TTree *covTreeMu1 = NULL; + covTreeMu1 = new TTree[knBinsMu]; Char_t hname[100], htitle[100]; COVMATRIX covmat; @@ -1911,32 +1906,32 @@ void AliTPCtrackerParam::MakeDataBase() { for(Int_t i=0; icd("/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; @@ -2090,8 +2085,8 @@ void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th, Int_t label; Int_t nTracks=(Int_t)th->GetEntries(); - cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<SetAddress(&tkRefArray); Int_t nTkRef = (Int_t)b->GetEntries(); - cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<GetEvent(l)) continue; nnn = tkRefArray->GetEntriesFast(); @@ -2182,7 +2176,8 @@ void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr, AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label); // reject if not at the inner part of TPC - if(TMath::Abs(ioseed->GetXL()-82.9701) > 0.1) { + if(TMath::Abs(ioseed->GetXL()-83.8) > 0.2) { + //cerr<<"not at TPC inner part: XL = "<GetXL()<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; } @@ -2545,16 +2609,83 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) { 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) { @@ -2563,10 +2694,66 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) { 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) { @@ -2575,6 +2762,34 @@ Int_t AliTPCtrackerParam::ReadRegParams(const Char_t* inName,Int_t pdg) { 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(); @@ -3218,15 +3433,25 @@ void AliTPCtrackerParam::SmearTrack(Double_t *xx,Double_t *xxsm,TMatrixD cov) //----------------------------------------------------------------------------- // 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]); + } + delete corgen; corgen = 0; - for(Int_t l=0;l<5;l++) { - xxsm[l] = xx[l]+corr[l]; - } + for(Int_t l=0;l<5;l++) xxsm[l] = xx[l]+corr[l]; return; }