#include "AliITSIOTrack.h"
#include <Riostream.h>
ClassImp(AliITSPid)
-
-Float_t AliITSPid::qcorr(Float_t xc)
-{
- assert(0);
- Float_t fcorr;
- fcorr=( 0.766 +0.9692*xc -1.267*xc*xc )*( 1.-TMath::Exp(-xc*64.75) );
- if(fcorr<=0.1)fcorr=0.1;
-return qtot/fcorr;
-}
-
+//
Float_t AliITSPid::qtrm(Int_t track)
{
TVector q(*( this->GetVec(track) ));
}
return qm;
}
-//-----------------------------------------------------------
-//inline
-Int_t AliITSPid::wpik(Int_t nc,Float_t q)
+
+Int_t AliITSPid::wpik(Float_t pm,Float_t q)
{
- // return pion();
+ Double_t par[6];
+ for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
+ fggpi->SetParameters(par);
- Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
- Float_t appi,apk;
- qmpi =cut[nc][1];
- sigpi=cut[nc][2];
- qmk =cut[nc][3];
- sigk =cut[nc][4];
- appi = aprob[0][nc-5];
- apk = aprob[1][nc-5];
-// cout<<"qmpi,sigpi,qmk,sigk="<<qmpi<<" "<<sigpi<<" "<<qmk<<" "<<sigk<<endl;
-// cout<<"appi,apk="<<appi<<","<<apk<<endl;
- Float_t dqpi=(q-qmpi)/sigpi;
- Float_t dqk =(q-qmk )/sigk;
- if( dqk<-1. )return pion();
- dpi =TMath::Abs(dqpi);
- dk =TMath::Abs(dqk);
- //ppi =1.- TMath::Erf(dpi); // +0.5;
- //pk =1.- TMath::Erf(dk); // +0.5;
- ppi=appi*TMath::Gaus(q,qmpi,sigpi)
- /(appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk));
- pk = apk*TMath::Gaus(q,qmk, sigk )
- /(appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk));
+ for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
+ fggka->SetParameters(par);
-// Float_t rpik=ppi/(pk+0.0000001);
-// cout<<"q,dqpi,dqk, wpik: ppi,pk,rpik="
-// <<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
+ Float_t ppi=fggpi->Eval(q);
+ Float_t pka=fggka->Eval(q);
+ Float_t p=ppi+pka;
+ /*
+ if(!fSilent){
+ fggka->Print();
+ fggpi->Print();
+ if(p>0)cout<<" ppi,pka="<<ppi/p<<" "<<pka/p<<endl;
+ }
+ */
- fWp=0.; fWpi=ppi; fWk=pk;
- if( pk>ppi){return kaon();}else{return pion();}
+ if(p>0){
+ ppi=ppi/p;
+ pka=pka/p;
+ fWp=0.; fWpi=ppi; fWk=pka;
+ if( pka>ppi){return fPcode=321;}else{return fPcode=211;}
+ }else{return 0;}
}
//-----------------------------------------------------------
-//inline
-Int_t AliITSPid::wpikp(Int_t nc,Float_t q)
+Int_t AliITSPid::wpikp(Float_t pm,Float_t q)
{
- Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pk,pp;
- Float_t appi,apk,app;
-
- qmpi =cut[nc][1];
- sigpi=cut[nc][2];
- qmk =cut[nc][3];
- sigk =cut[nc][4];
- qmp =cut[nc][5];
- sigp =cut[nc][6];
+ Double_t par[6];
+ for(int i=0;i<6;i++){par[i]=fGGpi[i]->Eval(pm);}
+ fggpi->SetParameters(par);
- //appi = apk = app = 1.;
- appi = aprob[0][nc-5];
- apk = aprob[1][nc-5];
- app = aprob[2][nc-5];
+ for(int i=0;i<3;i++){par[i]=fGGka[i]->Eval(pm);}
+ fggka->SetParameters(par);
- //ppi =TMath::Erf( TMath::Abs((q-qmpi)/sigpi) )+0.5;
- //pka =TMath::Erf( TMath::Abs((q-qmk )/sigk) )+0.5;
- //ppr =TMath::Erf( TMath::Abs((q-qmp )/sigp) )+0.5;
+ for(int i=0;i<3;i++){par[i]=fGGpr[i]->Eval(pm);}
+ fggpr->SetParameters(par);
- ppi=appi*TMath::Gaus(q,qmpi,sigpi)
- /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
- pk = apk*TMath::Gaus(q,qmk, sigk )
- /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
- pp = app*TMath::Gaus(q,qmp, sigp )
- /( appi*TMath::Gaus(q,qmpi,sigpi)+apk*TMath::Gaus(q,qmk,sigk)+app*TMath::Gaus(q,qmp,sigp) );
-
- fWp=pp; fWpi=ppi; fWk=pk;
+ Float_t p,ppi,pka,ppr;
+ if( q>(fggpr->GetParameter(1)+fggpr->GetParameter(2)) )
+ { p=1.0; ppr=1.0; ppi=pka=0.0;
+ }else{
+ ppi=fggpi->Eval(q);
+ pka=fggka->Eval(q);
+ ppr=fggpr->Eval(q);
+ p=ppi+pka+ppr;
+ }
+ if(p>0){
+ ppi=ppi/p;
+ pka=pka/p;
+ ppr=ppr/p;
+ fWp=ppr; fWpi=ppi; fWk=pka;
+ //if(!fSilent)cout<<" ppi,pka,ppr="<<ppi<<" "<<pka<<" "<<ppr<<endl;
-//cout<<" wpikp: mid,sig pi,k,p="<<qmpi<<" "<<sigpi<<"; "<<qmk<<" "<<sigk<<"; "
-// <<qmp<<" "<<sigp<<"; "<<endl;
-// cout<<" aprob: "<<appi<<" "<<apk<<" "<<app<<endl;
-//cout<<" ppi,pk,pp="<<ppi<<" "<<pk<<" "<<pp<<endl;
+ if( ppi>pka&&ppi>ppr )
+ {return fPcode=211;}
+ else{ if(pka>ppr){return fPcode=321;}else{return fPcode=2212;}
+ }
- if( ppi>pk&&ppi>pp ) { return pion(); }
- if(pk>pp){return kaon();}else{return proton();}
+ }else{return 0;}
}
//-----------------------------------------------------------
Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
{
fWpi=fWk=fWp=0.; fPcode=0;
-//1)---------------------- 0-120 MeV/c --------------
- if ( pm<=cut[1][0] )
- { return pion(); }
-//2)----------------------120-200 Mev/c -------------
- if ( pm<=cut[2][0] )
- { if( q<cut[2][2] ){ return pion(); } else { return kaon();} }
-//3)----------------------200-300 Mev/c -------------
- if ( pm<=cut[3][0] )
- if( q<cut[3][2])
- { return pion(); }
- else
- { if ( q<=cut[3][5] ) {return kaon();} else {return proton();}}
-//4)----------------------300-410 Mev/c -------------
- if ( pm<=cut[4][0] )
- if( q<cut[4][2] )
- { return pion(); }
- else
- { if( q<=cut[4][4] ) {return kaon();} else {return proton(); }}
-//5)----------------------410-470 Mev/c -------------
- if ( pm<=cut[5][0] )
- if ( q>cut[5][5] ) {return proton();} else {return wpik(5,q);};
-//6)----------------------470-530 Mev/c -------------
- if ( pm<=cut[6][0] )
- if ( q>cut[6][5] ) {return proton();} else {return wpik(6,q);};
-//7)----------------------530-590 Mev/c -------------
- if ( pm<=cut[7][0] )
- if ( q<=cut[7][5] ) {return wpik(7,q);} else {return proton();};
-//8)----------------------590-650 Mev/c -------------
- if ( pm<=cut[8][0] )
- if ( q<=cut[8][5] ) {return wpik(8,q);} else {return proton();};
-//9)----------------------650-730 Mev/c -------------
- if ( pm<=cut[9][0] )
- if ( q<=cut[9][5] ) {return wpik(9,q);} else {return proton();};
-//10)----------------------730-830 Mev/c -------------
- if( pm<=cut[10][0] ){ return wpikp(10,q); }
-//11)----------------------830-930 Mev/c -------------
- if( pm<=cut[11][0] ){ return wpikp(11,q); }
-//12)----------------------930-1030 Mev/c -------------
- if( pm<=cut[12][0] ){ return wpikp(12,q); }
+ if ( pm<=0.400 )
+ { if( q<fCutKa->Eval(pm) )
+ {return pion();}
+ else{ if( q<fCutPr->Eval(pm) )
+ {return kaon();}
+ else{return proton();}
+ }
+ }
+ if ( pm<=0.750 )
+ if ( q>fCutPr->Eval(pm) )
+ {return proton();} else {return wpik(pm,q);};
+ if( pm<=1.10 ){ return wpikp(pm,q); }
return fPcode;
}
//-----------------------------------------------------------
TClonesArray &arr=*trs;
for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
mxtrs=0;
+ //
+ fCutKa=new TF1("fcutka","pol4",0.05,0.4);
+ Double_t ka[5]={25.616, -161.59, 408.97, -462.17, 192.86};
+ fCutKa->SetParameters(ka);
+ //
+ fCutPr=new TF1("fcutpr","[0]/x/x+[1]",0.05,1.1);
+ Double_t pr[2]={0.70675,0.4455};
+ fCutPr->SetParameters(pr);
+ //
+ //---------- signal fit ----------
+{//Pions
+fGGpi[0]=new TF1("fp1pi","pol4",0.34,1.2);
+ Double_t parpi_0[10]={ -1.9096471071e+03, 4.5354331545e+04, -1.1860738840e+05,
+ 1.1405329025e+05, -3.8289694496e+04 };
+ fGGpi[0]->SetParameters(parpi_0);
+fGGpi[1]=new TF1("fp2pi","[0]/x/x+[1]",0.34,1.2);
+ Double_t parpi_1[10]={ 1.0791668283e-02, 9.7347716496e-01 };
+ fGGpi[1]->SetParameters(parpi_1);
+fGGpi[2]=new TF1("fp3pi","[0]/x/x+[1]",0.34,1.2);
+ Double_t parpi_2[10]={ 5.8191602279e-04, 9.7285601334e-02 };
+ fGGpi[2]->SetParameters(parpi_2);
+fGGpi[3]=new TF1("fp4pi","pol4",0.34,1.2);
+ Double_t parpi_3[10]={ 6.6267353195e+02, 7.1595101104e+02, -5.3095111914e+03,
+ 6.2900977606e+03, -2.2935862292e+03 };
+ fGGpi[3]->SetParameters(parpi_3);
+fGGpi[4]=new TF1("fp5pi","[0]/x/x+[1]",0.34,1.2);
+ Double_t parpi_4[10]={ 9.0419011783e-03, 1.1628922525e+00 };
+ fGGpi[4]->SetParameters(parpi_4);
+fGGpi[5]=new TF1("fp6pi","[0]/x/x+[1]",0.34,1.2);
+ Double_t parpi_5[10]={ 1.8324872519e-03, 2.1503968838e-01 };
+ fGGpi[5]->SetParameters(parpi_5);
+}//End Pions
+{//Kaons
+fGGka[0]=new TF1("fp1ka","pol4",0.24,1.2);
+ Double_t parka_0[20]={
+ -1.1204243395e+02,4.6716191428e+01,2.2584059281e+03,
+ -3.7123338009e+03,1.6003647641e+03 };
+ fGGka[0]->SetParameters(parka_0);
+fGGka[1]=new TF1("fp2ka","[0]/x/x+[1]",0.24,1.2);
+ Double_t parka_1[20]={
+ 2.5181172905e-01,8.7566001814e-01 };
+ fGGka[1]->SetParameters(parka_1);
+fGGka[2]=new TF1("fp3ka","pol6",0.24,1.2);
+ Double_t parka_2[20]={
+ 8.6236021573e+00,-7.0970427531e+01,2.4846827669e+02,
+ -4.6094401290e+02,4.7546751408e+02,-2.5807112462e+02,
+ 5.7545491696e+01 };
+ fGGka[2]->SetParameters(parka_2);
+}//End Kaons
+{//Protons
+fGGpr[0]=new TF1("fp1pr","pol4",0.4,1.2);
+ Double_t parpr_0[10]={
+ 6.0150106543e+01,-8.8176206410e+02,3.1222644604e+03,
+ -3.5269200901e+03,1.2859128345e+03 };
+ fGGpr[0]->SetParameters(parpr_0);
+fGGpr[1]=new TF1("fp2pr","[0]/x/x+[1]",0.4,1.2);
+ Double_t parpr_1[10]={
+ 9.4970837607e-01,7.3573504201e-01 };
+ fGGpr[1]->SetParameters(parpr_1);
+fGGpr[2]=new TF1("fp3pr","[0]/x/x+[1]",0.4,1.2);
+ Double_t parpr_2[10]={
+ 1.2498403757e-01,2.7845072306e-02 };
+ fGGpr[2]->SetParameters(parpr_2);
+}//End Protons
+ //----------- end fit -----------
+
+ fggpr=new TF1("ggpr","gaus",0.4,1.2);
+ fggpi=new TF1("ggpi","gaus+gaus(3)",0.4,1.2);
+ fggka=new TF1("ggka","gaus",0.4,1.2);
+
+ //-------------------------------------------------
const int inf=10;
// Ncut Pmom pilo pihi klo khi plo phi
// cut[j] [0] [1] [2] [3] [4] [5] [6]
aprob[0][7]= 355.; aprob[1][7]=
355.*(20.74/469.5); aprob[2][7]=34.08;
}
-//-----------------------------------------------------------
-
+//End AliITSPid.cxx
#include "../TPC/AliTPCtrack.h"
#include "AliITSIOTrack.h"
#include "AliITStrackV2.h"
+#include <TF1.h>
#include <assert.h>
//___________________________________________________________________________
class AliITSPid :
Int_t fPcode;
//private:
public:
- Float_t qcorr(Float_t);
int qcomp(Float_t* qa,Float_t* qb){return qa[0]>qb[0]?1:0;}
Float_t qtrm(Int_t track);
Float_t qtrm(Float_t qarr[6],Int_t narr);
Float_t fSigmin;
- Int_t wpik(Int_t,Float_t);
- Int_t wpikp(Int_t,Float_t);
+ Int_t wpik(Float_t,Float_t);
+ Int_t wpikp(Float_t,Float_t);
Int_t pion(){return fWpi=1.,fPcode=211;}
Int_t kaon(){return fWk=1.,fPcode=321;}
Int_t proton(){return fWp=1.,fPcode=2212;}
+ Int_t fSilent;
+ TF1* fCutKa;
+ TF1* fCutPr;
+ TF1* fGGpi[6];
+ TF1* fGGka[3];
+ TF1* fGGpr[3];
+ TF1* fggpi;
+ TF1* fggka;
+ TF1* fggpr;
public:
ClassDef(AliITSPid,1) // Class for ITS PID
};
void
AliITSSavePIDV2(Int_t evNumber1=0,Int_t evNumber2=0) {
const char *filename="AliITStracksV2.root";
-
- if (gClassTable->GetID("AliRun") < 0) {
- gROOT->LoadMacro("loadlibs.C"); loadlibs();
- } else { delete gAlice; gAlice=0;
- }
TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
if (!file) file = new TFile(filename);
TFile *fpid = new TFile("AliITStracksV2Pid.root","recreate");
+
+
Float_t factor=0;
- if(gAlice!=0)factor=gAlice->Field()->Factor();
- if(factor==0.)factor=1.;
- AliKalmanTrack::SetConvConst(100/0.299792458/0.2/factor);
+ if(gAlice!=0)factor=gAlice->Field()->SolenoidField();
+ if(factor==0.){
+ cout<<" ================ WARNING ====================================\n";
+ cout<<" The default magnetic field value of 0.4 T will be used\n";
+ cout<<" ==============================================================\n";
+ factor=4.; // Default mag. field = 0.4T
+ }
+ else {
+ cout<<"AliITSSavePIDV2.C: Magnetic field is "<<factor/10.<<" T\n";
+ }
+ AliKalmanTrack::SetConvConst(1000/0.299792458/factor);
+
+
+ AliITSPid *pid=new AliITSPid(100);
//
// Loop over events
//
TTree itspidTree(tpidname,"Tree with PID");
AliITStrackV2Pid *outpid=&pidtmp;
itspidTree.Branch("pids","AliITStrackV2Pid",&outpid,32000,1);
- AliITSPid pid;
-
AliITStrackV2 *iotrack=0;
for (Int_t i=0; i<nentr; i++) {
AliITStrackV2 *iotrack=new AliITStrackV2;
tbranch->SetAddress(&iotrack);
tracktree->GetEvent(i);
- Int_t pcode=pid->GetPcode(iotrack);
+ iotrack->CookdEdx();
Float_t signal=iotrack->GetdEdx();
- //iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
iotrack->PropagateTo(3.,0.0028,65.19);
-
iotrack->PropagateToVertex();
Double_t xk,par[5]; iotrack->GetExternalParameters(xk,par);
Float_t lam=TMath::ATan(par[3]);
Float_t pt_1=TMath::Abs(par[4]);
Float_t mom=0.;
if( (pt_1*TMath::Cos(lam))!=0. ){ mom=1./(pt_1*TMath::Cos(lam)); }else{mom=0.;};
- pidtmp.fPcode=pcode;
+ Float_t phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ signal=signal/35.;
pidtmp.fSignal=signal;
pidtmp.fMom=mom;
+ pidtmp.fPhi=phi;
+ pidtmp.fLam=lam;
+ pidtmp.fGlab=TMath::Abs(iotrack->GetLabel());
+ Int_t pcode=pid->GetPcode(signal,mom);
+ pidtmp.fPcode=pcode;
+ pidtmp.fWpi=pid->GetWpi();
+ pidtmp.fWk=pid->GetWk();
+ pidtmp.fWp=pid->GetWp();
+ //cout<<" pcode,sig,mom="<<pcode<<" "<<signal<<" "<<mom<<endl;
+ //cout<<" wpi,wka,wp="<<pidtmp.fWpi<<" "<<pidtmp.fWk<<" "<<pidtmp.fWp<<endl;
itspidTree.Fill();
delete iotrack;
}// End for i (tracks)
+ cout<<"n ev="<<nev<<endl;
fpid->Write();
}// End for nev (events)
file->Close();
cout<<"File AliITStracksV2Pid.root written"<<endl;
- return;
+ delete file;
+ cout<<"end of AliITSSavePIDV2.C "<<endl;
+ //return;
///////////////////////////////////////////////////////
}
///////////////////////////////////////////////////////////
-// Test macro for AliITStracksV1Pid.root file //
+// Test macro for AliITStracksV2Pid.root file //
// JINR Dubna Jan 2002 //
///////////////////////////////////////////////////////////
void
AliITSScanPIDV2(Int_t evNumber1=0,Int_t evNumber2=0) {
//................. Prepare histogramms ................
- TH2F *qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
- TH2F *qplotP= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
- TH2F *qplotKa= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
- TH2F *qplotPi= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
- TH2F *qplotE= new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
- qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlack); qplotP.SetMarkerSize(.3);
+ TH2F *qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotP= new TH2F("QtrmP","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotKa= new TH2F("QtrmKa","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotPi= new TH2F("QtrmPi","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ TH2F *qplotE= new TH2F("QtrmE","Qtrm vs Pmom",100,0,1.300,100,0,13);
+ qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlue); qplotP.SetMarkerSize(.3);
qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.3);
- qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); qplotPi.SetMarkerSize(.3);
+ qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlack); qplotPi.SetMarkerSize(.3);
qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.3);
//......................................................
TH1F *signal_mip = new TH1F("signal_mip","Signal (mips) for track",100,0.,15.);
TBranch *tbranch=tracktree->GetBranch("pids");
Int_t nentr=tracktree->GetEntries();
- cout<<"Found PID for "<<nentr<<" ITS V1 tracks on "<<tpidname<<endl;
+ cout<<"Found PID for "<<nentr<<" ITS V2 tracks on "<<tpidname<<endl;
AliITStrackV2Pid *iopid=0;
for(Int_t ii=0;ii<nentr;ii++)
signal_mip->Fill(iopid->fSignal);
- if(iopid->fPcode ==2212)qplotP.Fill(1000*iopid->fMom,iopid->fSignal);
- if(iopid->fPcode == 321)qplotKa.Fill(1000*iopid->fMom,iopid->fSignal );
- if(iopid->fPcode == 211)qplotPi.Fill(1000*iopid->fMom,iopid->fSignal );
- if(iopid->fPcode == 11)qplotE.Fill(1000*iopid->fMom,iopid->fSignal );
-
- cout<<"PID pcode,fsignal,fmom= "
- <<iopid->fPcode<<","<<iopid->fSignal<<","<<iopid->fMom<<endl;
+ if(iopid->fPcode ==2212)qplotP.Fill(iopid->fMom,iopid->fSignal);
+ if(iopid->fPcode == 321)qplotKa.Fill(iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 211)qplotPi.Fill(iopid->fMom,iopid->fSignal );
+ if(iopid->fPcode == 11)qplotE.Fill(iopid->fMom,iopid->fSignal );
+ /*
+
+if( (iopid->fWp<0.10)||(iopid->fWk<0.0)||(iopid->fWpi<0.0) ){
+ cout<<"PID pcode,fsignal,fmom= "<<iopid->fPcode<<","<<iopid->fSignal<<","<<iopid->fMom<<endl;
+ cout<<"wpi,wka,wp="<<iopid->fWpi<<" "<<iopid->fWk<<" "<<iopid->fWp<<endl;
+ }
+ */
delete iopid;
}// Enf for ii (tracks)
}// End for nev (events)
c1->cd(2); //gPad->SetFillColor(33);
qplot->Draw();
qplotP.Draw("same"); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same");
+
AliITSPid *pid =new AliITSPid(100);
- c1->Range(0,0,1300,10);
- gStyle->SetLineColor(kRed);
- gStyle->SetLineWidth(2);
- TLine *lj[3],*lk[3];
- for(Int_t j=0;j<3;j++){
- Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
- x1=pid->cut[j+1][0]; x2=pid->cut[j+2][0];
- y1=y2=pid->cut[j+2][2];
- lj[j]=new TLine(1000*x1,y1,1000*x2,y2); lj[j]->Draw();
- if(j==0){yy1=10.;}else{yy1=lj[j-1]->GetY1();}
- yy2=lj[j]->GetY1();
- xx1=xx2=x1;
- lk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); lk[j]->Draw();
- }
- //Draw pions-kaons cuts.
- TLine *mj[7],*mk[7];
- for(Int_t j=0;j<7;j++){
- Float_t x1,x2,y1,y2,xx1,xx2,yy1,yy2;
- x1=pid->cut[j+2][0]; x2=pid->cut[j+3][0];
- y1=y2=pid->cut[j+3][5];
- mj[j]=new TLine(1000*x1,y1,1000*x2,y2); mj[j]->Draw();
- if(j==0){yy1=10.;}else{yy1=mj[j-1]->GetY1();}
- yy2=mj[j]->GetY1();
- xx1=xx2=x1;
- mk[j]=new TLine(1000*xx1,yy1,1000*xx2,yy2); mk[j]->Draw();
- }
- cout<<"End of file AliITStracksV1Pid.root "<<endl;
+ fcutka.Draw("same"); fcutpr.Draw("same");
+ c1->Print("ITSPIDplot.ps");
+
+ cout<<"End of file AliITStracksV2Pid.root "<<endl;
return;
}
#!/bin/sh
-# 9 July 2002,Dubna
-# This script processes the following steps for n events (AliRoot v3-08-03)
+# 12 May 2003,Dubna
+# This script processes the following steps for n events (AliRoot v3-09-08,
+# root v3.05/04)
# (use the parameter n):
# - the TPC and ITS slow simulation (hits+digits+clusters),
# - the TPC+ITS tracking,
# - the TPC and ITS PID
-# (Dubnna version)
-#
+# - the AliITStracksPidV2.root is created with the following information:
+# fGlab - track number
+# fPcode - particle code after the TPC PID
+# fMom - particle momentum (from the AliTPCtracks.root)
+# fLam - particle lambder angle (from the AliTPCtracks.root)
+# fPhi - particle phi angle (from the AliTPCtracks.root)
+# fSignal- TPC signal (mips)
+# fWpi - the PID weight for the identified pion
+# fWk - the PID weight for the identified kaon
+# fWp - the PID weight for the identified proton
+# (the weghts are tuned for the HIJING events)
+# and the AliITSScanPIDV2.C is the example of macros to read the PID
+# information.
+# See ITSPIDplot.ps after run of this script.
+
if [ $# = 0 ]; then nev=1; else nev=$1; fi
if [ $nev = 0 ]; then nev=1; fi
+#
# delete eventual old files from the last run
echo "Start simulation for " $nev " event(s)"
$ALICE_ROOT/ITS/AliITSDeleteOldFiles.sh
#
# run the hit generation
-aliroot -q -b "$ALICE_ROOT/macros/grun.C($nev)"
-#
+aliroot -q -b "$ALICE_ROOT/macros/grun.C($nev)"
# digitize TPC
-aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($nev)"
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($nev)"
+# TPC tracking
+ln -s galice.root digits.root
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCFindClustersMI.C($nev)"
+aliroot -b <<EOI
+.L $ALICE_ROOT/TPC/AliTPCFindTracksMI.C
+AliTPCFindTracks($nev);
+.q
+EOI
#
# digitize ITS
-aliroot -q -b "$ALICE_ROOT/ITS/AliITSMakeSDigits.C($nev)"
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2SDigits.C"
aliroot -q -b "$ALICE_ROOT/ITS/AliITSSDigits2Digits.C"
-# create reconstructed points for the ITS
-aliroot -q -b "$ALICE_ROOT/ITS/AliITSDigits2RecPoints.C"
-
-# do the TPC+ITS tracking
-
-#bad: aliroot -q -b "$ALICE_ROOT/ITS/AliBarrelReconstructionV2.C($1)"
-aliroot -q -b "$ALICE_ROOT/ITS/AliBarrelReconstructionV2.C($nev)"
+# ITS-TPC tracking
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSFindClustersV2.C('s',$nev)"
+aliroot -b <<EOI
+TFile *inkin = TFile::Open("galice.root");
+if (!inkin->IsOpen()) cerr<<"Can't open galice.root !\n";
+if (gAlice) {delete gAlice; gAlice=0;}
+
+gAlice = (AliRun*)inkin->Get("gAlice");
+cout<<"AliRun object found on file "<<gAlice<<endl;
+cout<<"!!!! field ="<<gAlice->Field()->SolenoidField()<<endl;
+AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
+inkin->Close();
+.x $ALICE_ROOT/ITS/AliITSFindTracksV2.C($nev);
+.q
+EOI
#
# Do the PID procedure for the ITS. The output file AliITStrackV2Pid.root is
# created with the information for each track:
-# the ITS ADC trancated mean signal, the particle reconstructed momentum,
-# the probable weights for the different particle kinds (electron, pion, kaon,
-# proton). These weights are under study now.
-#
-aliroot -q -b "$ALICE_ROOT/ITS/AliITSSavePIDV2.C(0,$[$nev-1])"
-#
-# This last line checks the ITS PID only, i.e. creates and draws
-#the dEdx-momentum plot (comment if it's not necessary)
-aliroot "$ALICE_ROOT/ITS/AliITSScanPIDV2.C(0,$[$nev-1])"
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSSavePIDV2.C(0,$[$nev-1])"
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSScanPIDV2.C(0,$[$nev-1])"
+
ClassImp(AliITStrackV2Pid)
+
AliITStrackV2Pid::AliITStrackV2Pid()
{
fWpi=fWk=fWp=0.;
Int_t fGcode,fGlab,fFound;
Float_t fQ1,fQ2,fQ3,fQ4,fQ5,fQ6;
-
+ Float_t fD,fZ;
+
ClassDef(AliITStrackV2Pid,1) // ITS trackV2 PID
};
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff