* Output file contains sorted tracks, ready for matching with ITS. *
* *
* For details: *
- * Alice Internal Note (submitted - get it from andrea.dainese@pd.infn.it)*
- * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm *
+ * Alice Internal Note 2003-011 *
* *
* Test macro is: AliBarrelRec_TPCparam.C *
* *
* - 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
+//
+//
+// *
//------- Root headers --------
+#include <Riostream.h>
+#include <TCanvas.h>
#include <TChain.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
+#include <TH2F.h>
#include <TLegend.h>
#include <TLine.h>
#include <TMatrixD.h>
-#include <TNtuple.h>
+#include <TParticle.h>
+#include <TStyle.h>
#include <TSystem.h>
-#include <TH2F.h>
+#include <TFile.h>
+#include <TRandom.h>
//------ 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) {
ClassImp(AliTPCtrackerParam)
//-----------------------------------------------------------------------------
-AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz,
- const Int_t n, const char* evfoldname):
- fEvFolderName(evfoldname) {
+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
//-----------------------------------------------------------------------------
- fNevents = n; // events to be processed
- fBz = Bz; // value of the z component of L3 field (Tesla)
- fColl = coll; // 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"<<endl;
+ if(fBz!=0.4 && fBz!=0.5) {
+ Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid field!\n Available: 0.4 or 0.5");
}
- if(fColl!=0 && fColl!=1) {
- cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n";
- cerr<<" Available: 0 -> PbPb6000"<<endl;
- cerr<<" 1 -> pp"<<endl;
+ if(fColl!=0 && fColl!=1) {
+ Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n Available: 0 -> PbPb6000\n 1 -> pp");
}
fDBfileName = gSystem->Getenv("ALICE_ROOT");
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::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; l<nBinsPi; l++) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- str+=l;
- covTreePi[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsKa = fDBgridKa.GetTotBins();
- for(Int_t l=0; l<nBinsKa; l++) {
- str = "/CovMatrices/Kaons/CovTreeKa_bin";
- str+=l;
- covTreeKa[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsPr = fDBgridPr.GetTotBins();
- for(Int_t l=0; l<nBinsPr; l++) {
- if(fColl==0) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- } else {
- str = "/CovMatrices/Protons/CovTreePr_bin";
- }
- str+=l;
- covTreePr[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsEl = fDBgridEl.GetTotBins();
- for(Int_t l=0; l<nBinsEl; l++) {
- str = "/CovMatrices/Electrons/CovTreeEl_bin";
- str+=l;
- covTreeEl[l] = (TTree*)fileDB->Get(str.Data());
- }
- Int_t nBinsMu = fDBgridMu.GetTotBins();
- for(Int_t l=0; l<nBinsMu; l++) {
- if(fColl==0) {
- str = "/CovMatrices/Pions/CovTreePi_bin";
- } else {
- str = "/CovMatrices/Muons/CovTreeMu_bin";
- }
- str+=l;
- covTreeMu[l] = (TTree*)fileDB->Get(str.Data());
- }
- } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
+ // 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) {
+ 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;
+ }
+
+ // Set the conversion constant between curvature and Pt
+ AliTracker::SetFieldMap(fiel,kTRUE);
+
+ 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();
+ if(!event) return -1;
- // Get gAlice object from file
+ 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!"<<endl;
- cerr<<"Field selected is: "<<fBz<<" T\n";
- cerr<<"Field found on file is: "<<fieval<<" T\n";
- return 0;
- }
-
// Get TPC detector
- AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
- Int_t ver = TPC->IsVersion();
- cerr<<"+++ TPC version "<<ver<<" has been found !\n";
+ AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
rl->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);
+ tpc->SetParam(digp);
- // Set the conversion constant between curvature and Pt
- AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
-
- TParticle *Part=0;
+ 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; evt<fNevents; evt++){
- cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
- tracks=0;
-
- // tree for TPC tracks
- sprintf(tname,"TreeT_TPC_%d",evt);
- TTree *tracktree = new TTree(tname,"Tree with TPC tracks");
- tracktree->Branch("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; l<nParticles; l++) {
- Part = (TParticle*)gAlice->GetMCApp()->Particle(l);
- pdgCodes[l] = Part->GetPdgCode();
- ptkine[l] = Part->Pt();
- pzkine[l] = Part->Pz();
- done[l] = kFALSE;
- }
- cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nParticles<<
- "\n+++\n";
-
- cerr<<"\n ********** MAKING SEEDS *****************"<<endl<<endl;
-
-
- // 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
- TTree *TTR = rl->TreeTR();
- MakeSeedsFromRefs(TTR,sArray);
- }
-
-
- nSeeds = sArray.GetEntries();
- cerr<<"\n\n+++\n+++ Number of seeds: "<<nSeeds<<"\n+++\n";
-
-
- cerr<<"\n ********** BUILDING TRACKS **************"<<endl<<endl;
-
- // loop over entries in sArray
- for(Int_t l=0; l<nSeeds; l++) {
- //if(l%1000==0) cerr<<" --- Processing seed "
- // <<l<<" of "<<nSeeds<<" ---\r";
-
- seed = (AliTPCseedGeant*)sArray.At(l);
-
- // this is TEMPORARY: only for reconstruction of pp production for charm
- if(fColl==1) cl = CheckLabel(seed,nParticles,ptkine,pzkine);
- if(cl) continue;
+ 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; l<nParticles; l++) {
+ part = (TParticle*)gAlice->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; l<nSeeds; l++) {
+ //if(l%1==0) printf(" --- Processing seed %d of %d ---\n",l,nSeeds);
+
+ seed = (AliTPCseedGeant*)sArray.At(l);
+
+ // 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);
+
+ 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; l<arrentr; l++) {
- tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
- tracktree->Fill();
+ // 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<<endl;
-
- sArray.Delete();
- tArray.Delete();
+ // 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; l<arrentr; l++) {
+ tpctrack=(AliTPCtrack*)tArray.UncheckedAt(l);
+ AliESDtrack ioESDtrack;
+ ioESDtrack.UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
+ // rough PID
+ wgts[0]=0.; wgts[1]=0.; wgts[2]=0.; wgts[3]=0.; wgts[4]=0.;
+ if(TMath::Abs(tpctrack->GetMass()-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<<endl;
+
+ sArray.Delete();
+ tArray.Delete();
- infile->cd();
- } // loop on events
-
- if(fileDB) fileDB->Close();
-
return 0;
}
//-----------------------------------------------------------------------------
gStyle->SetOptStat(0);
gStyle->SetOptFit(10001);
- Char_t *part="PIONS";
+ const char *part="PIONS";
Double_t ymax=500.;
/*
SetParticle(pdg);
- const Int_t nTotBins = fDBgrid->GetTotBins();
+ const Int_t knTotBins = fDBgrid->GetTotBins();
- cerr<<" Fit bins: "<<nTotBins<<endl;
+ cerr<<" Fit bins: "<<knTotBins<<endl;
Int_t bin=0;
- Int_t *n = new Int_t[nTotBins];
- Double_t *p = new Double_t[nTotBins];
- Double_t *ep = new Double_t[nTotBins];
- Double_t *mean = new Double_t[nTotBins];
- Double_t *sigma = new Double_t[nTotBins];
+ 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];
- for(Int_t l=0; l<nTotBins; l++) {
+ for(Int_t l=0; l<knTotBins; l++) {
n[l] = 1; // set to 1 to avoid divisions by 0
p[l] = mean[l] = sigma[l] = ep[l] = 0.;
}
n[bin]++;
} // loop on chain entries
- for(Int_t l=0; l<nTotBins; l++) {
+ for(Int_t l=0; l<knTotBins; l++) {
p[l] /= n[l];
mean[l] /= n[l];
n[l] = 1; // set to 1 to avoid divisions by 0
sigma[bin] += (cmptrk.dEdx-mean[bin])*(cmptrk.dEdx-mean[bin]);
} // loop on chain entries
- for(Int_t l=0; l<nTotBins; l++) {
+ for(Int_t l=0; l<knTotBins; l++) {
sigma[l] = TMath::Sqrt(sigma[l]/n[l]);
}
TCanvas *canv = new TCanvas("canv","dEdx",0,0,900,700);
// create the graph for dEdx vs p
- TGraphErrors *gr = new TGraphErrors(nTotBins,p,mean,ep,sigma);
+ TGraphErrors *gr = new TGraphErrors(knTotBins,p,mean,ep,sigma);
TString title(" : dE/dx vs momentum"); title.Prepend(part);
TH2F *frame = new TH2F("frame1",title.Data(),2,0.1,50,2,0,ymax);
frame->SetXTitle("p [GeV/c]");
switch(pdg) {
case 211:
- for(Int_t i=0; i<nTotBins; i++) {
+ for(Int_t i=0; i<knTotBins; i++) {
fdEdxMeanPi.SetParam(i,mean[i]);
fdEdxRMSPi.SetParam(i,sigma[i]);
}
break;
case 321:
- for(Int_t i=0; i<nTotBins; i++) {
+ for(Int_t i=0; i<knTotBins; i++) {
fdEdxMeanKa.SetParam(i,mean[i]);
fdEdxRMSKa.SetParam(i,sigma[i]);
}
break;
case 2212:
- for(Int_t i=0; i<nTotBins; i++) {
+ for(Int_t i=0; i<knTotBins; i++) {
fdEdxMeanPr.SetParam(i,mean[i]);
fdEdxRMSPr.SetParam(i,sigma[i]);
}
break;
case 11:
- for(Int_t i=0; i<nTotBins; i++) {
+ for(Int_t i=0; i<knTotBins; i++) {
fdEdxMeanEl.SetParam(i,mean[i]);
fdEdxRMSEl.SetParam(i,sigma[i]);
}
break;
case 13:
- for(Int_t i=0; i<nTotBins; i++) {
+ for(Int_t i=0; i<knTotBins; i++) {
fdEdxMeanMu.SetParam(i,mean[i]);
fdEdxRMSMu.SetParam(i,sigma[i]);
}
cerr<<"nTotBins = "<<nTotBins<<endl;
// create histograms for the all the bins
- TH1F *hPulls0_=NULL;
- TH1F *hPulls1_=NULL;
- TH1F *hPulls2_=NULL;
- TH1F *hPulls3_=NULL;
- TH1F *hPulls4_=NULL;
+ TH1F *hPulls0=NULL;
+ TH1F *hPulls1=NULL;
+ TH1F *hPulls2=NULL;
+ TH1F *hPulls3=NULL;
+ TH1F *hPulls4=NULL;
- hPulls0_ = new TH1F[nTotBins];
- hPulls1_ = new TH1F[nTotBins];
- hPulls2_ = new TH1F[nTotBins];
- hPulls3_ = new TH1F[nTotBins];
- hPulls4_ = new TH1F[nTotBins];
+ hPulls0 = new TH1F[nTotBins];
+ hPulls1 = new TH1F[nTotBins];
+ hPulls2 = new TH1F[nTotBins];
+ hPulls3 = new TH1F[nTotBins];
+ hPulls4 = new TH1F[nTotBins];
for(Int_t i=0; i<nTotBins; i++) {
- sprintf(hname,"hPulls0_%d",i);
+ sprintf(hname,"hPulls0%d",i);
sprintf(htitle,"P0 pulls for bin %d",i);
hDum->SetName(hname); hDum->SetTitle(htitle);
- hPulls0_[i] = *hDum;
- sprintf(hname,"hPulls1_%d",i);
+ 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);
+ 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);
+ 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);
+ 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;
+ hPulls4[i] = *hDum;
}
// loop on chain entries
// fill histograms with the pulls
bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
//cerr<<" pt "<<cmptrk.pt<<" eta "<<cmptrk.eta<<" bin "<<bin<<endl;
- hPulls0_[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
- hPulls1_[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
- hPulls2_[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
- hPulls3_[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
- hPulls4_[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
+ hPulls0[bin].Fill(cmptrk.dP0/TMath::Sqrt(cmptrk.c00));
+ hPulls1[bin].Fill(cmptrk.dP1/TMath::Sqrt(cmptrk.c11));
+ hPulls2[bin].Fill(cmptrk.dP2/TMath::Sqrt(cmptrk.c22));
+ hPulls3[bin].Fill(cmptrk.dP3/TMath::Sqrt(cmptrk.c33));
+ hPulls4[bin].Fill(cmptrk.dP4/TMath::Sqrt(cmptrk.c44));
} // loop on chain entries
// compute the sigma of the distributions
for(Int_t i=0; i<nTotBins; i++) {
- if(hPulls0_[i].GetEntries()>10) {
- g->SetRange(-3.*hPulls0_[i].GetRMS(),3.*hPulls0_[i].GetRMS());
- hPulls0_[i].Fit("g","R,Q,N");
+ if(hPulls0[i].GetEntries()>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");
+ 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");
+ 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");
+ 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");
+ 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
break;
}
- delete [] hPulls0_;
- delete [] hPulls1_;
- delete [] hPulls2_;
- delete [] hPulls3_;
- delete [] hPulls4_;
+ delete [] hPulls0;
+ delete [] hPulls1;
+ delete [] hPulls2;
+ delete [] hPulls3;
+ delete [] hPulls4;
} // loop on particle species
SetParticle(pdg);
- const Int_t nPtBins = fEff->GetPointsPt();
- cerr<<"nPtBins = "<<nPtBins<<endl;
- Double_t *dP0 = new Double_t[nPtBins];
- Double_t *dP4 = new Double_t[nPtBins];
- Double_t *dPtToPt = new Double_t[nPtBins];
- Double_t *pt = new Double_t[nPtBins];
+ const Int_t knPtBins = fEff->GetPointsPt();
+ cerr<<"knPtBins = "<<knPtBins<<endl;
+ Double_t *dP0 = new Double_t[knPtBins];
+ Double_t *dP4 = new Double_t[knPtBins];
+ Double_t *dPtToPt = new Double_t[knPtBins];
+ Double_t *pt = new Double_t[knPtBins];
fEff->GetArrayPt(pt);
TF1 *g = new TF1("g","gaus");
// create histograms for the all the bins
- TH1F *hP0_=NULL;
- TH1F *hP4_=NULL;
- TH1F *hPt_=NULL;
+ TH1F *hP0=NULL;
+ TH1F *hP4=NULL;
+ TH1F *hPt=NULL;
- hP0_ = new TH1F[nPtBins];
- hP4_ = new TH1F[nPtBins];
- hPt_ = new TH1F[nPtBins];
+ hP0 = new TH1F[knPtBins];
+ hP4 = new TH1F[knPtBins];
+ hPt = new TH1F[knPtBins];
- for(Int_t i=0; i<nPtBins; i++) {
- hP0_[i] = *hDumP0;
- hP4_[i] = *hDumP4;
- hPt_[i] = *hDumPt;
+ for(Int_t i=0; i<knPtBins; i++) {
+ hP0[i] = *hDumP0;
+ hP4[i] = *hDumP4;
+ hPt[i] = *hDumPt;
}
// loop on chain entries
// fill histograms with the residuals
bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
//cerr<<" pt "<<cmptrk.pt<<" eta "<<cmptrk.eta<<" bin "<<bin<<endl;
- hP0_[bin].Fill(cmptrk.dP0);
- hP4_[bin].Fill(cmptrk.dP4);
- hPt_[bin].Fill(cmptrk.dpt/cmptrk.pt);
+ hP0[bin].Fill(cmptrk.dP0);
+ hP4[bin].Fill(cmptrk.dP4);
+ hPt[bin].Fill(cmptrk.dpt/cmptrk.pt);
} // loop on chain entries
cPtres->Divide(5,2);
// Draw histograms
- for(Int_t i=0; i<nPtBins; i++) {
- cP0res->cd(i+1); hP0_[i].Draw();
- cP4res->cd(i+1); hP4_[i].Draw();
- cPtres->cd(i+1); hPt_[i].Draw();
+ for(Int_t i=0; i<knPtBins; i++) {
+ cP0res->cd(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; i<nPtBins; i++) {
- if(hP0_[i].GetEntries()>10) {
- g->SetRange(-3.*hP0_[i].GetRMS(),3.*hP0_[i].GetRMS());
- hP0_[i].Fit("g","R,Q,N");
+ for(Int_t i=0; i<knPtBins; i++) {
+ if(hP0[i].GetEntries()>10) {
+ 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");
+ 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");
+ 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(nPtBins,pt,dP0);
- TGraph *grdP4 = new TGraph(nPtBins,pt,dP4);
- TGraph *grdPtToPt = new TGraph(nPtBins,pt,dPtToPt);
+ 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);
delete [] pt;
- delete [] hP0_;
- delete [] hP4_;
- delete [] hPt_;
+ delete [] hP0;
+ delete [] hP4;
+ delete [] hPt;
return;
}
//-----------------------------------------------------------------------------
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;
// 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 &&
- TMath::Abs(sPt-ptkine[sLabel])*
- TMath::Abs(sPz-pzkine[sLabel])<0.001) return 0;
-
- if((sLabel-30)>=nPart) return 1;
-
- Double_t diff=0,mindiff=1000.;
- Int_t bestLabel=0;
-
- for(Int_t i=sLabel-30; i<sLabel; i++) {
- if(i<0 || i>=nPart) continue;
- diff = TMath::Abs(sPt-ptkine[i])*TMath::Abs(sPz-pzkine[i]);
- if(diff<mindiff) { mindiff = diff; bestLabel = i; }
- }
-
- if(mindiff>0.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,
Double_t *pt = new Double_t[effBinsPt];
fEffPi.GetArrayPt(pt);
- TParticle *Part;
+ TParticle *part;
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];
Char_t tname[100];
// loop on events in file
- for(Int_t evt=0; evt<fNevents; evt++) {
+ for(Int_t evt=0; evt<nEvents; evt++) {
cerr<<"\n --- Reading tracks for event "<<evt<<" ---\n\n";
sprintf(tname,"TreeT_TPC_%d",evt);
// particles from TreeK
- const Int_t nparticles = gAlice->GetEvent(evt);
+ const Int_t knparticles = gAlice->GetEvent(evt);
- Int_t *kalLab = new Int_t[nparticles];
- for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
+ Int_t *kalLab = new Int_t[knparticles];
+ for(Int_t i=0; i<knparticles; i++) kalLab[i] = -1;
// tracks from Kalman
// Read the labels of the seeds
char sname[100];
Int_t sLabel,ncol;
- Bool_t *hasSeed = new Bool_t[nparticles];
- for(Int_t i=0; i<nparticles; i++) hasSeed[i] = kFALSE;
+ Bool_t *hasSeed = new Bool_t[knparticles];
+ for(Int_t i=0; i<knparticles; i++) hasSeed[i] = kFALSE;
sprintf(sname,"seedLabels.%d.dat",evt);
FILE *seedFile = fopen(sname,"r");
while(1) {
geatree->GetEvent(j);
label = geatrack->GetLabel();
- Part = (TParticle*)gAlice->GetMCApp()->Particle(label);
+ part = (TParticle*)gAlice->GetMCApp()->Particle(label);
// use only injected tracks with fixed values of pT
- ptgener = Part->Pt();
+ ptgener = part->Pt();
usethis = kFALSE;
for(Int_t l=0; l<fEffPi.GetPointsPt(); l++) {
if(TMath::Abs(ptgener-pt[l])<0.01) usethis = kTRUE;
if(out) continue;
*/
- cmptrk.pdg = Part->GetPdgCode();
- cmptrk.eta = Part->Eta();
- cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
+ 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.pt = geatrack->Pt();
cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
cmptrk.p = cmptrk.pt/cmptrk.cosl;
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();
//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) {
// 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;
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
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;
ReadEffs(inName);
SetParticle(pdg);
- const Int_t n = fEff->GetPointsPt();
- Double_t *effsA = new Double_t[n];
- Double_t *effsB = new Double_t[n];
- Double_t *effsC = new Double_t[n];
- Double_t *pt = new Double_t[n];
+ 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;i<n;i++) {
+ for(Int_t i=0;i<kn;i++) {
effsA[i] = fEff->GetParam(i,0);
effsB[i] = fEff->GetParam(i,1);
effsC[i] = fEff->GetParam(i,2);
}
- TGraph *grA = new TGraph(n,pt,effsA);
- TGraph *grB = new TGraph(n,pt,effsB);
- TGraph *grC = new TGraph(n,pt,effsC);
+ 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);
ReadPulls(inName);
SetParticle(pdg);
- const Int_t n = (fPulls+par)->GetPointsPt();
- Double_t *pullsA = new Double_t[n];
- Double_t *pullsB = new Double_t[n];
- Double_t *pullsC = new Double_t[n];
- Double_t *pt = new Double_t[n];
+ 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;i<n;i++) {
+ for(Int_t i=0;i<kn;i++) {
pullsA[i] = (fPulls+par)->GetParam(i,0);
pullsB[i] = (fPulls+par)->GetParam(i,1);
pullsC[i] = (fPulls+par)->GetParam(i,2);
}
- TGraph *grA = new TGraph(n,pt,pullsA);
- TGraph *grB = new TGraph(n,pt,pullsB);
- TGraph *grC = new TGraph(n,pt,pullsC);
+ 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);
// "dEdx" -> initialize fdEdx... members
//-----------------------------------------------------------------------------
- const char *DB = strstr(which,"DB");
+ const char *db = strstr(which,"DB");
const char *eff = strstr(which,"eff");
const char *pulls = strstr(which,"pulls");
const char *dEdx = strstr(which,"dEdx");
Double_t *eta=0,*pt=0;
- if(DB) {
+ if(db) {
nEta = 2;
nPt = 9;
eta = etaPoints;
AliTPCkineGrid *dummy=0;
- if(DB) {
+ if(db) {
dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
new(&fDBgridPi) AliTPCkineGrid(*dummy);
new(&fDBgridKa) AliTPCkineGrid(*dummy);
//-----------------------------------------------------------------------------
// define some file names
- Char_t *effFile ="TPCeff.root";
- Char_t *pullsFile ="pulls.root";
- Char_t *regPiFile ="regPi.root";
- Char_t *regKaFile ="regKa.root";
- Char_t *regPrFile ="regPr.root";
- Char_t *regElFile ="regEl.root";
- Char_t *regMuFile ="regMu.root";
- Char_t *dEdxPiFile="dEdxPi.root";
- Char_t *dEdxKaFile="dEdxKa.root";
- Char_t *dEdxPrFile="dEdxPr.root";
- Char_t *dEdxElFile="dEdxEl.root";
- Char_t *dEdxMuFile="dEdxMu.root";
- Char_t *cmFile ="CovMatrix_AllEvts.root";
+ 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";
//
InitializeKineGrid("DB");
- const Int_t nBinsPi = fDBgridPi.GetTotBins();
- const Int_t nBinsKa = fDBgridKa.GetTotBins();
- const Int_t nBinsPr = fDBgridPr.GetTotBins();
- const Int_t nBinsEl = fDBgridEl.GetTotBins();
- const Int_t nBinsMu = fDBgridMu.GetTotBins();
+ 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 *CovTreePi_ = NULL;
- CovTreePi_ = new TTree[nBinsPi];
+ TTree *covTreePi1 = NULL;
+ covTreePi1 = new TTree[knBinsPi];
// trees for kaons
- TTree *CovTreeKa_ = NULL;
- CovTreeKa_ = new TTree[nBinsKa];
+ TTree *covTreeKa1 = NULL;
+ covTreeKa1 = new TTree[knBinsKa];
// trees for protons
- TTree *CovTreePr_ = NULL;
- CovTreePr_ = new TTree[nBinsPr];
+ TTree *covTreePr1 = NULL;
+ covTreePr1 = new TTree[knBinsPr];
// trees for electrons
- TTree *CovTreeEl_ = NULL;
- CovTreeEl_ = new TTree[nBinsEl];
+ TTree *covTreeEl1 = NULL;
+ covTreeEl1 = new TTree[knBinsEl];
// trees for muons
- TTree *CovTreeMu_ = NULL;
- CovTreeMu_ = new TTree[nBinsMu];
+ TTree *covTreeMu1 = NULL;
+ covTreeMu1 = new TTree[knBinsMu];
Char_t hname[100], htitle[100];
COVMATRIX covmat;
- for(Int_t i=0; i<nBinsPi; i++) {
+ for(Int_t i=0; i<knBinsPi; i++) {
sprintf(hname,"CovTreePi_bin%d",i);
sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
- CovTreePi_[i].SetName(hname); CovTreePi_[i].SetTitle(htitle);
- CovTreePi_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
+ covTreePi1[i].SetName(hname); covTreePi1[i].SetTitle(htitle);
+ covTreePi1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",5000000);
}
- for(Int_t i=0; i<nBinsKa; i++) {
+ for(Int_t i=0; i<knBinsKa; i++) {
sprintf(hname,"CovTreeKa_bin%d",i);
sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
- CovTreeKa_[i].SetName(hname); CovTreeKa_[i].SetTitle(htitle);
- CovTreeKa_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
+ covTreeKa1[i].SetName(hname); covTreeKa1[i].SetTitle(htitle);
+ covTreeKa1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
}
- for(Int_t i=0; i<nBinsPr; i++) {
+ for(Int_t i=0; i<knBinsPr; i++) {
sprintf(hname,"CovTreePr_bin%d",i);
sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
- CovTreePr_[i].SetName(hname); CovTreePr_[i].SetTitle(htitle);
- CovTreePr_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
+ covTreePr1[i].SetName(hname); covTreePr1[i].SetTitle(htitle);
+ covTreePr1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
}
- for(Int_t i=0; i<nBinsEl; i++) {
+ for(Int_t i=0; i<knBinsEl; i++) {
sprintf(hname,"CovTreeEl_bin%d",i);
sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
- CovTreeEl_[i].SetName(hname); CovTreeEl_[i].SetTitle(htitle);
- CovTreeEl_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
+ covTreeEl1[i].SetName(hname); covTreeEl1[i].SetTitle(htitle);
+ covTreeEl1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
}
- for(Int_t i=0; i<nBinsMu; i++) {
+ for(Int_t i=0; i<knBinsMu; i++) {
sprintf(hname,"CovTreeMu_bin%d",i);
sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
- CovTreeMu_[i].SetName(hname); CovTreeMu_[i].SetTitle(htitle);
- CovTreeMu_[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
+ covTreeMu1[i].SetName(hname); covTreeMu1[i].SetTitle(htitle);
+ covTreeMu1[i].Branch("matrix",&covmat,"c00/D:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44",1000000);
}
/*
Int_t trkPdg,trkBin;
Double_t trkKine[1],trkRegPar[3];
- Int_t *nPerBinPi = new Int_t[nBinsPi];
- for(Int_t k=0;k<nBinsPi;k++) nPerBinPi[k]=0;
- Int_t *nPerBinKa = new Int_t[nBinsKa];
- for(Int_t k=0;k<nBinsKa;k++) nPerBinKa[k]=0;
- Int_t *nPerBinMu = new Int_t[nBinsMu];
- for(Int_t k=0;k<nBinsMu;k++) nPerBinMu[k]=0;
- Int_t *nPerBinEl = new Int_t[nBinsEl];
- for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
- Int_t *nPerBinPr = new Int_t[nBinsPr];
- for(Int_t k=0;k<nBinsPr;k++) nPerBinPr[k]=0;
+ Int_t *nPerBinPi = new Int_t[knBinsPi];
+ for(Int_t k=0;k<knBinsPi;k++) nPerBinPi[k]=0;
+ Int_t *nPerBinKa = new Int_t[knBinsKa];
+ for(Int_t k=0;k<knBinsKa;k++) nPerBinKa[k]=0;
+ Int_t *nPerBinMu = new Int_t[knBinsMu];
+ for(Int_t k=0;k<knBinsMu;k++) nPerBinMu[k]=0;
+ Int_t *nPerBinEl = new Int_t[knBinsEl];
+ for(Int_t k=0;k<knBinsEl;k++) nPerBinEl[k]=0;
+ Int_t *nPerBinPr = new Int_t[knBinsPr];
+ for(Int_t k=0;k<knBinsPr;k++) nPerBinPr[k]=0;
// loop on chain entries
for(Int_t l=0; l<entries; l++) {
// fill the tree
switch (trkPdg) {
case 211: // pions
- CovTreePi_[trkBin].Fill();
+ covTreePi1[trkBin].Fill();
nPerBinPi[trkBin]++;
break;
case 321: // kaons
- CovTreeKa_[trkBin].Fill();
+ covTreeKa1[trkBin].Fill();
nPerBinKa[trkBin]++;
break;
case 2212: // protons
- CovTreePr_[trkBin].Fill();
+ covTreePr1[trkBin].Fill();
nPerBinPr[trkBin]++;
break;
case 11: // electrons
- CovTreeEl_[trkBin].Fill();
+ covTreeEl1[trkBin].Fill();
nPerBinEl[trkBin]++;
break;
case 13: // muons
- CovTreeMu_[trkBin].Fill();
+ 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");
+ TFile *dbfile = new TFile(fDBfileName.Data(),"update");
+ dbfile->mkdir("CovMatrices");
gDirectory->cd("/CovMatrices");
gDirectory->mkdir("Pions");
gDirectory->mkdir("Kaons");
// store pions
gDirectory->cd("/CovMatrices/Pions");
fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
- for(Int_t i=0;i<nBinsPi;i++) CovTreePi_[i].Write();
+ for(Int_t i=0;i<knBinsPi;i++) covTreePi1[i].Write();
// store kaons
gDirectory->cd("/CovMatrices/Kaons");
fDBgridKa.SetName("DBgridKa"); fDBgridKa.Write();
- for(Int_t i=0;i<nBinsKa;i++) CovTreeKa_[i].Write();
+ for(Int_t i=0;i<knBinsKa;i++) covTreeKa1[i].Write();
// store kaons
gDirectory->cd("/CovMatrices/Protons");
fDBgridPr.SetName("DBgridPr"); fDBgridPr.Write();
- for(Int_t i=0;i<nBinsPr;i++) CovTreePr_[i].Write();
+ for(Int_t i=0;i<knBinsPr;i++) covTreePr1[i].Write();
// store electrons
gDirectory->cd("/CovMatrices/Electrons");
fDBgridEl.SetName("DBgridEl"); fDBgridEl.Write();
- for(Int_t i=0;i<nBinsEl;i++) CovTreeEl_[i].Write();
+ for(Int_t i=0;i<knBinsEl;i++) covTreeEl1[i].Write();
// store kaons
gDirectory->cd("/CovMatrices/Muons");
fDBgridMu.SetName("DBgridMu"); fDBgridMu.Write();
- for(Int_t i=0;i<nBinsMu;i++) CovTreeMu_[i].Write();
+ for(Int_t i=0;i<knBinsMu;i++) covTreeMu1[i].Write();
- DBfile->Close();
+ dbfile->Close();
delete [] nPerBinPi;
delete [] nPerBinKa;
delete [] nPerBinPr;
return;
}
//-----------------------------------------------------------------------------
-void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *TPC,TTree *TH,
+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();
+ Int_t nTracks=(Int_t)th->GetEntries();
- cerr<<"+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
- "\n+++\n\n";
+ cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
+ "\n";
AliTPChit *tpcHit=0;
for(Int_t l=0; l<nTracks; l++) {
if(l%1000==0) cerr<<" --- Processing primary track "
<<l<<" of "<<nTracks<<" ---\r";
- TPC->ResetHits();
- TH->GetEvent(l);
+ tpc->ResetHits();
+ th->GetEvent(l);
// Get FirstHit
- tpcHit=(AliTPChit*)TPC->FirstHit(-1);
- for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
+ 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();
// Get track label
label=tpcHit->Track();
- if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
+ 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();
return;
}
//-----------------------------------------------------------------------------
-void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *TTR,
+void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
TObjArray &seedArray) const {
//-----------------------------------------------------------------------------
// This function makes the seeds for tracks from the track references
TClonesArray *tkRefArray = new TClonesArray("AliTrackReference");
- TBranch *b =(TBranch*)TTR->GetBranch("TPC");
+ TBranch *b =(TBranch*)ttr->GetBranch("TPC");
if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
b->SetAddress(&tkRefArray);
Int_t nTkRef = (Int_t)b->GetEntries();
- cerr<<"+++\n+++ Number of entries in TreeTR(TPC): "<<nTkRef<<
- "\n+++\n\n";
+ cerr<<"+++ Number of entries in TreeTR(TPC): "<<nTkRef<<"\n";
// loop on track references
for(Int_t l=0; l<nTkRef; l++){
- if(l%1000==0) cerr<<" --- Processing primary track "
- <<l<<" of "<<nTkRef<<" ---\r";
+ //if(l%1000==0) cerr<<" --- Processing primary track "
+ // <<l<<" of "<<nTkRef<<" ---\r";
if(!b->GetEvent(l)) continue;
nnn = tkRefArray->GetEntriesFast();
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 = "<<ioseed->GetXL()<<endl;
delete ioseed; continue;
}
// 2) computes the average TPC efficiencies
//-----------------------------------------------------------------------------
- Char_t *outName="TPCeff.root";
+ const char *outName="TPCeff.root";
// merge files with tracks
cerr<<" ******** MERGING FILES **********\n\n";
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: "<<inName<<" not found\n";
+ return 0; }
+
+ TString str;
+
+ TFile *inFile = TFile::Open(inName);
+
+
+ Int_t nBinsPi = fDBgridPi.GetTotBins();
+ for(Int_t l=0; l<nBinsPi; l++) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ str+=l;
+ fCovTreePi[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreePi[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsKa = fDBgridKa.GetTotBins();
+ for(Int_t l=0; l<nBinsKa; l++) {
+ str = "/CovMatrices/Kaons/CovTreeKa_bin";
+ str+=l;
+ fCovTreeKa[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeKa[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsPr = fDBgridPr.GetTotBins();
+ for(Int_t l=0; l<nBinsPr; l++) {
+ if(fColl==0) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ } else {
+ str = "/CovMatrices/Protons/CovTreePr_bin";
+ }
+ str+=l;
+ fCovTreePr[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreePr[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsEl = fDBgridEl.GetTotBins();
+ for(Int_t l=0; l<nBinsEl; l++) {
+ str = "/CovMatrices/Electrons/CovTreeEl_bin";
+ str+=l;
+ fCovTreeEl[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeEl[l] = new TTree(*fCovTree);
+ }
+ Int_t nBinsMu = fDBgridMu.GetTotBins();
+ for(Int_t l=0; l<nBinsMu; l++) {
+ if(fColl==0) {
+ str = "/CovMatrices/Pions/CovTreePi_bin";
+ } else {
+ str = "/CovMatrices/Muons/CovTreeMu_bin";
+ }
+ str+=l;
+ fCovTreeMu[l] = (TTree*)inFile->Get(str.Data());
+ //fCovTree = (TTree*)inFile->Get(str.Data());
+ //fCovTreeMu[l] = new TTree(*fCovTree);
+ }
+
+ //inFile->Close();
return 1;
}
if(gSystem->AccessPathName(inName,kFileExists)) {
cerr<<"AliTPCtrackerParam::ReadRegParams: "<<inName<<" not found\n";
return 0; }
+ // Introduced change to "reverse" the TMatrixD read from file.
+ // Needed because storage mode of TMatrixD changed from column-wise
+ // to rwo-wise in ROOT.
+ //
+ // A.Dainese 03/06/05
+
+ TMatrixD dummyMatrix(9,3);
TFile *inFile = TFile::Open(inName);
switch (pdg) {
case 211:
fRegPar = (TMatrixD*)inFile->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/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/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();
gStyle->SetOptFit(10001);
Int_t thisPdg=211;
- Char_t *part="Pions - ";
+ const char *part="Pions - ";
InitializeKineGrid("DB");
SetParticle(pdg);
- const Int_t fitbins = fDBgrid->GetBinsPt();
- cerr<<" Fit bins: "<<fitbins<<endl;
+ const Int_t kfitbins = fDBgrid->GetBinsPt();
+ cerr<<" Fit bins: "<<kfitbins<<endl;
switch (pdg) {
case 211: // pions
Int_t pbin;
- Int_t *n = new Int_t[fitbins];
- Int_t *n00 = new Int_t[fitbins];
- Int_t *n11 = new Int_t[fitbins];
- Int_t *n20 = new Int_t[fitbins];
- Int_t *n22 = new Int_t[fitbins];
- Int_t *n31 = new Int_t[fitbins];
- Int_t *n33 = new Int_t[fitbins];
- Int_t *n40 = new Int_t[fitbins];
- Int_t *n42 = new Int_t[fitbins];
- Int_t *n44 = new Int_t[fitbins];
- Double_t *p = new Double_t[fitbins];
- Double_t *ep = new Double_t[fitbins];
- Double_t *mean00 = new Double_t[fitbins];
- Double_t *mean11 = new Double_t[fitbins];
- Double_t *mean20 = new Double_t[fitbins];
- Double_t *mean22 = new Double_t[fitbins];
- Double_t *mean31 = new Double_t[fitbins];
- Double_t *mean33 = new Double_t[fitbins];
- Double_t *mean40 = new Double_t[fitbins];
- Double_t *mean42 = new Double_t[fitbins];
- Double_t *mean44 = new Double_t[fitbins];
- Double_t *sigma00 = new Double_t[fitbins];
- Double_t *sigma11 = new Double_t[fitbins];
- Double_t *sigma20 = new Double_t[fitbins];
- Double_t *sigma22 = new Double_t[fitbins];
- Double_t *sigma31 = new Double_t[fitbins];
- Double_t *sigma33 = new Double_t[fitbins];
- Double_t *sigma40 = new Double_t[fitbins];
- Double_t *sigma42 = new Double_t[fitbins];
- Double_t *sigma44 = new Double_t[fitbins];
- Double_t *rmean = new Double_t[fitbins];
- Double_t *rsigma = new Double_t[fitbins];
+ 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; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
n[l]=1;
n00[l]=n11[l]=n20[l]=n22[l]=n31[l]=n33[l]=n40[l]=n42[l]=n44[l]=1;
p[l ]=ep[l]=0.;
mean44[pbin]+=cmptrk.c44;
} // loop on chain entries
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
p[l]/=n[l];
mean00[l]/=n[l];
mean11[l]/=n[l];
sigma44[pbin]+=(cmptrk.c44-mean44[pbin])*(cmptrk.c44-mean44[pbin]); }
} // loop on chain entries
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
sigma00[l] = TMath::Sqrt(sigma00[l]/n00[l]);
sigma11[l] = TMath::Sqrt(sigma11[l]/n11[l]);
sigma20[l] = TMath::Sqrt(sigma20[l]/n20[l]);
TCanvas *canv00 = new TCanvas("canv00","c00",0,0,700,900);
canv00->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr00 = new TGraphErrors(fitbins,p,mean00,ep,sigma00);
+ 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]");
gr00->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean00[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma00[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr00reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr00reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B})");
regtitle00.Prepend(part);
TH2F *frame00reg = new TH2F("frame00reg",regtitle00.Data(),2,0.1,50,2,0,2);
TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
canv11->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
+ 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]");
gr11->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean11[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma11[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr11reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr11reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B})");
regtitle11.Prepend(part);
TH2F *frame11reg = new TH2F("frame11reg",regtitle11.Data(),2,0.1,50,2,0,2);
TCanvas *canv20 = new TCanvas("canv20","c20",0,0,700,900);
canv20->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr20 = new TGraphErrors(fitbins,p,mean20,ep,sigma20);
+ 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]");
gr20->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean20[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma20[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr20reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr20reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B})");
regtitle20.Prepend(part);
TH2F *frame20reg = new TH2F("frame20reg",regtitle20.Data(),2,0.1,50,2,0,2);
TCanvas *canv22 = new TCanvas("canv22","c22",0,0,700,900);
canv22->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr22 = new TGraphErrors(fitbins,p,mean22,ep,sigma22);
+ 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]");
gr22->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean22[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma22[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr22reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr22reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B})");
regtitle22.Prepend(part);
TH2F *frame22reg = new TH2F("frame22reg",regtitle22.Data(),2,0.1,50,2,0,2);
TCanvas *canv31 = new TCanvas("canv31","c31",0,0,700,900);
canv31->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr31 = new TGraphErrors(fitbins,p,mean31,ep,sigma31);
+ 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]");
gr31->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean31[l]/RegFunc(&p[l],fitpar);
rsigma[l] = -sigma31[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr31reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr31reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B})");
regtitle31.Prepend(part);
TH2F *frame31reg = new TH2F("frame31reg",regtitle31.Data(),2,0.1,50,2,0,2);
TCanvas *canv33 = new TCanvas("canv33","c33",0,0,700,900);
canv33->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr33 = new TGraphErrors(fitbins,p,mean33,ep,sigma33);
+ 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]");
gr33->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean33[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma33[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr33reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr33reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B})");
regtitle33.Prepend(part);
TH2F *frame33reg = new TH2F("frame33reg",regtitle33.Data(),2,0.1,50,2,0,2);
TCanvas *canv40 = new TCanvas("canv40","c40",0,0,700,900);
canv40->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr40 = new TGraphErrors(fitbins,p,mean40,ep,sigma40);
+ 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]");
gr40->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean40[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma40[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr40reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr40reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B})");
regtitle40.Prepend(part);
TH2F *frame40reg = new TH2F("frame40reg",regtitle40.Data(),2,0.1,50,2,0,2);
TCanvas *canv42 = new TCanvas("canv42","c42",0,0,700,900);
canv42->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr42 = new TGraphErrors(fitbins,p,mean42,ep,sigma42);
+ 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]");
gr42->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean42[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma42[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr42reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr42reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B})");
regtitle42.Prepend(part);
TH2F *frame42reg = new TH2F("frame42reg",regtitle42.Data(),2,0.1,50,2,0,2);
TCanvas *canv44 = new TCanvas("canv44","c44",0,0,700,900);
canv44->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr44 = new TGraphErrors(fitbins,p,mean44,ep,sigma44);
+ 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]");
gr44->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
+ for(Int_t l=0; l<kfitbins; l++) {
rmean[l] = mean44[l]/RegFunc(&p[l],fitpar);
rsigma[l] = sigma44[l]/RegFunc(&p[l],fitpar);
}
// create the graph the regularized cov. matrix
- TGraphErrors *gr44reg = new TGraphErrors(fitbins,p,rmean,ep,rsigma);
+ TGraphErrors *gr44reg = new TGraphErrors(kfitbins,p,rmean,ep,rsigma);
TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B})");
regtitle44.Prepend(part);
TH2F *frame44reg = new TH2F("frame44reg",regtitle44.Data(),2,0.1,50,2,0,2);
//-----------------------------------------------------------------------------
// 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;
}
//-----------------------------------------------------------------------------
Option_t *opt;
- Char_t *dirName="Pions";
- Char_t *meanName="dEdxMeanPi";
- Char_t *RMSName="dEdxRMSPi";
+ const char *dirName="Pions";
+ const char *meanName="dEdxMeanPi";
+ const char *rmsName="dEdxRMSPi";
SetParticle(pdg);
case 211:
dirName="Pions";
meanName="dEdxMeanPi";
- RMSName="dEdxRMSPi";
+ rmsName="dEdxRMSPi";
break;
case 321:
dirName="Kaons";
meanName="dEdxMeanKa";
- RMSName="dEdxRMSKa";
+ rmsName="dEdxRMSKa";
break;
case 2212:
dirName="Protons";
meanName="dEdxMeanPr";
- RMSName="dEdxRMSPr";
+ rmsName="dEdxRMSPr";
break;
case 11:
dirName="Electrons";
meanName="dEdxMeanEl";
- RMSName="dEdxRMSEl";
+ rmsName="dEdxRMSEl";
break;
case 13:
dirName="Muons";
meanName="dEdxMeanMu";
- RMSName="dEdxRMSMu";
+ rmsName="dEdxRMSMu";
break;
}
TDirectory *dir2 = gDirectory->mkdir(dirName);
dir2->cd();
fdEdxMean->SetName(meanName); fdEdxMean->Write();
- fdEdxRMS->SetName(RMSName); fdEdxRMS->Write();
+ fdEdxRMS->SetName(rmsName); fdEdxRMS->Write();
outFile->Close();
delete outFile;
//-----------------------------------------------------------------------------
Option_t *opt;
- Char_t *dirName="Pions";
- Char_t *keyName="RegPions";
+ const char *dirName="Pions";
+ const char *keyName="RegPions";
SetParticle(pdg);
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff