* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.8 2002/05/08 18:21:40 kowal2
-Now the code is blind to the binning used for pulls, efficiencies etc.
-
-Revision 1.7 2002/04/10 16:30:07 kowal2
-logs added
-
-*/
-
+/* $Id$ */
/**************************************************************************
* *
* Output file contains sorted tracks, ready for matching with ITS. *
* *
* For details: *
- * http://www.pd.infn.it/alipd/talks/soft/adIII02/TPCtrackingParam.htm *
+ * Alice Internal Note 2003-011 *
* *
* Test macro is: AliBarrelRec_TPCparam.C *
* *
+ * 2002/10/01: Introduction of the database for pp collisions (B=0.4 T) *
+ * - Everything done separately for pions, kaons, protons, electrons and *
+ * muons. *
+ * - Now (only for pp) the tracks are built from the AliTrackReferences *
+ * which contain the position and momentum of all tracks at R = 83 cm; *
+ * This eliminates the loss of tracks due the dead zone of the TPC *
+ * where the 1st hit is not created. *
+ * - In AliBarrelRec_TPCparam.C there many possible ways of determining *
+ * the z coordinate of the primary vertex in pp events (from pixels, *
+ * from ITS tracks, smearing according to resolution given by tracks. *
+ * *
* 2002/04/28: Major upgrade of the class *
* - Covariance matrices and pulls are now separeted for pions, kaons *
* and electrons. *
* - 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 <TDatabasePDG.h>
+#include <TStyle.h>
#include <TSystem.h>
-#include "alles.h"
+#include <TFile.h>
+#include <TRandom.h>
+//------ AliRoot headers ------
#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 "AliTPC.h"
+#include "AliTPCParamSR.h"
#include "AliTPCkineGrid.h"
#include "AliTPCtrack.h"
#include "AliTPCtrackerParam.h"
+#include "AliTrackReference.h"
+#include "AliESDtrack.h"
+//-----------------------------
Double_t RegFunc(Double_t *x,Double_t *par) {
// This is the function used to regularize the covariance matrix
- Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
- TMath::Power(x[1],par[3]);
- return value;
-}
-Double_t FitRegFunc(Double_t *x,Double_t *par) {
-// This is the function used for the fit of the covariance
-// matrix as a function of the momentum.
-// For cosl the average value 0.87 is used.
- Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2])/
- TMath::Power(0.87,par[3]);
+ Double_t value = par[0]+par[1]/TMath::Power(x[0],par[2]);
return value;
}
ClassImp(AliTPCtrackerParam)
//-----------------------------------------------------------------------------
-AliTPCtrackerParam::AliTPCtrackerParam(const Int_t coll,const Double_t Bz)
+AliTPCtrackerParam::AliTPCtrackerParam(Int_t kcoll, Double_t kBz,
+ const char* evfoldname):TObject(),
+ fEvFolderName(evfoldname),
+ fBz(kBz),
+ fColl(kcoll),
+ fSelAndSmear(kTRUE),
+ fDBfileName(""),
+ fTrack(),
+ fCovTree(0),
+ fDBgrid(0),
+ fDBgridPi(),
+ fDBgridKa(),
+ fDBgridPr(),
+ fDBgridEl(),
+ fDBgridMu(),
+ fEff(0),
+ fEffPi(),
+ fEffKa(),
+ fEffPr(),
+ fEffEl(),
+ fEffMu(),
+ fPulls(0),
+ fRegPar(0),
+ fRegParPi(),
+ fRegParKa(),
+ fRegParPr(),
+ fRegParEl(),
+ fRegParMu(),
+ fdEdxMean(0),
+ fdEdxMeanPi(),
+ fdEdxMeanKa(),
+ fdEdxMeanPr(),
+ fdEdxMeanEl(),
+ fdEdxMeanMu(),
+ fdEdxRMS(0),
+ fdEdxRMSPi(),
+ fdEdxRMSKa(),
+ fdEdxRMSPr(),
+ fdEdxRMSEl(),
+ fdEdxRMSMu()
{
//-----------------------------------------------------------------------------
// This is the class conctructor
//-----------------------------------------------------------------------------
- fBz = Bz; // value of the z component of L3 field (Tesla)
- fColl = coll; // collision code (0: PbPb6000)
- 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) {
- cerr<<"AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n";
- cerr<<" Available: 0 -> PbPb6000"<<endl;
+ if(fColl!=0 && fColl!=1) {
+ Fatal("AliTPCtrackerParam","AliTPCtrackerParam::AliTPCtrackerParam: Invalid collision!\n Available: 0 -> PbPb6000\n 1 -> pp");
}
fDBfileName = gSystem->Getenv("ALICE_ROOT");
fDBfileName.Append("/TPC/CovMatrixDB_");
+ //fDBfileName = "CovMatrixDB_";
if(fColl==0) fDBfileName.Append("PbPb6000");
+ if(fColl==1) fDBfileName.Append("pp");
if(fBz==0.4) fDBfileName.Append("_B0.4T.root");
+ // use same DB for 0.4 and 0.5 T; for 0.5 T, correction done in CookTrack()
+ if(fBz==0.5) fDBfileName.Append("_B0.4T.root");
}
//-----------------------------------------------------------------------------
-AliTPCtrackerParam::~AliTPCtrackerParam()
-{}
-//-----------------------------------------------------------------------------
-
-Int_t AliTPCtrackerParam::BuildTPCtracks(const TFile *inp, TFile *out, Int_t n)
+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();
+}
//-----------------------------------------------------------------------------
-// This function creates the TPC parameterized tracks
+Int_t AliTPCtrackerParam::Init() {
+//-----------------------------------------------------------------------------
+// This function reads the DB from file
//-----------------------------------------------------------------------------
-
- TFile *fileDB=0;
- TTree *covTreePi[50];
- TTree *covTreeKa[50];
- TTree *covTreeEl[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;
- }
- // 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 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());
+ 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");
- } else cerr<<"\n ! Creating ALL TRUE tracks at TPC 1st hit !\n\n";
-
- TFile *infile=(TFile*)inp;
- infile->cd();
-
- // Get gAlice object from file
- if(!(gAlice=(AliRun*)infile->Get("gAlice"))) {
- cerr<<"gAlice has not been found on galice.root !\n";
+ 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;
}
- // Check if value in the galice file is equal to selected one (fBz)
- AliMagFCM *fiel = (AliMagFCM*)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;
+ 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
+//-----------------------------------------------------------------------------
+
+
+ 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;
+ }
+
+ // Get gAlice object from file
+ if(!gAlice) rl->LoadgAlice();
+ //rl->LoadHeader();
+ rl->LoadKinematics();
+ tpcloader->LoadHits("read");
+
+ if(!(gAlice=rl->GetAliRun())) {
+ printf("Cannot get gAlice from Run Loader !\n");
+ return 1;
}
// Get TPC detector
- AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
- Int_t ver = TPC->IsVersion();
- cerr<<"+++ TPC version "<<ver<<" has been found !\n";
- AliTPCParam *digp=(AliTPCParam*)infile->Get("75x40_100x60");
+ AliTPC *tpc=(AliTPC*)gAlice->GetDetector("TPC");
+
+ rl->CdGAFile();
+ AliTPCParam *digp=(AliTPCParam*)gDirectory->Get("75x40_100x60");
if(digp){
delete digp;
digp = new AliTPCParamSR();
}
- else digp=(AliTPCParam*)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;
- AliTPChit *tpcHit=0;
- AliTPCtrack *tpctrack=0;
- Double_t hPx,hPy,hPz,hPt,hEta,xg,yg,zg,xl,yl,zl;
- Double_t alpha;
- Float_t cosAlpha,sinAlpha;
+ TParticle *part=0;
+ AliTPCseedGeant *seed=0;
+ AliTPCtrack *tpctrack=0;
+ Double_t sPt,sEta;
Int_t bin,label,pdg,charge;
- Int_t tracks=0;
- Int_t nParticles,nTracks,arrentr;
- Char_t tname[100];
+ Int_t tracks;
+ Int_t nParticles,nSeeds,arrentr;
//Int_t nSel=0,nAcc=0;
- // loop over first n events in file
- for(Int_t evt=0; evt<n; evt++){
- cerr<<"+++\n+++ Processing event "<<evt<<"\n+++\n";
+ 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;
- // 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);
-
- // get the particles stack
- nParticles = gAlice->GetEvent(evt);
- Bool_t * done = new Bool_t[nParticles];
- Int_t * pdgCodes = new Int_t[nParticles];
-
- // loop on particles and store pdg codes
- for(Int_t l=0; l<nParticles; l++) {
- Part = gAlice->Particle(l);
- pdgCodes[l] = Part->GetPdgCode();
- done[l] = kFALSE;
- }
+ // array for TPC tracks
+ TObjArray tArray(20000);
+
+ // array for TPC seeds with geant info
+ TObjArray sArray(20000);
+
+ // get the particles stack
+ nParticles = (Int_t)gAlice->GetEvent(evt);
+
+ Bool_t *done = new Bool_t[nParticles];
+ Int_t *pdgCodes = new Int_t[nParticles];
+
+ // loop on particles and store pdg codes
+ for(Int_t l=0; l<nParticles; l++) {
+ part = (TParticle*)gAlice->GetMCApp()->Particle(l);
+ pdgCodes[l] = part->GetPdgCode();
+ done[l] = kFALSE;
+ }
+
+ 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=gAlice->TreeH();
- nTracks=(Int_t)TH->GetEntries();
- cerr<<"+++\n+++ Number of particles in event "<<evt<<": "<<nParticles<<
- "\n+++\n+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
- "\n+++\n\n";
-
- // loop over entries in TreeH
- 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);
- // Get FirstHit
- tpcHit=(AliTPChit*)TPC->FirstHit(-1);
- for( ; tpcHit; tpcHit=(AliTPChit*)TPC->NextHit() ) {
- if(tpcHit->fQ !=0.) continue;
- // Get particle momentum at hit
- hPx=tpcHit->X(); hPy=tpcHit->Y(); hPz=tpcHit->Z();
- hPt=TMath::Sqrt(hPx*hPx+hPy*hPy);
- // reject hits with Pt<mag*0.45 GeV/c
- if(hPt<(fBz*0.45)) continue;
-
- // Get track label
- label=tpcHit->Track();
- // 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);
-
- if((tpcHit=(AliTPChit*)TPC->NextHit())==0) break;
- if(tpcHit->fQ != 0.) continue;
- // Get global coordinates of hit
- xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
- if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
-
- // Get TPC sector, Alpha angle and local coordinates
- // cerr<<"Sector "<<tpcHit->fSector<<endl;
- digp->AdjustCosSin(tpcHit->fSector,cosAlpha,sinAlpha);
- alpha = TMath::ATan2(sinAlpha,cosAlpha);
- xl = xg*cosAlpha + yg*sinAlpha;
- yl =-xg*sinAlpha + yg*cosAlpha;
- zl = zg;
- //printf("Alpha %f xl %f yl %f zl %f\n",Alpha,xl,yl,zl);
-
- // reject tracks which are not in the TPC acceptance
- if(TMath::Abs(zl+(244.-xl)*hPz/hPt)>252.) continue;
-
- hEta = -TMath::Log(TMath::Tan(0.25*TMath::Pi()-0.5*TMath::ATan(hPz/hPt)));
-
- // Apply selection according to TPC efficiency
- //if(TMath::Abs(pdg)==211) nAcc++;
- if(fSelAndSmear && !SelectedTrack(hPt,hEta)) continue;
- //if(TMath::Abs(pdg)==211) nSel++;
-
- // create AliTPCtrack object
- BuildTrack(alpha,xl,yl,zl,hPx,hPy,hPz,hPt,charge);
-
- if(fSelAndSmear) {
- bin = fDBgrid->GetBin(hPt,hEta);
- switch (pdg) {
- case 211:
- fCovTree = covTreePi[bin];
- break;
- case 321:
- fCovTree = covTreeKa[bin];
- break;
- case 2212:
- fCovTree = covTreePi[bin];
- break;
- case 11:
- fCovTree = covTreeEl[bin];
- break;
- case 13:
- fCovTree = covTreePi[bin];
- break;
- }
- // deal with covariance matrix and smearing of parameters
- CookTrack(hPt,hEta);
-
- // assign the track a dE/dx and make a rough PID
- CookdEdx(hPt,hEta);
- }
-
- // 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++;
+ TTree *th = tpcloader->TreeH();
+ MakeSeedsFromHits(tpc,th,sArray);
+ } else {
+ // Get TreeTR with track references
+ rl->LoadTrackRefs();
+ TTree *ttr = rl->TreeTR();
+ if (!ttr) return -3;
+ MakeSeedsFromRefs(ttr,sArray);
+ }
+
+ nSeeds = sArray.GetEntries();
+ printf("+++ Number of seeds: %d\n",nSeeds);
+
+ // loop over entries in sArray
+ for(Int_t l=0; 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;
}
-
- } // loop over entries in TreeH
+ // deal with covariance matrix and smearing of parameters
+ CookTrack(sPt,sEta);
- // 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;
-
- printf("\n\n+++\n+++ Number of TPC tracks: %d\n+++\n",tracks);
- //cerr<<"Average Eff: "<<(Float_t)nSel/nAcc<<endl;
-
- } // loop on events
-
- if(fileDB) fileDB->Close();
-
+
+ // 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();
+
return 0;
}
//-----------------------------------------------------------------------------
gStyle->SetOptStat(0);
gStyle->SetOptFit(10001);
- Char_t *part="PIONS";
+ const char *part="PIONS";
Double_t ymax=500.;
+ /*
// create a chain with compared tracks
- TChain cmptrkchain("cmptrktree");
- cmptrkchain.Add("CovMatrix_AllEvts_1.root");
- cmptrkchain.Add("CovMatrix_AllEvts_2.root");
- cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ TChain *cmptrkchain = new ("cmptrktree");
+ cmptrkchain.Add("CovMatrix_AllEvts.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ */
+
+
+ TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
COMPTRACK cmptrk;
- cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
- Int_t entries = (Int_t)cmptrkchain.GetEntries();
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
cerr<<" Number of entries: "<<entries<<endl;
- InitializeKineGrid("DB","points");
- InitializeKineGrid("dEdx","bins");
+ InitializeKineGrid("DB");
+ InitializeKineGrid("dEdx");
switch(pdg) {
case 211:
part = "ELECTRONS";
ymax = 100.;
break;
+ case 13:
+ part = "MUONS";
+ ymax = 100.;
+ break;
}
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.;
}
// loop on chain entries for the mean
for(Int_t l=0; l<entries; l++) {
- cmptrkchain.GetEvent(l);
+ cmptrktree->GetEvent(l);
if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
p[bin] += cmptrk.p;
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
// loop on chain entries for the sigma
for(Int_t l=0; l<entries; l++) {
- cmptrkchain.GetEvent(l);
+ cmptrktree->GetEvent(l);
if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
if(cmptrk.p<1. && TMath::Abs(cmptrk.p-p[bin])>0.025) continue;
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<knTotBins; i++) {
+ fdEdxMeanMu.SetParam(i,mean[i]);
+ fdEdxRMSMu.SetParam(i,sigma[i]);
+ }
+ break;
}
// write results to file
// Output file is pulls.root.
//-----------------------------------------------------------------------------
+ /*
// create a chain with compared tracks
- TChain cmptrkchain("cmptrktree");
- cmptrkchain.Add("CovMatrix_AllEvts_1.root");
- cmptrkchain.Add("CovMatrix_AllEvts_2.root");
- cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ TChain *cmptrkchain = new ("cmptrktree");
+ cmptrkchain.Add("CovMatrix_AllEvts.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ */
+
+
+ TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
COMPTRACK cmptrk;
- cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
- Int_t entries = (Int_t)cmptrkchain.GetEntries();
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
cerr<<" Number of entries: "<<entries<<endl;
Int_t thisPdg=0;
TH1F *hDum = new TH1F("name","title",100,-7.,7.);
TF1 *g = new TF1("g","gaus");
- InitializeKineGrid("pulls","bins");
- InitializeKineGrid("DB","points");
+ InitializeKineGrid("pulls");
+ InitializeKineGrid("DB");
- // loop on the particles Pi,Ka,El
- for(Int_t part=0; part<3; part++) {
+ // loop on the particles Pi,Ka,Pr,El,Mu
+ for(Int_t part=0; part<5; part++) {
switch (part) {
case 0: // pions
thisPdg=321;
cerr<<" Processing kaons ...\n";
break;
-
- case 2: // electrons
+ case 2: // protons
+ thisPdg=2212;
+ cerr<<" Processing protons ...\n";
+ break;
+ case 3: // electrons
thisPdg=11;
cerr<<" Processing electrons ...\n";
break;
+ case 4: // muons
+ thisPdg=13;
+ cerr<<" Processing muons ...\n";
+ break;
}
SetParticle(thisPdg);
new(&pulls[i]) AliTPCkineGrid(*(fPulls+i));
}
nTotBins = fDBgrid->GetTotBins();
-
+ 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
for(Int_t i=0; i<entries; i++) {
- cmptrkchain.GetEvent(i);
+ cmptrktree->GetEvent(i);
if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
// fill histograms with the pulls
- bin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
- 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));
+ 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));
} // loop on chain entries
// compute the sigma of the distributions
for(Int_t i=0; i<nTotBins; i++) {
- if(hPulls0_[i].GetEntries()>1000) {
- 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()>1000) {
- 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()>1000) {
- 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()>1000) {
- 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()>1000) {
- 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
new(&fPullsKa[i]) AliTPCkineGrid(pulls[i]);
}
break;
- case 2: // electrons
+ case 2: // protons
+ for(Int_t i=0;i<5;i++) {
+ fPullsPr[i].~AliTPCkineGrid();
+ new(&fPullsPr[i]) AliTPCkineGrid(pulls[i]);
+ }
+ break;
+ case 3: // electrons
for(Int_t i=0;i<5;i++) {
fPullsEl[i].~AliTPCkineGrid();
new(&fPullsEl[i]) AliTPCkineGrid(pulls[i]);
}
break;
+ case 4: // muons
+ for(Int_t i=0;i<5;i++) {
+ fPullsMu[i].~AliTPCkineGrid();
+ new(&fPullsMu[i]) AliTPCkineGrid(pulls[i]);
+ //cerr<<" mu pulls "<<i<<" "<<fPullsMu[i].GetParam(0)<<endl;
+ }
+ 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
return;
}
//-----------------------------------------------------------------------------
-void AliTPCtrackerParam::BuildTrack(Double_t alpha,
- Double_t x,Double_t y,Double_t z,
- Double_t px,Double_t py,
- Double_t pz,Double_t pt,
- Int_t ch) {
+void AliTPCtrackerParam::AnalyzeResolutions(Int_t pdg) {
+//-----------------------------------------------------------------------------
+// This function computes the resolutions:
+// - dy
+// - dC
+// - dPt/Pt
+// as a function of Pt
+// Input file is CovMatrix_AllEvts.root.
+//-----------------------------------------------------------------------------
+
+ /*
+ // create a chain with compared tracks
+ TChain *cmptrkchain = new ("cmptrktree");
+ cmptrkchain.Add("CovMatrix_AllEvts.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ */
+
+
+ TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
+
+ COMPTRACK cmptrk;
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
+ cerr<<" Number of entries: "<<entries<<endl;
+
+
+ Int_t bin = 0;
+
+ InitializeKineGrid("DB");
+ InitializeKineGrid("eff");
+
+ SetParticle(pdg);
+
+ 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);
+
+
+ TH1F *hDumP0 = new TH1F("nameP0","dy",100,-0.3,0.3);
+ TH1F *hDumP4 = new TH1F("nameP4","dC",100,-0.0005,0.0005);
+ TH1F *hDumPt = new TH1F("namePt","dp_{T}/p_{T}",100,-0.5,0.5);
+
+ TF1 *g = new TF1("g","gaus");
+
+ // create histograms for the all the bins
+ TH1F *hP0=NULL;
+ TH1F *hP4=NULL;
+ TH1F *hPt=NULL;
+
+ hP0 = new TH1F[knPtBins];
+ hP4 = new TH1F[knPtBins];
+ hPt = new TH1F[knPtBins];
+
+ for(Int_t i=0; i<knPtBins; i++) {
+ hP0[i] = *hDumP0;
+ hP4[i] = *hDumP4;
+ hPt[i] = *hDumPt;
+ }
+
+ // loop on chain entries
+ for(Int_t i=0; i<entries; i++) {
+ cmptrktree->GetEvent(i);
+ if(TMath::Abs(cmptrk.pdg)!=pdg) continue;
+ // fill histograms with the residuals
+ bin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
+ //cerr<<" pt "<<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);
+ } // loop on chain entries
+
+
+ TCanvas *cP0res = new TCanvas("cP0res","cP0res",0,0,1200,700);
+ cP0res->Divide(5,2);
+ TCanvas *cP4res = new TCanvas("cP4res","cP4res",0,0,1200,700);
+ cP4res->Divide(5,2);
+ TCanvas *cPtres = new TCanvas("cPtres","cPtres",0,0,1200,700);
+ cPtres->Divide(5,2);
+
+ // Draw histograms
+ for(Int_t i=0; 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<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");
+ dP4[i] = g->GetParameter(2);
+ } else dP4[i] = 0.;
+ if(hPt[i].GetEntries()>10) {
+ g->SetRange(-3.*hPt[i].GetRMS(),3.*hPt[i].GetRMS());
+ hPt[i].Fit("g","R,Q,N");
+ dPtToPt[i] = 100.*g->GetParameter(2);
+ } else dPtToPt[i] = 0.;
+ } // loop on bins
+
+
+ TGraph *grdP0 = new TGraph(knPtBins,pt,dP0);
+ TGraph *grdP4 = new TGraph(knPtBins,pt,dP4);
+ TGraph *grdPtToPt = new TGraph(knPtBins,pt,dPtToPt);
+
+ grdP0->SetMarkerStyle(20); grdP0->SetMarkerColor(2); grdP0->SetMarkerSize(1.5);
+ grdP4->SetMarkerStyle(21); grdP4->SetMarkerColor(3); grdP4->SetMarkerSize(1.5);
+ grdPtToPt->SetMarkerStyle(22); grdPtToPt->SetMarkerColor(4); grdPtToPt->SetMarkerSize(1.5);
+
+ // Plot Results
+ gStyle->SetOptStat(0);
+ TCanvas *c1 = new TCanvas("c1","dP0",0,0,900,700);
+ c1->SetLogx();
+ c1->SetGridx();
+ c1->SetGridy();
+
+ TH2F *frame1 = new TH2F("frame1","y resolution VS p_{T} in TPC",2,0.1,30,2,0,0.1);
+ frame1->SetXTitle("p_{T} [GeV/c]");
+ frame1->SetYTitle("#sigma(y) [cm]");
+ frame1->Draw();
+ grdP0->Draw("P");
+
+
+ TCanvas *c2 = new TCanvas("c2","dP4",0,0,900,700);
+ c2->SetLogx();
+ c2->SetGridx();
+ c2->SetGridy();
+
+ TH2F *frame2 = new TH2F("frame2","C resolution VS p_{T} in TPC",2,0.1,30,2,0,0.0001);
+ frame2->SetXTitle("p_{T} [GeV/c]");
+ frame2->SetYTitle("#sigma(C) [1/cm]");
+ frame2->Draw();
+ grdP4->Draw("P");
+
+ TCanvas *c3 = new TCanvas("c3","dPtToPt",0,0,900,700);
+ c3->SetLogx();
+ c3->SetLogy();
+ c3->SetGridx();
+ c3->SetGridy();
+
+ TH2F *frame3 = new TH2F("frame3","Relative p_{T} resolution VS p_{T} in TPC",2,0.1,30,2,0.1,30.);
+ frame3->SetXTitle("p_{T} [GeV/c]");
+ frame3->SetYTitle("dp_{T}/p_{T} (%)");
+ frame3->Draw();
+ grdPtToPt->Draw("P");
+
+
+ delete [] dP0;
+ delete [] dP4;
+ delete [] dPtToPt;
+ delete [] pt;
+
+
+ delete [] hP0;
+ delete [] hP4;
+ delete [] hPt;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::BuildTrack(AliTPCseedGeant *s,Int_t ch) {
//-----------------------------------------------------------------------------
// This function uses GEANT info to set true track parameters
//-----------------------------------------------------------------------------
- Double_t xref = x;
+ Double_t xref = s->GetXL();
Double_t xx[5],cc[15];
- cc[0]=cc[2]=cc[5]=cc[9]=cc[14]=10.;
+ cc[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 = pt*100./0.299792458/bfield.Z();
+ Double_t rho = s->GetPt()*100./0.299792458/TMath::Abs(bfield.Z());
// center of track projection in local reference frame
- TVector3 hmom,hpos;
+ TVector3 sMom,sPos;
- // position (local) and momentum (local) at the hit
+ // position (local) and momentum (local) at the seed
// in the bending plane (z=0)
- hpos.SetXYZ(x,y,0.);
- hmom.SetXYZ(px*TMath::Cos(alpha)+py*TMath::Sin(alpha),-px*TMath::Sin(alpha)+py*TMath::Cos(alpha),0.);
- TVector3 vrho = hmom.Cross(bfield);
+ sPos.SetXYZ(s->GetXL(),s->GetYL(),0.);
+ sMom.SetXYZ(s->GetPx()*TMath::Cos(s->GetAlpha())+s->GetPy()*TMath::Sin(s->GetAlpha()),-s->GetPx()*TMath::Sin(s->GetAlpha())+s->GetPy()*TMath::Cos(s->GetAlpha()),0.);
+ TVector3 vrho = sMom.Cross(bfield);
vrho *= ch;
vrho.SetMag(rho);
- TVector3 vcenter = hpos+vrho;
+ TVector3 vcenter = sPos+vrho;
Double_t x0 = vcenter.X();
// fAlpha = Alpha Rotation angle the local (TPC sector)
// fP0 = YL Y-coordinate of a track
// fP1 = ZG Z-coordinate of a track
- // fP2 = C*x0 x0 is center x in rotated frame
+ // fP2 = sin(phi) sine of the (local) azimuthal angle
// fP3 = Tgl tangent of the track momentum dip angle
// fP4 = C track curvature
- xx[0] = y;
- xx[1] = z;
- xx[3] = pz/pt;
- xx[4] = -ch/rho;
- xx[2] = xx[4]*x0;
+ xx[0] = s->GetYL();
+ xx[1] = s->GetZL();
+ xx[2] = ch/rho*(xref-x0);
+ xx[3] = s->GetPz()/s->GetPt();
+ xx[4] = ch/rho;
// create the object AliTPCtrack
- AliTPCtrack track(0,xx,cc,xref,alpha);
+ AliTPCtrack track(xref,s->GetAlpha(),xx,cc,0);
new(&fTrack) AliTPCtrack(track);
return;
}
//-----------------------------------------------------------------------------
void AliTPCtrackerParam::CompareTPCtracks(
+ Int_t nEvents,
const Char_t* galiceName,
const Char_t* trkGeaName,
const Char_t* trkKalName,
const Char_t* covmatName,
- const Char_t* tpceffName) const {
+ const Char_t* tpceffasciiName,
+ const Char_t* tpceffrootName) {
//-----------------------------------------------------------------------------
// This function compares tracks from TPC Kalman Filter V2 with
-// true tracks at TPC 1st hit. It gives:
+// geant tracks at TPC 1st hit. It gives:
// - a tree with Kalman cov. matrix elements, resolutions, dEdx
// - the efficiencies as a function of particle type, pT, eta
//-----------------------------------------------------------------------------
TFile *geaFile = TFile::Open(trkGeaName);
TFile *galiceFile = TFile::Open(galiceName);
- // particles from TreeK
- AliRun *gAlice = (AliRun*)galiceFile->Get("gAlice");
- const Int_t nparticles = gAlice->GetEvent(0);
+ // get the AliRun object
+ AliRun *lAlice = (AliRun*)galiceFile->Get("gAlice");
- Int_t * kalLab = new Int_t[nparticles];
- for(Int_t i=0; i<nparticles; i++) kalLab[i] = -1;
-
- // tracks from Kalman
- TTree *kaltree=(TTree*)kalFile->Get("TreeT_TPC_0");
- AliTPCtrack *kaltrack=new AliTPCtrack;
- kaltree->SetBranchAddress("tracks",&kaltrack);
- Int_t kalEntries = (Int_t)kaltree->GetEntries();
-
- // tracks from 1st hit
- TTree *geatree=(TTree*)geaFile->Get("TreeT_TPC_0");
- AliTPCtrack *geatrack=new AliTPCtrack;
- geatree->SetBranchAddress("tracks",&geatrack);
- Int_t geaEntries = (Int_t)geatree->GetEntries();
-
- cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<<kalEntries<<endl<<"+++ TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
-
- // set pointers for TPC tracks:
- // the entry number of the track labelled l is stored in kalLab[l]
- Int_t fake=0,mult=0;
- for (Int_t j=0; j<kalEntries; j++) {
- kaltree->GetEvent(j);
- if(kaltrack->GetLabel()>=0) {
- if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
- kalLab[kaltrack->GetLabel()] = j;
- }
- else fake++;
- }
-
- cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
- cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
+ // create the tree for comparison results
+ COMPTRACK cmptrk;
+ TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
+ cmptrktree->Branch("comptracks",&cmptrk,"pdg/I:bin:r/D:p:pt:cosl:eta:dpt:dP0:dP1:dP2:dP3:dP4:c00:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44:dEdx");
- TParticle *Part;
+ InitializeKineGrid("eff");
+ InitializeKineGrid("DB");
+ Int_t effBins = fEffPi.GetTotPoints();
+ Int_t effBinsPt = fEffPi.GetPointsPt();
+ Double_t *pt = new Double_t[effBinsPt];
+ fEffPi.GetArrayPt(pt);
+
+ TParticle *part;
+ Double_t ptgener;
+ Bool_t usethis;
Int_t label;
- Double_t cc[15],dAlpha;
- // --- [Pt,eta] binning ---
- //
- // |eta|<0.3 0.3<=|eta|<0.6 0.6<=|eta|
- // Pt<0.3 0 1 2
- // 0.3<=Pt<0.5 3 4 5
- // 0.5<=Pt<1.0 6 7 8
- // 1.0<=Pt<1.5 9 10 11
- // 1.5<=Pt<3.0 12 13 14
- // 3.0<=Pt<5.0 15 16 17
- // 5.0<=Pt<7.0 18 19 20
- // 7.0<=Pt<15. 21 22 23
- // 15.<=Pt 24 25 26
- //
+ Double_t dAlpha;
Int_t pi=0,ka=0,mu=0,el=0,pr=0;
- Int_t geaPi[27],geaKa[27],geaPr[27],geaEl[27],geaMu[27];
- Int_t kalPi[27],kalKa[27],kalPr[27],kalEl[27],kalMu[27];
- Float_t effPi[27],effKa[27],effPr[27],effEl[27],effMu[27];
+ Int_t *geaPi = new Int_t[effBins];
+ Int_t *geaKa = new Int_t[effBins];
+ Int_t *geaPr = new Int_t[effBins];
+ Int_t *geaEl = new Int_t[effBins];
+ Int_t *geaMu = new Int_t[effBins];
+ Int_t *kalPi = new Int_t[effBins];
+ Int_t *kalKa = new Int_t[effBins];
+ Int_t *kalPr = new Int_t[effBins];
+ Int_t *kalEl = new Int_t[effBins];
+ Int_t *kalMu = new Int_t[effBins];
+ Float_t *effPi = new Float_t[effBins];
+ Float_t *effKa = new Float_t[effBins];
+ Float_t *effPr = new Float_t[effBins];
+ Float_t *effEl = new Float_t[effBins];
+ Float_t *effMu = new Float_t[effBins];
Int_t bin;
- for(Int_t j=0; j<27; j++) {
+ for(Int_t j=0; j<effBins; j++) {
geaPi[j]=geaKa[j]=geaPr[j]=geaEl[j]=geaMu[j]=0;
kalPi[j]=kalKa[j]=kalPr[j]=kalEl[j]=kalMu[j]=0;
effPi[j]=effKa[j]=effPr[j]=effEl[j]=effMu[j]=-1.;
}
- COMPTRACK cmptrk;
- // create the tree for comparison results
- TTree *cmptrktree = new TTree("cmptrktree","results of track comparison");
- cmptrktree->Branch("comptracks",&cmptrk,"pdg/I:bin:r/D:p:pt:cosl:eta:dpt:dP0:dP1:dP2:dP3:dP4:c00:c10:c11:c20:c21:c22:c30:c31:c32:c33:c40:c41:c42:c43:c44:dEdx");
-
- cerr<<"Doing track comparison...\n";
- // loop on tracks at 1st hit
- for(Int_t j=0; j<geaEntries; j++) {
- geatree->GetEvent(j);
+ Char_t tname[100];
+
+ // loop on events in file
+ 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 knparticles = lAlice->GetEvent(evt);
+
+ Int_t *kalLab = new Int_t[knparticles];
+ for(Int_t i=0; i<knparticles; i++) kalLab[i] = -1;
+
+
+ // tracks from Kalman
+ TTree *kaltree=(TTree*)kalFile->Get(tname);
+ if(!kaltree) continue;
+ AliTPCtrack *kaltrack=new AliTPCtrack;
+ kaltree->SetBranchAddress("tracks",&kaltrack);
+ Int_t kalEntries = (Int_t)kaltree->GetEntries();
+
+ // tracks from 1st hit
+ TTree *geatree=(TTree*)geaFile->Get(tname);
+ if(!geatree) continue;
+ AliTPCtrack *geatrack=new AliTPCtrack;
+ geatree->SetBranchAddress("tracks",&geatrack);
+ Int_t geaEntries = (Int_t)geatree->GetEntries();
- label = geatrack->GetLabel();
- Part = (TParticle*)gAlice->Particle(label);
- cmptrk.pdg = Part->GetPdgCode();
- cmptrk.eta = Part->Eta();
- cmptrk.r = TMath::Sqrt(Part->Vx()*Part->Vx()+Part->Vy()*Part->Vy());
-
- cmptrk.pt = 1/TMath::Abs(geatrack->Get1Pt());
- cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
- cmptrk.p = cmptrk.pt/cmptrk.cosl;
-
-
- bin = GetBin(cmptrk.pt,cmptrk.eta);
- cmptrk.bin = bin;
- if(abs(cmptrk.pdg)==211) geaPi[bin]++;
- if(abs(cmptrk.pdg)==321) geaKa[bin]++;
- if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
- if(abs(cmptrk.pdg)==11) geaEl[bin]++;
- if(abs(cmptrk.pdg)==13) geaMu[bin]++;
-
-
- // check if there is the corresponding track in the TPC kalman and get it
- if(kalLab[label]==-1) continue;
- kaltree->GetEvent(kalLab[label]);
-
- // and go on only if it has xref = 84.57 cm (inner pad row)
- if(kaltrack->GetX()>90.) continue;
-
- if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; }
- if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; }
- if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
- if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; }
- if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; }
-
- kaltrack->PropagateTo(geatrack->GetX());
-
- cmptrk.dEdx = kaltrack->GetdEdx();
-
- // compute errors on parameters
- dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
- if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
-
- cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
- cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
- cmptrk.dP2 = kaltrack->GetEta()-geatrack->GetEta();
- cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
- cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
- cmptrk.dpt = 1/kaltrack->Get1Pt()-1/geatrack->Get1Pt();
-
- // get covariance matrix
- // beware: lines 3 and 4 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];
-
- // fill tree
- cmptrktree->Fill();
-
- } // loop on tracks at TPC 1st hit
+ cerr<<"+++\n+++ Number of tracks: TPC Kalman = "<<kalEntries<<endl<<"+++ TPC 1st hit = "<<geaEntries<<endl<<"+++\n";
+
+ // set pointers for TPC tracks:
+ // the entry number of the track labelled l is stored in kalLab[l]
+ Int_t fake=0,mult=0;
+ for (Int_t j=0; j<kalEntries; j++) {
+ kaltree->GetEvent(j);
+ if(kaltrack->GetLabel()>=0) {
+ if(kalLab[kaltrack->GetLabel()]!=-1) mult++;
+ kalLab[kaltrack->GetLabel()] = j;
+ }
+ else fake++;
+ }
+
+ cerr<<"+++ Number of fake tracks in TPC Kalman: "<<fake<<endl;
+ cerr<<"+++ Number of multiply found tracks in TPC Kalman: "<<mult<<endl;
+
+ /*
+ // Read the labels of the seeds
+ char sname[100];
+ Int_t sLabel,ncol;
+ 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) {
+ ncol = fscanf(seedFile,"%d",&sLabel);
+ if(ncol<1) break;
+ if(sLabel>0) hasSeed[sLabel]=kTRUE;
+ }
+ fclose(seedFile);
+ */
+
+ cerr<<"Doing track comparison...\n";
+ // loop on tracks at 1st hit
+ for(Int_t j=0; j<geaEntries; j++) {
+ geatree->GetEvent(j);
+
+ label = geatrack->GetLabel();
+ part = (TParticle*)lAlice->GetMCApp()->Particle(label);
+
+ // use only injected tracks with fixed values of 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(!usethis) continue;
+
+ // check if it has the seed
+ //if(!hasSeed[label]) continue;
+
+ /*
+ // check if track is entirely contained in a TPC sector
+ Bool_t out = kFALSE;
+ for(Int_t l=0; l<10; l++) {
+ Double_t x = 85. + (250.-85.)*(Double_t)l/9.;
+ //cerr<<"x "<<x<<" X "<<geatrack->GetX()<<endl;
+ Double_t y = geatrack->GetY() + (
+ TMath::Sqrt(1-(geatrack->GetC()*geatrack->GetX()-geatrack->GetEta())*
+ (geatrack->GetC()*geatrack->GetX()-geatrack->GetEta()))
+ -TMath::Sqrt(1-(geatrack->GetC()*x-geatrack->GetEta())*
+ (geatrack->GetC()*x-geatrack->GetEta()))
+ )/geatrack->GetC();
+ y = TMath::Abs(y);
+ //cerr<<"y "<<y<<" Y "<<geatrack->GetY()<<endl;
+ if(y > 0.8*x*TMath::Tan(TMath::Pi()/18)) { out = kTRUE; break; }
+ }
+ if(out) continue;
+ */
+
+ cmptrk.pdg = part->GetPdgCode();
+ cmptrk.eta = part->Eta();
+ cmptrk.r = TMath::Sqrt(part->Vx()*part->Vx()+part->Vy()*part->Vy());
+
+ cmptrk.pt = geatrack->Pt();
+ cmptrk.cosl = TMath::Cos(TMath::ATan(geatrack->GetTgl()));
+ cmptrk.p = cmptrk.pt/cmptrk.cosl;
+
+
+ bin = fDBgridPi.GetBin(cmptrk.pt,cmptrk.eta);
+ cmptrk.bin = bin;
+ if(abs(cmptrk.pdg)==211) geaPi[bin]++;
+ if(abs(cmptrk.pdg)==321) geaKa[bin]++;
+ if(abs(cmptrk.pdg)==2212) geaPr[bin]++;
+ if(abs(cmptrk.pdg)==11) geaEl[bin]++;
+ if(abs(cmptrk.pdg)==13) geaMu[bin]++;
+
+
+ // check if there is the corresponding track in the TPC kalman and get it
+ if(kalLab[label]==-1) continue;
+ kaltree->GetEvent(kalLab[label]);
+
+ // and go on only if it has xref = 84.57 cm (inner pad row)
+ if(kaltrack->GetX()>90.) continue;
+
+ if(abs(cmptrk.pdg)==211) { kalPi[bin]++; pi++; }
+ if(abs(cmptrk.pdg)==321) { kalKa[bin]++; ka++; }
+ if(abs(cmptrk.pdg)==2212) { kalPr[bin]++; pr++; }
+ if(abs(cmptrk.pdg)==11) { kalEl[bin]++; el++; }
+ if(abs(cmptrk.pdg)==13) { kalMu[bin]++; mu++; }
+
+ kaltrack->PropagateTo(geatrack->GetX(),1.e-9,0.);
+
+ cmptrk.dEdx = kaltrack->GetdEdx();
+
+ // compute errors on parameters
+ dAlpha = kaltrack->GetAlpha()-geatrack->GetAlpha();
+ if(TMath::Abs(dAlpha)>0.1) { cerr<<" ! WRONG SECTOR !\n"; continue; }
+
+ cmptrk.dP0 = kaltrack->GetY()-geatrack->GetY();
+ cmptrk.dP1 = kaltrack->GetZ()-geatrack->GetZ();
+ cmptrk.dP2 = kaltrack->GetSnp()-geatrack->GetSnp();
+ cmptrk.dP3 = kaltrack->GetTgl()-geatrack->GetTgl();
+ cmptrk.dP4 = kaltrack->GetC()-geatrack->GetC();
+ cmptrk.dpt = 1/kaltrack->GetSigned1Pt()-1/geatrack->GetSigned1Pt();
+
+ // get covariance matrix
+ // beware: lines 3 and 4 in the matrix are inverted!
+ //kaltrack->GetCovariance(cc);
+
+ cmptrk.c00 = kaltrack->GetSigmaY2();
+ cmptrk.c10 = kaltrack->GetSigmaZY();
+ cmptrk.c11 = kaltrack->GetSigmaZ2();
+ cmptrk.c20 = kaltrack->GetSigmaSnpY();
+ cmptrk.c21 = kaltrack->GetSigmaSnpY();
+ cmptrk.c22 = kaltrack->GetSigmaSnp2();
+ cmptrk.c30 = kaltrack->GetSigmaTglY();
+ cmptrk.c31 = kaltrack->GetSigmaTglZ();
+ cmptrk.c32 = kaltrack->GetSigmaTglSnp();
+ cmptrk.c33 = kaltrack->GetSigmaTgl2();
+ cmptrk.c40 = kaltrack->GetSigma1PtY();
+ cmptrk.c41 = kaltrack->GetSigma1PtZ();
+ cmptrk.c42 = kaltrack->GetSigma1PtSnp();
+ cmptrk.c43 = kaltrack->GetSigma1PtTgl();
+ cmptrk.c44 = kaltrack->GetSigma1Pt2();
+
+ // fill tree
+ cmptrktree->Fill();
+
+ } // loop on tracks at TPC 1st hit
+
+ delete [] kalLab;
+ //delete [] hasSeed;
+
+ } // end loop on events in file
+
cerr<<"+++\n+++ Number of compared tracks: "<<pi+ka+el+mu+pr<<endl;
cerr<<"+++ Pions: "<<pi<<", Kaons: "<<ka<<", Protons : "<<pr<<", Electrons: "<<el<<", Muons: "<<mu<<endl;
outfile->Close();
- // Write efficiencies to file
- FILE *effFile = fopen(tpceffName,"w");
+ // Write efficiencies to ascii file
+ FILE *effFile = fopen(tpceffasciiName,"w");
//fprintf(effFile,"%d\n",kalEntries);
- for(Int_t j=0; j<27; j++) {
- if(geaPi[j]>=5) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
- if(geaKa[j]>=5) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
- if(geaPr[j]>=5) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
- if(geaEl[j]>=5) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
- if(geaMu[j]>=5) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
+ for(Int_t j=0; j<effBins; j++) {
+ if(geaPi[j]>=100) effPi[j]=(Float_t)kalPi[j]/geaPi[j];
+ if(geaKa[j]>=100) effKa[j]=(Float_t)kalKa[j]/geaKa[j];
+ if(geaPr[j]>=100) effPr[j]=(Float_t)kalPr[j]/geaPr[j];
+ if(geaEl[j]>=500) effEl[j]=(Float_t)kalEl[j]/geaEl[j];
+ if(geaMu[j]>=100) effMu[j]=(Float_t)kalMu[j]/geaMu[j];
fprintf(effFile,"%f %f %f %f %f\n",effPi[j],effKa[j],effPr[j],effEl[j],effMu[j]);
}
- for(Int_t j=0; j<27; j++) {
+ for(Int_t j=0; j<effBins; j++) {
fprintf(effFile,"%d %d %d %d %d\n",geaPi[j],geaKa[j],geaPr[j],geaEl[j],geaMu[j]);
}
- for(Int_t j=0; j<27; j++) {
+ for(Int_t j=0; j<effBins; j++) {
fprintf(effFile,"%d %d %d %d %d\n",kalPi[j],kalKa[j],kalPr[j],kalEl[j],kalMu[j]);
}
fclose(effFile);
+ // Write efficiencies to root file
+ for(Int_t j=0; j<effBins; j++) {
+ fEffPi.SetParam(j,(Double_t)effPi[j]);
+ fEffKa.SetParam(j,(Double_t)effKa[j]);
+ fEffPr.SetParam(j,(Double_t)effPr[j]);
+ fEffEl.SetParam(j,(Double_t)effEl[j]);
+ fEffMu.SetParam(j,(Double_t)effMu[j]);
+ }
+ WriteEffs(tpceffrootName);
+
// delete AliRun object
- delete gAlice; gAlice=0;
+ delete lAlice; lAlice=0;
// close all input files
kalFile->Close();
geaFile->Close();
galiceFile->Close();
- delete [] kalLab;
+ delete [] pt;
+ delete [] effPi;
+ delete [] effKa;
+ delete [] effPr;
+ delete [] effEl;
+ delete [] effMu;
+ delete [] geaPi;
+ delete [] geaKa;
+ delete [] geaPr;
+ delete [] geaEl;
+ delete [] geaMu;
+ delete [] kalPi;
+ delete [] kalKa;
+ delete [] kalPr;
+ delete [] kalEl;
+ delete [] kalMu;
return;
}
//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) {
COVMATRIX covmat;
Double_t p,cosl;
- Double_t trkKine[2],trkRegPar[4];
+ Double_t trkKine[1],trkRegPar[3];
Double_t xref,alpha,xx[5],xxsm[5],cc[15];
Int_t treeEntries;
// 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;
- trkKine[1] = cosl;
// get covariance matrix from regularized matrix
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(0,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(0,l);
cc[0] = covmat.c00*RegFunc(trkKine,trkRegPar);
cc[1] = covmat.c10;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(1,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(1,l);
cc[2] = covmat.c11*RegFunc(trkKine,trkRegPar);
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(2,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(2,l);
cc[3] = covmat.c20*RegFunc(trkKine,trkRegPar);
cc[4] = covmat.c21;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(3,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(3,l);
cc[5] = covmat.c22*RegFunc(trkKine,trkRegPar);
cc[6] = covmat.c30;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(4,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(4,l);
cc[7] = covmat.c31*RegFunc(trkKine,trkRegPar);
cc[8] = covmat.c32;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(5,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(5,l);
cc[9] = covmat.c33*RegFunc(trkKine,trkRegPar);
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(6,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(6,l);
cc[10]= covmat.c40*RegFunc(trkKine,trkRegPar);
cc[11]= covmat.c41;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(7,l);
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(7,l);
cc[12]= covmat.c42*RegFunc(trkKine,trkRegPar);
cc[13]= covmat.c43;
- for(Int_t l=0;l<4;l++) trkRegPar[l] = (*fRegPar)(8,l);
- cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
-
- TMatrixD covMatSmear(5,5);
-
+ for(Int_t l=0;l<3;l++) trkRegPar[l] = (*fRegPar)(8,l);
+ cc[14]= covmat.c44*RegFunc(trkKine,trkRegPar);
+
+
+ TMatrixD covMatSmear(5,5);
covMatSmear = GetSmearingMatrix(cc,pt,eta);
// get track original parameters
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);
case 321:
title.Prepend("KAONS - ");
break;
+ case 2212:
+ title.Prepend("PROTONS - ");
+ break;
case 11:
title.Prepend("ELECTRONS - ");
break;
+ case 13:
+ title.Prepend("MUONS - ");
+ break;
}
switch (par) {
case 0:
return;
}
//-----------------------------------------------------------------------------
-Int_t AliTPCtrackerParam::GetBin(Double_t pt,Double_t eta) const {
-//-----------------------------------------------------------------------------
-// This function tells bin number in a grid (pT,eta)
-//-----------------------------------------------------------------------------
- if(TMath::Abs(eta)<0.3) {
- if(pt<0.3) return 0;
- if(pt>=0.3 && pt<0.5) return 3;
- if(pt>=0.5 && pt<1.) return 6;
- if(pt>=1. && pt<1.5) return 9;
- if(pt>=1.5 && pt<3.) return 12;
- if(pt>=3. && pt<5.) return 15;
- if(pt>=5. && pt<7.) return 18;
- if(pt>=7. && pt<15.) return 21;
- if(pt>=15.) return 24;
- }
- if(TMath::Abs(eta)>=0.3 && TMath::Abs(eta)<0.6) {
- if(pt<0.3) return 1;
- if(pt>=0.3 && pt<0.5) return 4;
- if(pt>=0.5 && pt<1.) return 7;
- if(pt>=1. && pt<1.5) return 10;
- if(pt>=1.5 && pt<3.) return 13;
- if(pt>=3. && pt<5.) return 16;
- if(pt>=5. && pt<7.) return 19;
- if(pt>=7. && pt<15.) return 22;
- if(pt>=15.) return 25;
- }
- if(TMath::Abs(eta)>=0.6) {
- if(pt<0.3) return 2;
- if(pt>=0.3 && pt<0.5) return 5;
- if(pt>=0.5 && pt<1.) return 8;
- if(pt>=1. && pt<1.5) return 11;
- if(pt>=1.5 && pt<3.) return 14;
- if(pt>=3. && pt<5.) return 17;
- if(pt>=5. && pt<7.) return 20;
- if(pt>=7. && pt<15.) return 23;
- if(pt>=15.) return 26;
- }
-
- return -1;
-
-}
-//-----------------------------------------------------------------------------
TMatrixD AliTPCtrackerParam::GetSmearingMatrix(Double_t* cc,Double_t pt,
Double_t eta) const {
//-----------------------------------------------------------------------------
covMat(4,3)=cc[13]; covMat(3,4)=covMat(4,3);
covMat(4,4)=cc[14];
+
TMatrixD stretchMat(5,5);
for(Int_t k=0;k<5;k++) {
for(Int_t l=0;l<5;l++) {
}
}
- for(Int_t i=0;i<5;i++) stretchMat(i,i) =
- TMath::Sqrt((fPulls+i)->GetValueAt(pt,eta));
+ for(Int_t i=0;i<5;i++) {
+ stretchMat(i,i) = (fPulls+i)->GetValueAt(pt,eta);
+ if(stretchMat(i,i)==0.) stretchMat(i,i) = 1.;
+ }
TMatrixD mat(stretchMat,TMatrixD::kMult,covMat);
TMatrixD covMatSmear(mat,TMatrixD::kMult,stretchMat);
return covMatSmear;
}
//-----------------------------------------------------------------------------
-void AliTPCtrackerParam::InitializeKineGrid(Option_t* which,Option_t* how) {
+void AliTPCtrackerParam::InitializeKineGrid(Option_t* which) {
//-----------------------------------------------------------------------------
// This function initializes ([pt,eta] points) the data members AliTPCkineGrid
// which = "DB" -> initialize fDBgrid... members
// "eff" -> initialize fEff... members
// "pulls" -> initialize fPulls... members
// "dEdx" -> initialize fdEdx... members
-// how = "points" -> initialize with the points of the grid
-// "bins" -> initialize with the bins of the grid
//-----------------------------------------------------------------------------
- const char *points = strstr(how,"points");
- const char *bins = strstr(how,"bins");
- 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");
- Int_t nEta=0,nPtPi=0,nPtKa=0,nPtPr=0,nPtEl=0,nPtMu=0;
+ Int_t nEta=0, nPt=0;
Double_t etaPoints[2] = {0.3,0.6};
Double_t etaBins[3] = {0.15,0.45,0.75};
- Double_t ptPoints8[8] = {0.3,0.5,1.,1.5,3.,5.,7.,15.};
- Double_t ptPoints5[5] = {0.3,0.5,1.,1.5,5.};
- Double_t ptBins9[9] = {0.244,0.390,0.676,1.190,2.36,4.,6.,10.,20.};
- Double_t ptBins6[6] = {0.244,0.390,0.676,1.190,2.36,10.};
-
- Double_t *eta=0,*ptPi=0,*ptKa=0,*ptPr=0,*ptEl=0,*ptMu=0;
-
- if(points) {
- nEta = 2;
- nPtPi = 8;
- nPtKa = 5;
- nPtPr = 4;
- nPtEl = 8;
- nPtMu = 2;
+
+ Double_t ptPoints[9] = {0.4,0.6,0.8,1.2,1.7,3.,5.,8.,15.};
+ Double_t ptBins[10] = {0.3,0.5,0.7,1.,1.5,2.,4.,6.,10.,20.};
+
+
+ Double_t *eta=0,*pt=0;
+
+ if(db) {
+ nEta = 2;
+ nPt = 9;
eta = etaPoints;
- ptPi = ptPoints8;
- ptKa = ptPoints5;
- ptPr = ptPoints8;
- ptEl = ptPoints8;
- ptMu = ptPoints8;
- }
- if(bins) {
- nEta = 3;
- nPtPi = 9;
- nPtKa = 6;
- nPtPr = 5;
- nPtEl = 9;
- nPtMu = 3;
+ pt = ptPoints;
+ } else {
+ nEta = 3;
+ nPt = 10;
eta = etaBins;
- ptPi = ptBins9;
- ptKa = ptBins6;
- ptPr = ptBins9;
- ptEl = ptBins9;
- ptMu = ptBins9;
+ pt = ptBins;
}
AliTPCkineGrid *dummy=0;
- if(DB) {
- dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
+ if(db) {
+ dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
new(&fDBgridPi) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
new(&fDBgridKa) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
+ new(&fDBgridPr) AliTPCkineGrid(*dummy);
new(&fDBgridEl) AliTPCkineGrid(*dummy);
+ new(&fDBgridMu) AliTPCkineGrid(*dummy);
delete dummy;
}
if(eff) {
- dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
+ dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
new(&fEffPi) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
new(&fEffKa) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtPr,nEta,ptPr,eta);
new(&fEffPr) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
new(&fEffEl) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtMu,nEta,ptMu,eta);
new(&fEffMu) AliTPCkineGrid(*dummy);
delete dummy;
}
if(pulls) {
- dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
+ dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
for(Int_t i=0;i<5;i++) new(&fPullsPi[i]) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
for(Int_t i=0;i<5;i++) new(&fPullsKa[i]) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
+ for(Int_t i=0;i<5;i++) new(&fPullsPr[i]) AliTPCkineGrid(*dummy);
for(Int_t i=0;i<5;i++) new(&fPullsEl[i]) AliTPCkineGrid(*dummy);
+ for(Int_t i=0;i<5;i++) new(&fPullsMu[i]) AliTPCkineGrid(*dummy);
delete dummy;
}
if(dEdx) {
- dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
+ dummy = new AliTPCkineGrid(nPt,nEta,pt,eta);
new(&fdEdxMeanPi) AliTPCkineGrid(*dummy);
new(&fdEdxRMSPi) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtKa,nEta,ptKa,eta);
new(&fdEdxMeanKa) AliTPCkineGrid(*dummy);
new(&fdEdxRMSKa) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtPi,nEta,ptPi,eta);
new(&fdEdxMeanPr) AliTPCkineGrid(*dummy);
new(&fdEdxRMSPr) AliTPCkineGrid(*dummy);
- delete dummy;
- dummy = new AliTPCkineGrid(nPtEl,nEta,ptEl,eta);
new(&fdEdxMeanEl) AliTPCkineGrid(*dummy);
new(&fdEdxRMSEl) AliTPCkineGrid(*dummy);
- delete dummy;
+ new(&fdEdxMeanMu) AliTPCkineGrid(*dummy);
+ new(&fdEdxRMSMu) AliTPCkineGrid(*dummy);
+ delete dummy;
}
return;
//-----------------------------------------------------------------------------
// define some file names
- Char_t *effFile ="CovMatrixDB_partial.root";
- Char_t *pullsFile="CovMatrixDB_partial.root";
- Char_t *regPiFile="CovMatrixDB_partial.root";
- Char_t *regKaFile="CovMatrixDB_partial.root";
- Char_t *regElFile="CovMatrixDB_partial.root";
- Char_t *dEdxPiFile="dEdxPi.root";
- Char_t *dEdxKaFile="dEdxKa.root";
- Char_t *dEdxPrFile="dEdxPr.root";
- Char_t *dEdxElFile="dEdxEl.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";
Char_t *cmFile3 ="CovMatrix_AllEvts_3.root";
+ */
// store the effieciencies
ReadEffs(effFile);
WriteRegParams(fDBfileName.Data(),211);
ReadRegParams(regKaFile,321);
WriteRegParams(fDBfileName.Data(),321);
+ ReadRegParams(regPrFile,2212);
+ WriteRegParams(fDBfileName.Data(),2212);
ReadRegParams(regElFile,11);
WriteRegParams(fDBfileName.Data(),11);
+ ReadRegParams(regMuFile,13);
+ WriteRegParams(fDBfileName.Data(),13);
// store the dEdx parameters
ReaddEdx(dEdxPiFile,211);
WritedEdx(fDBfileName.Data(),2212);
ReaddEdx(dEdxElFile,11);
WritedEdx(fDBfileName.Data(),11);
+ ReaddEdx(dEdxMuFile,13);
+ WritedEdx(fDBfileName.Data(),13);
//
// store the regularized covariance matrices
//
- InitializeKineGrid("DB","points");
+ InitializeKineGrid("DB");
- const Int_t nBinsPi = fDBgridPi.GetTotBins();
- const Int_t nBinsKa = fDBgridKa.GetTotBins();
- const Int_t nBinsEl = fDBgridEl.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 *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 *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<knBinsPr; i++) {
+ sprintf(hname,"CovTreePr_bin%d",i);
+ sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
+ 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<knBinsMu; i++) {
+ sprintf(hname,"CovTreeMu_bin%d",i);
+ sprintf(htitle,"Tree with cov matrix elements for bin %d",i);
+ 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);
+ }
+
+ /*
// create the chain with the compared tracks
- TChain cmptrkchain("cmptrktree");
+ TChain *cmptrktree = new TChain("cmptrktree");
cmptrkchain.Add(cmFile1);
cmptrkchain.Add(cmFile2);
cmptrkchain.Add(cmFile3);
+ */
+ TFile *infile = TFile::Open(cmFile);
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
+
COMPTRACK cmptrk;
- cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
- Int_t entries = (Int_t)cmptrkchain.GetEntries();
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
cerr<<" Number of entries: "<<entries<<endl;
Int_t trkPdg,trkBin;
- Double_t trkKine[2],trkRegPar[4];
- 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 * nPerBinEl = new Int_t[nBinsEl];
- for(Int_t k=0;k<nBinsEl;k++) nPerBinEl[k]=0;
+ Double_t trkKine[1],trkRegPar[3];
+ 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++) {
if(l % 10000 == 0) cerr<<"--- Processing track "<<l<<" of "<<entries<<" ---"<<endl;
// get the event
- cmptrkchain.GetEvent(l);
+ cmptrktree->GetEvent(l);
// get the pdg code
trkPdg = TMath::Abs(cmptrk.pdg);
- // use only pions, kaons, electrons
- if(trkPdg!=211 && trkPdg!=321 && trkPdg!=11) continue;
+ // use only pions, kaons, protons, electrons, muons
+ if(trkPdg!=211 && trkPdg!=321 && trkPdg!=2212 && trkPdg!=11 && trkPdg!=13) continue;
SetParticle(trkPdg);
trkBin = fDBgrid->GetBin(cmptrk.pt,cmptrk.eta);
//cerr<<cmptrk.pt<<" "<<cmptrk.eta<<" "<<trkBin<<endl;
if(trkPdg==211 && nPerBinPi[trkBin]>=5000) continue;
if(trkPdg==321 && nPerBinKa[trkBin]>=5000) continue;
+ if(trkPdg==2212 && nPerBinPr[trkBin]>=5000) continue;
if(trkPdg==11 && nPerBinEl[trkBin]>=5000) continue;
+ if(trkPdg==13 && nPerBinMu[trkBin]>=5000) continue;
trkKine[0] = cmptrk.p;
- trkKine[1] = cmptrk.cosl;
// get regularized covariance matrix
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(0,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(0,k);
covmat.c00 = cmptrk.c00/RegFunc(trkKine,trkRegPar);
covmat.c10 = cmptrk.c10;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(1,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(1,k);
covmat.c11 = cmptrk.c11/RegFunc(trkKine,trkRegPar);
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(2,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(2,k);
covmat.c20 = cmptrk.c20/RegFunc(trkKine,trkRegPar);
covmat.c21 = cmptrk.c21;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(3,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(3,k);
covmat.c22 = cmptrk.c22/RegFunc(trkKine,trkRegPar);
covmat.c30 = cmptrk.c30;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(4,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(4,k);
covmat.c31 = cmptrk.c31/RegFunc(trkKine,trkRegPar);
covmat.c32 = cmptrk.c32;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(5,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(5,k);
covmat.c33 = cmptrk.c33/RegFunc(trkKine,trkRegPar);
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(6,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(6,k);
covmat.c40 = cmptrk.c40/RegFunc(trkKine,trkRegPar);
covmat.c41 = cmptrk.c41;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(7,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(7,k);
covmat.c42 = cmptrk.c42/RegFunc(trkKine,trkRegPar);
covmat.c43 = cmptrk.c43;
- for(Int_t k=0;k<4;k++) trkRegPar[k] = (*fRegPar)(8,k);
+ for(Int_t k=0;k<3;k++) trkRegPar[k] = (*fRegPar)(8,k);
covmat.c44 = cmptrk.c44/RegFunc(trkKine,trkRegPar);
// fill the tree
switch (trkPdg) {
case 211: // pions
- CovTreePi_[trkBin].Fill();
+ covTreePi1[trkBin].Fill();
nPerBinPi[trkBin]++;
break;
case 321: // kaons
- CovTreeKa_[trkBin].Fill();
+ covTreeKa1[trkBin].Fill();
nPerBinKa[trkBin]++;
break;
+ case 2212: // protons
+ covTreePr1[trkBin].Fill();
+ nPerBinPr[trkBin]++;
+ break;
case 11: // electrons
- CovTreeEl_[trkBin].Fill();
+ covTreeEl1[trkBin].Fill();
nPerBinEl[trkBin]++;
break;
+ case 13: // muons
+ covTreeMu1[trkBin].Fill();
+ nPerBinMu[trkBin]++;
+ break;
}
} // loop on chain entries
// store all trees the DB file
- TFile *DBfile = new TFile(fDBfileName.Data(),"update");
- DBfile->mkdir("CovMatrices");
+ TFile *dbfile = new TFile(fDBfileName.Data(),"update");
+ dbfile->mkdir("CovMatrices");
gDirectory->cd("/CovMatrices");
gDirectory->mkdir("Pions");
gDirectory->mkdir("Kaons");
+ gDirectory->mkdir("Protons");
gDirectory->mkdir("Electrons");
+ gDirectory->mkdir("Muons");
// store pions
gDirectory->cd("/CovMatrices/Pions");
fDBgridPi.SetName("DBgridPi"); fDBgridPi.Write();
- for(Int_t i=0;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<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<knBinsMu;i++) covTreeMu1[i].Write();
- DBfile->Close();
+ dbfile->Close();
delete [] nPerBinPi;
delete [] nPerBinKa;
+ delete [] nPerBinPr;
delete [] nPerBinEl;
-
+ delete [] nPerBinMu;
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::MakeSeedsFromHits(AliTPC *tpc,TTree *th,
+ TObjArray &seedArray) const {
+//-----------------------------------------------------------------------------
+// This function makes the seeds for tracks from the 1st hits in the TPC
+//-----------------------------------------------------------------------------
+
+ Double_t xg,yg,zg,px,py,pz,pt;
+ Int_t label;
+ Int_t nTracks=(Int_t)th->GetEntries();
+
+ cout<<"+++ Number of \"primary tracks\"(entries in TreeH): "<<nTracks<<
+ "\n";
+
+ AliTPChit *tpcHit=0;
+
+ // loop over entries in TreeH
+ 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);
+ // Get FirstHit
+ tpcHit=(AliTPChit*)tpc->FirstHit(-1);
+ for( ; tpcHit; tpcHit=(AliTPChit*)tpc->NextHit() ) {
+ if(tpcHit->fQ !=0.) continue;
+ // Get particle momentum at hit
+ px=tpcHit->X(); py=tpcHit->Y(); pz=tpcHit->Z();
+
+ pt=TMath::Sqrt(px*px+py*py);
+ // reject hits with Pt<mag*0.45 GeV/c
+ if(pt<(fBz*0.45)) continue;
+
+ // Get track label
+ label=tpcHit->Track();
+
+ if((tpcHit=(AliTPChit*)tpc->NextHit())==0) break;
+ if(tpcHit->fQ != 0.) continue;
+ // Get global coordinates of hit
+ xg=tpcHit->X(); yg=tpcHit->Y(); zg=tpcHit->Z();
+ if(TMath::Sqrt(xg*xg+yg*yg)>90.) continue;
+
+ // create the seed
+ AliTPCseedGeant *ioseed = new AliTPCseedGeant(xg,yg,zg,px,py,pz,label);
+
+ // reject tracks which are not in the TPC acceptance
+ if(!ioseed->InTPCAcceptance()) { delete ioseed; continue; }
+
+ // put seed in array
+ seedArray.AddLast(ioseed);
+ }
+
+ } // loop over entries in TreeH
+
+ return;
+}
+//-----------------------------------------------------------------------------
+void AliTPCtrackerParam::MakeSeedsFromRefs(TTree *ttr,
+ TObjArray &seedArray) const {
+//-----------------------------------------------------------------------------
+// This function makes the seeds for tracks from the track references
+//-----------------------------------------------------------------------------
+
+ Double_t xg,yg,zg,px,py,pz,pt;
+ Int_t label;
+ Int_t nnn,k;
+
+ TClonesArray *tkRefArray = new TClonesArray("AliTrackReference");
+
+ TBranch *b =(TBranch*)ttr->GetBranch("TPC");
+ if(!b) {cerr<<"TPC branch of TreeTR not found"<<endl; return; }
+ b->SetAddress(&tkRefArray);
+ Int_t nTkRef = (Int_t)b->GetEntries();
+ 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(!b->GetEvent(l)) continue;
+ nnn = tkRefArray->GetEntriesFast();
+
+ if(nnn <= 0) continue;
+ for(k=0; k<nnn; k++) {
+ AliTrackReference *tkref = (AliTrackReference*)tkRefArray->UncheckedAt(k);
+
+ xg = tkref->X();
+ yg = tkref->Y();
+ zg = tkref->Z();
+ px = tkref->Px();
+ py = tkref->Py();
+ pz = tkref->Pz();
+ label = tkref->GetTrack();
+
+ pt=TMath::Sqrt(px*px+py*py);
+ // reject hits with Pt<mag*0.45 GeV/c
+ if(pt<(fBz*0.45)) continue;
+
+ // create the seed
+ 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()-83.8) > 0.2) {
+ //cerr<<"not at TPC inner part: XL = "<<ioseed->GetXL()<<endl;
+ delete ioseed; continue;
+ }
+
+ // reject tracks which are not in the TPC acceptance
+ if(!ioseed->InTPCAcceptance()) {
+ delete ioseed; continue;
+ }
+
+ // put seed in array
+ seedArray.AddLast(ioseed);
+
+ }
+
+
+ } // loop on track references
+
+ delete tkRefArray;
+
+
return;
}
//-----------------------------------------------------------------------------
// 2) computes the average TPC efficiencies
//-----------------------------------------------------------------------------
- Char_t *outName="TPCeff.root";
+ const char *outName="TPCeff.root";
- Int_t nCol;
- Float_t effPi,effKa,effPr,effEl,effMu;
- Float_t avEffPi[27],avEffKa[27],avEffPr[27],avEffEl[27],avEffMu[27];
- Int_t evtsPi[27],evtsKa[27],evtsPr[27],evtsEl[27],evtsMu[27];
- for(Int_t j=0; j<27; j++) {
- avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
- evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
- }
+ // merge files with tracks
+ cerr<<" ******** MERGING FILES **********\n\n";
// create the chain for the tree of compared tracks
TChain ch1("cmptrktree");
TChain ch2("cmptrktree");
TChain ch3("cmptrktree");
-
for(Int_t j=evFirst; j<=evLast; j++) {
cerr<<"Processing event "<<j<<endl;
TString covName("CovMatrix.");
- TString effName("TPCeff.");
covName+=j;
- effName+=j;
covName.Append(".root");
- effName.Append(".dat");
- if(gSystem->AccessPathName(covName.Data(),kFileExists) ||
- gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
+ if(gSystem->AccessPathName(covName.Data(),kFileExists)) continue;
+
if(j<=100) ch1.Add(covName.Data());
if(j>100 && j<=200) ch2.Add(covName.Data());
if(j>200) ch3.Add(covName.Data());
-
- FILE *effIn = fopen(effName.Data(),"r");
- for(Int_t k=0; k<27; k++) {
- nCol = fscanf(effIn,"%f %f %f %f %f",&effPi,&effKa,&effPr,&effEl,&effMu);
- if(effPi>=0.) {avEffPi[k]+=effPi; evtsPi[k]++;}
- if(effKa>=0.) {avEffKa[k]+=effKa; evtsKa[k]++;}
- if(effPr>=0.) {avEffPr[k]+=effPr; evtsPr[k]++;}
- if(effEl>=0.) {avEffEl[k]+=effEl; evtsEl[k]++;}
- if(effMu>=0.) {avEffMu[k]+=effMu; evtsMu[k]++;}
- }
- fclose(effIn);
-
} // loop on events
// merge chain in one file
// efficiencies
+ cerr<<" ***** EFFICIENCIES ******\n\n";
+
+ ReadEffs("TPCeff.1.root");
+
+ Int_t n = fEffPi.GetTotPoints();
+ Double_t *avEffPi = new Double_t[n];
+ Double_t *avEffKa = new Double_t[n];
+ Double_t *avEffPr = new Double_t[n];
+ Double_t *avEffEl = new Double_t[n];
+ Double_t *avEffMu = new Double_t[n];
+ Int_t *evtsPi = new Int_t[n];
+ Int_t *evtsKa = new Int_t[n];
+ Int_t *evtsPr = new Int_t[n];
+ Int_t *evtsEl = new Int_t[n];
+ Int_t *evtsMu = new Int_t[n];
+
+ for(Int_t j=0; j<n; j++) {
+ avEffPi[j]=avEffKa[j]=avEffPr[j]=avEffEl[j]=avEffMu[j]=0.;
+ evtsPi[j]=evtsKa[j]=evtsPr[j]=evtsEl[j]=evtsMu[j]=0;
+ }
+
+ for(Int_t j=evFirst; j<=evLast; j++) {
+ cerr<<"Processing event "<<j<<endl;
- InitializeKineGrid("eff","bins");
+ TString effName("TPCeff.");
+ effName+=j;
+ effName.Append(".root");
+
+ if(gSystem->AccessPathName(effName.Data(),kFileExists)) continue;
+
+ ReadEffs(effName.Data());
+
+ for(Int_t k=0; k<n; k++) {
+ if(fEffPi.GetParam(k)>=0.) {avEffPi[k]+=fEffPi.GetParam(k); evtsPi[k]++;}
+ if(fEffKa.GetParam(k)>=0.) {avEffKa[k]+=fEffKa.GetParam(k); evtsKa[k]++;}
+ if(fEffPr.GetParam(k)>=0.) {avEffPr[k]+=fEffPr.GetParam(k); evtsPr[k]++;}
+ if(fEffEl.GetParam(k)>=0.) {avEffEl[k]+=fEffEl.GetParam(k); evtsEl[k]++;}
+ if(fEffMu.GetParam(k)>=0.) {avEffMu[k]+=fEffMu.GetParam(k); evtsMu[k]++;}
+ }
- // pions
- for(Int_t j=0; j<27; j++) {
+ } // loop on events
+
+ // compute average efficiencies
+ for(Int_t j=0; j<n; j++) {
if(evtsPi[j]==0) evtsPi[j]++;
fEffPi.SetParam(j,(Double_t)avEffPi[j]/evtsPi[j]);
- }
-
- // kaons
- Int_t j=0;
- for(Int_t k=0; k<27; k++) {
- if(k>14 && k!=21 && k!=22 && k!=23) continue;
- if(evtsKa[k]==0) evtsKa[k]++;
- fEffKa.SetParam(j,(Double_t)avEffKa[k]/evtsKa[k]);
- j++;
- }
-
- // protons
- for(Int_t j=0; j<15; j++) {
+ if(evtsKa[j]==0) evtsKa[j]++;
+ fEffKa.SetParam(j,(Double_t)avEffKa[j]/evtsKa[j]);
if(evtsPr[j]==0) evtsPr[j]++;
fEffPr.SetParam(j,(Double_t)avEffPr[j]/evtsPr[j]);
- }
-
- // electrons
- for(Int_t j=0; j<27; j++) {
if(evtsEl[j]==0) evtsEl[j]++;
fEffEl.SetParam(j,(Double_t)avEffEl[j]/evtsEl[j]);
- }
-
- // muons
- for(Int_t j=0; j<9; j++) {
if(evtsMu[j]==0) evtsMu[j]++;
fEffMu.SetParam(j,(Double_t)avEffMu[j]/evtsMu[j]);
}
-
+
// write efficiencies to a file
WriteEffs(outName);
+ delete [] avEffPi;
+ delete [] avEffKa;
+ delete [] avEffPr;
+ delete [] avEffEl;
+ delete [] avEffMu;
+ delete [] evtsPi;
+ delete [] evtsKa;
+ delete [] evtsPr;
+ delete [] evtsEl;
+ delete [] evtsMu;
+
return;
}
//-----------------------------------------------------------------------------
if(!ReadPulls(inName)) return 0;
if(!ReadRegParams(inName,211)) return 0;
if(!ReadRegParams(inName,321)) return 0;
+ if(!ReadRegParams(inName,2212)) return 0;
if(!ReadRegParams(inName,11)) return 0;
+ if(!ReadRegParams(inName,13)) return 0;
if(!ReaddEdx(inName,211)) return 0;
if(!ReaddEdx(inName,321)) return 0;
if(!ReaddEdx(inName,2212)) return 0;
if(!ReaddEdx(inName,11)) return 0;
+ if(!ReaddEdx(inName,13)) return 0;
if(!ReadDBgrid(inName)) return 0;
+ if(!ReadCovTrees(inName)) return 0;
+
+ return 1;
+}
+//-----------------------------------------------------------------------------
+Int_t AliTPCtrackerParam::ReadCovTrees(const Char_t* inName) {
+//-----------------------------------------------------------------------------
+// This function reads the covariance matrix trees from the DB
+//-----------------------------------------------------------------------------
+
+ if(gSystem->AccessPathName(inName,kFileExists)) {
+ cerr<<"AliTPCtrackerParam::ReaddEdx: "<<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;
}
fdEdxRMSEl.~AliTPCkineGrid();
new(&fdEdxRMSEl) AliTPCkineGrid(*fdEdxRMS);
break;
+ case 13:
+ if(fColl==0) {
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxMeanPi");
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Pions/dEdxRMSPi");
+ } else {
+ fdEdxMean = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxMeanMu");
+ fdEdxRMS = (AliTPCkineGrid*)inFile->Get("/dEdx/Muons/dEdxRMSMu");
+ }
+ fdEdxMeanMu.~AliTPCkineGrid();
+ new(&fdEdxMeanMu) AliTPCkineGrid(*fdEdxMean);
+ fdEdxRMSMu.~AliTPCkineGrid();
+ new(&fdEdxRMSMu) AliTPCkineGrid(*fdEdxRMS);
+ break;
}
inFile->Close();
fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Kaons/DBgridKa");
fDBgridKa.~AliTPCkineGrid();
new(&fDBgridKa) AliTPCkineGrid(*fDBgrid);
+ if(fColl==0) {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
+ } else {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Protons/DBgridPr");
+ }
+ fDBgridPr.~AliTPCkineGrid();
+ new(&fDBgridPr) AliTPCkineGrid(*fDBgrid);
fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Electrons/DBgridEl");
- fDBgridEl.~AliTPCkineGrid();
+ fDBgridEl.~AliTPCkineGrid();
new(&fDBgridEl) AliTPCkineGrid(*fDBgrid);
+ if(fColl==0) {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Pions/DBgridPi");
+ } else {
+ fDBgrid = (AliTPCkineGrid*)inFile->Get("/CovMatrices/Muons/DBgridMu");
+ }
+ fDBgridMu.~AliTPCkineGrid();
+ new(&fDBgridMu) AliTPCkineGrid(*fDBgrid);
inFile->Close();
for(Int_t i=0; i<5; i++) {
TString pi("/Pulls/Pions/PullsPi_"); pi+=i;
TString ka("/Pulls/Kaons/PullsKa_"); ka+=i;
+ TString pr("/Pulls/Protons/PullsPr_"); pr+=i;
TString el("/Pulls/Electrons/PullsEl_"); el+=i;
+ TString mu("/Pulls/Muons/PullsMu_"); mu+=i;
+
fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
fPullsPi[i].~AliTPCkineGrid();
new(&fPullsPi[i]) AliTPCkineGrid(*fPulls);
+
fPulls = (AliTPCkineGrid*)inFile->Get(ka.Data());
fPullsKa[i].~AliTPCkineGrid();
new(&fPullsKa[i]) AliTPCkineGrid(*fPulls);
+
+ if(fColl==0) {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
+ } else {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pr.Data());
+ }
+ fPullsPr[i].~AliTPCkineGrid();
+ new(&fPullsPr[i]) AliTPCkineGrid(*fPulls);
+
fPulls = (AliTPCkineGrid*)inFile->Get(el.Data());
fPullsEl[i].~AliTPCkineGrid();
new(&fPullsEl[i]) AliTPCkineGrid(*fPulls);
+
+ if(fColl==0) {
+ fPulls = (AliTPCkineGrid*)inFile->Get(pi.Data());
+ } else {
+ fPulls = (AliTPCkineGrid*)inFile->Get(mu.Data());
+ }
+ fPullsMu[i].~AliTPCkineGrid();
+ new(&fPullsMu[i]) AliTPCkineGrid(*fPulls);
}
inFile->Close();
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/Pions/RegPions");
+ } else {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Protons/RegProtons");
+ }
+ new(&fRegParPr) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParPr;
+ fRegParPr(0,0) = dummyMatrix(0,0);
+ fRegParPr(0,1) = dummyMatrix(0,1);
+ fRegParPr(0,2) = dummyMatrix(0,2);
+ fRegParPr(1,0) = dummyMatrix(3,0);
+ fRegParPr(1,1) = dummyMatrix(1,1);
+ fRegParPr(1,2) = dummyMatrix(1,2);
+ fRegParPr(2,0) = dummyMatrix(6,0);
+ fRegParPr(2,1) = dummyMatrix(3,2);
+ fRegParPr(2,2) = dummyMatrix(2,2);
+ fRegParPr(3,0) = dummyMatrix(1,0);
+ fRegParPr(3,1) = dummyMatrix(4,0);
+ fRegParPr(3,2) = dummyMatrix(7,0);
+ fRegParPr(4,0) = dummyMatrix(3,1);
+ fRegParPr(4,1) = dummyMatrix(4,1);
+ fRegParPr(4,2) = dummyMatrix(7,1);
+ fRegParPr(5,0) = dummyMatrix(6,1);
+ fRegParPr(5,1) = dummyMatrix(4,2);
+ fRegParPr(5,2) = dummyMatrix(7,2);
+ fRegParPr(6,0) = dummyMatrix(2,0);
+ fRegParPr(6,1) = dummyMatrix(5,0);
+ fRegParPr(6,2) = dummyMatrix(8,0);
+ fRegParPr(7,0) = dummyMatrix(2,1);
+ fRegParPr(7,1) = dummyMatrix(5,1);
+ fRegParPr(7,2) = dummyMatrix(8,1);
+ fRegParPr(8,0) = dummyMatrix(6,2);
+ fRegParPr(8,1) = dummyMatrix(5,2);
+ fRegParPr(8,2) = dummyMatrix(8,2);
break;
case 11:
fRegPar = (TMatrixD*)inFile->Get("/RegParams/Electrons/RegElectrons");
new(&fRegParEl) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParEl;
+ fRegParEl(0,0) = dummyMatrix(0,0);
+ fRegParEl(0,1) = dummyMatrix(0,1);
+ fRegParEl(0,2) = dummyMatrix(0,2);
+ fRegParEl(1,0) = dummyMatrix(3,0);
+ fRegParEl(1,1) = dummyMatrix(1,1);
+ fRegParEl(1,2) = dummyMatrix(1,2);
+ fRegParEl(2,0) = dummyMatrix(6,0);
+ fRegParEl(2,1) = dummyMatrix(3,2);
+ fRegParEl(2,2) = dummyMatrix(2,2);
+ fRegParEl(3,0) = dummyMatrix(1,0);
+ fRegParEl(3,1) = dummyMatrix(4,0);
+ fRegParEl(3,2) = dummyMatrix(7,0);
+ fRegParEl(4,0) = dummyMatrix(3,1);
+ fRegParEl(4,1) = dummyMatrix(4,1);
+ fRegParEl(4,2) = dummyMatrix(7,1);
+ fRegParEl(5,0) = dummyMatrix(6,1);
+ fRegParEl(5,1) = dummyMatrix(4,2);
+ fRegParEl(5,2) = dummyMatrix(7,2);
+ fRegParEl(6,0) = dummyMatrix(2,0);
+ fRegParEl(6,1) = dummyMatrix(5,0);
+ fRegParEl(6,2) = dummyMatrix(8,0);
+ fRegParEl(7,0) = dummyMatrix(2,1);
+ fRegParEl(7,1) = dummyMatrix(5,1);
+ fRegParEl(7,2) = dummyMatrix(8,1);
+ fRegParEl(8,0) = dummyMatrix(6,2);
+ fRegParEl(8,1) = dummyMatrix(5,2);
+ fRegParEl(8,2) = dummyMatrix(8,2);
+ break;
+ case 13:
+ if(fColl==0) {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Pions/RegPions");
+ } else {
+ fRegPar = (TMatrixD*)inFile->Get("/RegParams/Muons/RegMuons");
+ }
+ new(&fRegParMu) TMatrixD(*fRegPar);
+ dummyMatrix = fRegParMu;
+ fRegParMu(0,0) = dummyMatrix(0,0);
+ fRegParMu(0,1) = dummyMatrix(0,1);
+ fRegParMu(0,2) = dummyMatrix(0,2);
+ fRegParMu(1,0) = dummyMatrix(3,0);
+ fRegParMu(1,1) = dummyMatrix(1,1);
+ fRegParMu(1,2) = dummyMatrix(1,2);
+ fRegParMu(2,0) = dummyMatrix(6,0);
+ fRegParMu(2,1) = dummyMatrix(3,2);
+ fRegParMu(2,2) = dummyMatrix(2,2);
+ fRegParMu(3,0) = dummyMatrix(1,0);
+ fRegParMu(3,1) = dummyMatrix(4,0);
+ fRegParMu(3,2) = dummyMatrix(7,0);
+ fRegParMu(4,0) = dummyMatrix(3,1);
+ fRegParMu(4,1) = dummyMatrix(4,1);
+ fRegParMu(4,2) = dummyMatrix(7,1);
+ fRegParMu(5,0) = dummyMatrix(6,1);
+ fRegParMu(5,1) = dummyMatrix(4,2);
+ fRegParMu(5,2) = dummyMatrix(7,2);
+ fRegParMu(6,0) = dummyMatrix(2,0);
+ fRegParMu(6,1) = dummyMatrix(5,0);
+ fRegParMu(6,2) = dummyMatrix(8,0);
+ fRegParMu(7,0) = dummyMatrix(2,1);
+ fRegParMu(7,1) = dummyMatrix(5,1);
+ fRegParMu(7,2) = dummyMatrix(8,1);
+ fRegParMu(8,0) = dummyMatrix(6,2);
+ fRegParMu(8,1) = dummyMatrix(5,2);
+ fRegParMu(8,2) = dummyMatrix(8,2);
break;
}
inFile->Close();
gStyle->SetOptFit(10001);
Int_t thisPdg=211;
- Char_t *part="Pions - ";
+ const char *part="Pions - ";
- InitializeKineGrid("DB","points");
+ 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
part="Kaons - ";
cerr<<" Processing kaons ...\n";
break;
+ case 2212: //protons
+ thisPdg=2212;
+ part="Protons - ";
+ cerr<<" Processing protons ...\n";
+ break;
case 11: // electrons
thisPdg= 11;
part="Electrons - ";
cerr<<" Processing electrons ...\n";
break;
+ case 13: // muons
+ thisPdg= 13;
+ part="Muons - ";
+ cerr<<" Processing muons ...\n";
+ break;
}
- // create the chain with all compared tracks
- TChain cmptrkchain("cmptrktree");
- cmptrkchain.Add("CovMatrix_AllEvts_1.root");
- cmptrkchain.Add("CovMatrix_AllEvts_2.root");
- cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+
+ /*
+ // create a chain with compared tracks
+ TChain *cmptrkchain = new ("cmptrktree");
+ cmptrkchain.Add("CovMatrix_AllEvts.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_1.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_2.root");
+ //cmptrkchain.Add("CovMatrix_AllEvts_3.root");
+ */
+
+
+ TFile *infile = TFile::Open("CovMatrix_AllEvts.root");
+ TTree *cmptrktree = (TTree*)infile->Get("cmptrktree");
COMPTRACK cmptrk;
- cmptrkchain.SetBranchAddress("comptracks",&cmptrk);
- Int_t entries = (Int_t)cmptrkchain.GetEntries();
+ cmptrktree->SetBranchAddress("comptracks",&cmptrk);
+ Int_t entries = (Int_t)cmptrktree->GetEntries();
Int_t pbin;
- Int_t * n = new Int_t[fitbins];
- Int_t * n00 = new Int_t[fitbins];
- Int_t * n11 = new Int_t[fitbins];
- Int_t * n20 = new Int_t[fitbins];
- Int_t * n22 = new Int_t[fitbins];
- Int_t * n31 = new Int_t[fitbins];
- Int_t * n33 = new Int_t[fitbins];
- Int_t * n40 = new Int_t[fitbins];
- Int_t * n42 = new Int_t[fitbins];
- Int_t * n44 = new Int_t[fitbins];
- Double_t * p = new Double_t[fitbins];
- Double_t * ep = new Double_t[fitbins];
- Double_t * mean00 = new Double_t[fitbins];
- Double_t * mean11 = new Double_t[fitbins];
- Double_t * mean20 = new Double_t[fitbins];
- Double_t * mean22 = new Double_t[fitbins];
- Double_t * mean31 = new Double_t[fitbins];
- Double_t * mean33 = new Double_t[fitbins];
- Double_t * mean40 = new Double_t[fitbins];
- Double_t * mean42 = new Double_t[fitbins];
- Double_t * mean44 = new Double_t[fitbins];
- Double_t * sigma00 = new Double_t[fitbins];
- Double_t * sigma11 = new Double_t[fitbins];
- Double_t * sigma20 = new Double_t[fitbins];
- Double_t * sigma22 = new Double_t[fitbins];
- Double_t * sigma31 = new Double_t[fitbins];
- Double_t * sigma33 = new Double_t[fitbins];
- Double_t * sigma40 = new Double_t[fitbins];
- Double_t * sigma42 = new Double_t[fitbins];
- Double_t * sigma44 = new Double_t[fitbins];
- Double_t * rmean = new Double_t[fitbins];
- Double_t * rsigma = new Double_t[fitbins];
- Double_t fitpar[4];
-
- for(Int_t l=0; l<fitbins; l++) {
+ 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<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.;
// loop on chain entries for mean
for(Int_t l=0; l<entries; l++) {
- cmptrkchain.GetEvent(l);
+ cmptrktree->GetEvent(l);
if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
n[pbin]++;
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];
// loop on chain entries for sigma
for(Int_t l=0; l<entries; l++) {
- cmptrkchain.GetEvent(l);
+ cmptrktree->GetEvent(l);
if(TMath::Abs(cmptrk.pdg)!=thisPdg) continue;
pbin = (Int_t)fDBgrid->GetBin(cmptrk.pt,cmptrk.eta)/fDBgrid->GetBinsEta();
if(TMath::Abs(cmptrk.c00-mean00[pbin])<0.4*mean00[pbin]) { n00[pbin]++;
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]);
// Fit function
- TF1 *func = new TF1("FitRegFunc",FitRegFunc,0.23,50.,4);
- func->SetParNames("A_meas","A_scatt","B1","B2");
+ TF1 *func = new TF1("RegFunc",RegFunc,0.23,50.,3);
+ func->SetParNames("A_meas","A_scatt","B");
// line to draw on the plots
TLine *lin = new TLine(-1,1,1.69,1);
lin->SetLineWidth(2);
// matrix used to store fit results
- TMatrixD fitRes(9,4);
+ TMatrixD fitRes(9,3);
// --- c00 ---
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]");
frame00->Draw();
gr00->Draw("P");
// Sets initial values for parameters
- func->SetParameters(1.6e-3,1.9e-4,1.5,0.);
+ func->SetParameters(1.6e-3,1.9e-4,1.5);
// Fit points in range defined by function
- gr00->Fit("FitRegFunc","R,Q");
+ gr00->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(0,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean00[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma00[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(0,i)=fitpar[i];
+ 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);
- TString regtitle00("C(y,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame00reg->SetXTitle("p [GeV/c]");
// --- c11 ---
-
+
// create the canvas
TCanvas *canv11 = new TCanvas("canv11","c11",0,0,700,900);
canv11->Divide(1,2);
// create the graph for cov matrix
- TGraphErrors *gr11 = new TGraphErrors(fitbins,p,mean11,ep,sigma11);
+ 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]");
frame11->Draw();
gr11->Draw("P");
// Sets initial values for parameters
- func->SetParameters(1.2e-3,8.1e-4,1.,0.);
+ func->SetParameters(1.2e-3,8.1e-4,1.);
// Fit points in range defined by function
- gr11->Fit("FitRegFunc","R,Q");
+ gr11->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(1,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean11[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma11[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(1,i)=fitpar[i];
+ 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);
- TString regtitle11("C(z,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame11reg->SetXTitle("p [GeV/c]");
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]");
frame20->Draw();
gr20->Draw("P");
// Sets initial values for parameters
- func->SetParameters(7.3e-5,1.2e-5,1.5,0.);
+ func->SetParameters(7.3e-5,1.2e-5,1.5);
// Fit points in range defined by function
- gr20->Fit("FitRegFunc","R,Q");
+ gr20->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(2,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean20[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma20[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(2,i)=fitpar[i];
+ 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);
- TString regtitle20("C(#eta, y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame20reg->SetXTitle("p [GeV/c]");
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]");
frame22->Draw();
gr22->Draw("P");
// Sets initial values for parameters
- func->SetParameters(5.2e-6,1.1e-6,2.,1.);
+ func->SetParameters(5.2e-6,1.1e-6,2.);
// Fit points in range defined by function
- gr22->Fit("FitRegFunc","R,Q");
+ gr22->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(3,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean22[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma22[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(3,i)=fitpar[i];
+ 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);
- TString regtitle22("C(#eta, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame22reg->SetXTitle("p [GeV/c]");
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]");
frame31->Draw();
gr31->Draw("P");
// Sets initial values for parameters
- func->SetParameters(-1.2e-5,-1.2e-5,1.5,3.);
+ func->SetParameters(-1.2e-5,-1.2e-5,1.5);
// Fit points in range defined by function
- gr31->Fit("FitRegFunc","R,Q");
+ gr31->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(4,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean31[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = -sigma31[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(4,i)=fitpar[i];
+ 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);
- TString regtitle31("C(tg #lambda,z)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame31reg->SetXTitle("p [GeV/c]");
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]");
frame33->Draw();
gr33->Draw("P");
// Sets initial values for parameters
- func->SetParameters(1.3e-7,4.6e-7,1.7,4.);
+ func->SetParameters(1.3e-7,4.6e-7,1.7);
// Fit points in range defined by function
- gr33->Fit("FitRegFunc","R,Q");
+ gr33->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(5,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean33[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma33[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(5,i)=fitpar[i];
+ 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);
- TString regtitle33("C(tg #lambda,tg #lambda)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame33reg->SetXTitle("p [GeV/c]");
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]");
frame40->Draw();
gr40->Draw("P");
// Sets initial values for parameters
- func->SetParameters(4.e-7,4.4e-8,1.5,0.);
+ func->SetParameters(4.e-7,4.4e-8,1.5);
// Fit points in range defined by function
- gr40->Fit("FitRegFunc","R,Q");
+ gr40->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(6,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean40[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma40[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(6,i)=fitpar[i];
+ 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);
- TString regtitle40("C(C,y)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame40reg->SetXTitle("p [GeV/c]");
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]");
frame42->Draw();
gr42->Draw("P");
// Sets initial values for parameters
- func->SetParameters(3.e-8,8.2e-9,2.,1.);
+ func->SetParameters(3.e-8,8.2e-9,2.);
// Fit points in range defined by function
- gr42->Fit("FitRegFunc","R,Q");
+ gr42->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(7,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean42[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma42[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(7,i)=fitpar[i];
+ 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);
- TString regtitle42("C(C, #eta)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame42reg->SetXTitle("p [GeV/c]");
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]");
frame44->Draw();
gr44->Draw("P");
// Sets initial values for parameters
- func->SetParameters(1.8e-10,5.8e-11,2.,3.);
+ func->SetParameters(1.8e-10,5.8e-11,2.);
// Fit points in range defined by function
- gr44->Fit("FitRegFunc","R,Q");
+ gr44->Fit("RegFunc","R,Q");
func->GetParameters(fitpar);
- for(Int_t i=0; i<4; i++) fitRes(8,i)=fitpar[i];
- for(Int_t l=0; l<fitbins; l++) {
- rmean[l] = mean44[l]/FitRegFunc(&p[l],fitpar);
- rsigma[l] = sigma44[l]/FitRegFunc(&p[l],fitpar);
+ for(Int_t i=0; i<3; i++) fitRes(8,i)=fitpar[i];
+ 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);
- TString regtitle44("C(C,C)/(A_meas+A_scatt/p^{B1}/cos^{B2} #lambda");
+ 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);
frame44reg->SetXTitle("p [GeV/c]");
case 321:
new(&fRegParKa) TMatrixD(fitRes);
break;
+ case 2212:
+ new(&fRegParPr) TMatrixD(fitRes);
+ break;
case 11:
new(&fRegParEl) TMatrixD(fitRes);
break;
+ case 13:
+ new(&fRegParMu) TMatrixD(fitRes);
+ break;
}
// write fit parameters to file
fdEdxRMS = &fdEdxRMSKa;
break;
case 2212:
- fDBgrid = &fDBgridPi;
+ fDBgrid = &fDBgridPr;
fEff = &fEffPr;
- fPulls = fPullsPi;
- fRegPar = &fRegParPi;
+ fPulls = fPullsPr;
+ fRegPar = &fRegParPr;
fdEdxMean = &fdEdxMeanPr;
fdEdxRMS = &fdEdxRMSPr;
break;
fdEdxRMS = &fdEdxRMSEl;
break;
case 13:
- fDBgrid = &fDBgridPi;
+ fDBgrid = &fDBgridMu;
fEff = &fEffMu;
- fPulls = fPullsPi;
- fRegPar = &fRegParPi;
- fdEdxMean = &fdEdxMeanPi;
- fdEdxRMS = &fdEdxRMSPi;
+ fPulls = fPullsMu;
+ fRegPar = &fRegParMu;
+ fdEdxMean = &fdEdxMeanMu;
+ fdEdxRMS = &fdEdxRMSMu;
break;
}
//-----------------------------------------------------------------------------
// 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";
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;
gDirectory->cd("/Pulls");
gDirectory->mkdir("Pions");
gDirectory->mkdir("Kaons");
+ gDirectory->mkdir("Protons");
gDirectory->mkdir("Electrons");
+ gDirectory->mkdir("Muons");
for(Int_t i=0;i<5;i++) {
TString pi("PullsPi_"); pi+=i;
TString ka("PullsKa_"); ka+=i;
+ TString pr("PullsPr_"); pr+=i;
TString el("PullsEl_"); el+=i;
+ TString mu("PullsMu_"); mu+=i;
fPullsPi[i].SetName(pi.Data());
fPullsKa[i].SetName(ka.Data());
+ fPullsPr[i].SetName(pr.Data());
fPullsEl[i].SetName(el.Data());
+ fPullsMu[i].SetName(mu.Data());
gDirectory->cd("/Pulls/Pions");
fPullsPi[i].Write();
gDirectory->cd("/Pulls/Kaons");
fPullsKa[i].Write();
+ gDirectory->cd("/Pulls/Protons");
+ fPullsPr[i].Write();
gDirectory->cd("/Pulls/Electrons");
fPullsEl[i].Write();
+ gDirectory->cd("/Pulls/Muons");
+ fPullsMu[i].Write();
}
outFile->Close();
delete outFile;
//-----------------------------------------------------------------------------
Option_t *opt;
- Char_t *dirName="Pions";
- Char_t *keyName="RegPions";
+ const char *dirName="Pions";
+ const char *keyName="RegPions";
SetParticle(pdg);
dirName="Kaons";
keyName="RegKaons";
break;
+ case 2212:
+ dirName="Protons";
+ keyName="RegProtons";
+ break;
case 11:
dirName="Electrons";
keyName="RegElectrons";
break;
+ case 13:
+ dirName="Muons";
+ keyName="RegMuons";
+ break;
}
TFile *outFile = new TFile(outName,opt);
return 1;
}
+//-----------------------------------------------------------------------------
+//*****************************************************************************
+//-----------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff