--- /dev/null
+#include "AliITSPid.h"
+#include "TMath.h"
+#include <iostream.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;
+}
+
+//-----------------------------------------------------------
+#include "vector.h"
+#include <algorithm>
+Float_t AliITSPid::qtrm(Int_t track)
+{
+ TVector q(*( this->GetVec(track) ));
+ Int_t nl=(Int_t)q(0);
+ if(nl<1)return 0.;
+ vector<float> vf(nl);
+ switch(nl){
+ case 1:q(6)=q(1); break;
+ case 2:q(6)=(q(1)+q(2))/2.; break;
+ default:
+ for(Int_t i=0;i<nl;i++){vf[i]=q(i+1);}
+ sort(vf.begin(),vf.end());
+ q(6)= (vf[0]+vf[1])/2.;
+ break;
+ }
+ for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
+ return (q(6));
+}
+//-----------------------------------------------------------
+inline
+Int_t AliITSPid::wpik(Int_t nc,Float_t q)
+{
+ return pion();
+
+ Float_t qmpi,qmk,sigpi,sigk,dpi,dk,ppi,pk;
+ qmpi =cut[nc][1];
+ sigpi=cut[nc][2];
+ qmk =cut[nc][3];
+ sigk =cut[nc][4];
+ Float_t dqpi=(q-qmpi)/sigpi;
+ Float_t dqk =(q-qmk )/sigk;
+ if( dqk<-1. )return pion();
+ dpi =fabs(dqpi);
+ dk =fabs(dqk);
+ ppi =1.- TMath::Erf(dpi); // +0.5;
+ pk =1.- TMath::Erf(dk); // +0.5;
+ Float_t rpik=ppi/(pk+0.0000001);
+ cout<<"q,dqpi,dqk, wpik: ppi,pk,rpik="
+ <<q<<" "<<dqpi<<" "<<dqk<<" "<<ppi<<" "<<pk<<" "<<rpik<<endl;
+ if( rpik>=1000. ){return pion();}
+ if( rpik>0.001 ){fWk=1./(rpik+1.); fWpi=1.-fWk;return 0;}
+ return pion();
+}
+//-----------------------------------------------------------
+inline
+Int_t AliITSPid::wpikp(Int_t nc,Float_t q)
+{
+ return pion();
+
+ Float_t qmpi,qmk,qmp,sigpi,sigk,sigp,ppi,pka,ppr,rppi,rpka;
+
+ qmpi =cut[nc][1];
+ sigpi=cut[nc][2];
+ qmk =cut[nc][3];
+ sigk =cut[nc][4];
+ qmp =cut[nc][5];
+ sigp =cut[nc][6];
+
+ 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;
+
+ rppi=ppr/ppi;
+ rpka=ppr/pka;
+ Float_t rp;
+ if( rppi>rpka ) { rp=rpka; }
+ else{ rp=rppi; }
+ if( rp>=1000. ){ return proton(); }
+ if( rp>0.001 ){ fWp=rp/(rp+1.);return proton(); }
+
+ return 0;
+}
+//-----------------------------------------------------------
+Int_t AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
+{
+ return 0;
+}
+//-----------------------------------------------------------
+Int_t AliITSPid::GetPcode(Float_t q,Float_t pm)
+{
+ fWpi=fWk=fWp=-1.; fPcode=0;
+//1)
+ if ( pm<=cut[1][0] )
+ { return pion(); }
+//2)
+ if ( pm<=cut[2][0] )
+ { if( q<cut[2][2] ){ return pion(); } else { return kaon();} }
+//3)
+ if ( pm<=cut[3][0] )
+ if( q<cut[3][2])
+ { return pion(); }
+ else
+ { if ( q<=cut[3][5] ) {return kaon();} else {return proton();}}
+//4)
+ if ( pm<=cut[4][0] )
+ if( q<cut[4][2] )
+ { return pion(); }
+ else
+ { if( q<=cut[4][4] ) {return kaon();} else {return proton(); }}
+//5)
+ if ( pm<=cut[5][0] )
+ if ( q>cut[5][5] ) {return proton();} else {return wpik(5,q);};
+//6)
+ if ( pm<=cut[6][0] )
+ if ( q>cut[6][5] ) {return proton();} else {return wpik(6,q);};
+//7)
+ if ( pm<=cut[7][0] )
+ if ( q<=cut[7][5] ) {return fPcode;} else {return proton();};
+//8)
+ if ( pm<=cut[8][0] )
+ if ( q<=cut[8][5] ) {return fPcode;} else {return proton();};
+//9)
+ if ( pm<=cut[9][0] )
+ if ( q<=cut[9][5] ) {return fPcode;} else {return proton();};
+//10)
+ if( pm<=cut[10][0] ){ return wpikp(10,q); }
+//11)
+ if( pm<=cut[11][0] ){ return wpikp(11,q); }
+//12)
+ if( pm<=cut[12][0] ){ return wpikp(12,q); }
+
+ return fPcode;
+}
+//-----------------------------------------------------------
+void AliITSPid::SetCut(Int_t n,Float_t pm,Float_t pilo,Float_t pihi,
+ Float_t klo,Float_t khi,Float_t plo,Float_t phi)
+{
+ cut[n][0]=pm;
+ cut[n][1]=pilo;
+ cut[n][2]=pihi;
+ cut[n][3]=klo;
+ cut[n][4]=khi;
+ cut[n][5]=plo;
+ cut[n][6]=phi;
+ return ;
+}
+//-----------------------------------------------------------
+void AliITSPid::SetVec(Int_t ntrack,TVector info)
+{
+TClonesArray& arr=*trs;
+ new( arr[ntrack] ) TVector(info);
+}
+//-----------------------------------------------------------
+TVector* AliITSPid::GetVec(Int_t ntrack)
+{
+TClonesArray& arr=*trs;
+ return (TVector*)arr[ntrack];
+}
+//-----------------------------------------------------------
+void AliITSPid::SetEdep(Int_t track,Float_t Edep)
+{
+ TVector xx(0,11);
+ if( ((TVector*)trs->At(track))->IsValid() )
+ {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
+ Int_t j=(Int_t)xx(0); if(j>4)return;
+ xx(++j)=Edep;xx(0)=j;
+ TClonesArray &arr=*trs;
+ new(arr[track])TVector(xx);
+}
+//-----------------------------------------------------------
+void AliITSPid::SetPmom(Int_t track,Float_t Pmom)
+{
+ TVector xx(0,11);
+ if( ((TVector*)trs->At(track))->IsValid() )
+ {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
+ xx(10)=Pmom;
+ TClonesArray &arr=*trs;
+ new(arr[track])TVector(xx);
+}
+//-----------------------------------------------------------
+void AliITSPid::SetPcod(Int_t track,Int_t partcode)
+{
+ TVector xx(0,11);
+ if( ((TVector*)trs->At(track))->IsValid() )
+ {TVector yy( *((TVector*)trs->At(track)) );xx=yy; }
+ if(xx(11)==0)
+ {xx(11)=partcode; mxtrs++;
+ TClonesArray &arr=*trs;
+ new(arr[track])TVector(xx);
+ }
+}
+//-----------------------------------------------------------
+void AliITSPid::Print(Int_t track)
+{cout<<mxtrs<<" tracks in AliITSPid obj."<<endl;
+ if( ((TVector*)trs->At(track))->IsValid() )
+ {TVector xx( *((TVector*)trs->At(track)) );
+ xx.Print();
+ }
+ else
+ {cout<<"No data for track "<<track<<endl;return;}
+}
+//-----------------------------------------------------------
+void AliITSPid::Tab(void)
+{
+if(trs->GetEntries()==0){cout<<"No entries in TAB"<<endl;return;}
+cout<<"------------------------------------------------------------------------"<<endl;
+cout<<"Nq"<<" q1 "<<" q2 "<<" q3 "<<" q4 "<<" q5 "<<
+ " Qtrm " <<" Wpi "<<" Wk "<<" Wp "<<"Pmom "<<endl;
+cout<<"------------------------------------------------------------------------"<<endl;
+for(Int_t i=0;i<trs->GetEntries();i++)
+{
+ TVector xx( *((TVector*)trs->At(i)) );
+ if( xx.IsValid() && xx(0)>0 )
+ {
+ TVector xx( *((TVector*)trs->At(i)) );
+ if(xx(0)>=2)
+ {
+// 1)Calculate Qtrm
+ xx(6)=(this->qtrm(i));
+
+ }else{
+ xx(6)=xx(1);
+ }
+// 2)Calculate Wpi,Wk,Wp
+ this->GetPcode(xx(6),xx(10)/1000.);
+ xx(7)=GetWpi();
+ xx(8)=GetWk();
+ xx(9)=GetWp();
+// 3)Print table
+ if(xx(0)>0){
+ cout<<xx(0)<<" ";
+ for(Int_t j=1;j<11;j++){
+ cout.width(7);cout.precision(5);cout<<xx(j);
+ }
+ cout<<endl;
+ }
+// 4)Update data in TVector
+ TClonesArray &arr=*trs;
+ new(arr[i])TVector(xx);
+ }
+ else
+ {/*cout<<"No data for track "<<i<<endl;*/}
+}// End loop for tracks
+}
+void AliITSPid::Reset(void)
+{
+ for(Int_t i=0;i<trs->GetEntries();i++){
+ TVector xx(0,11);
+ TClonesArray &arr=*trs;
+ new(arr[i])TVector(xx);
+ }
+}
+//-----------------------------------------------------------
+AliITSPid::AliITSPid(Int_t ntrack)
+{
+ trs = new TClonesArray("TVector",ntrack);
+ TClonesArray &arr=*trs;
+ for(Int_t i=0;i<ntrack;i++)new(arr[i])TVector(0,11);
+ mxtrs=0;
+const int inf=10;
+// Ncut Pmom pilo pihi klo khi plo phi
+// cut[j] [0] [1] [2] [3] [4] [5] [6]
+//----------------------------------------------------------------
+ SetCut( 1, 0.12 , 0. , 0. , inf , inf , inf , inf );
+ SetCut( 2, 0.20 , 0. , 6.0 , 6.0 , inf , inf , inf );
+ SetCut( 3, 0.30 , 0. , 3.5 , 3.5 , 9.0 , 9.0 , inf );
+ SetCut( 4, 0.41 , 0. , 1.9 , 1.9 , 4.0 , 4.0 , inf );
+//----------------------------------------------------------------
+ SetCut( 5, 0.47 , 1. , 0.12, 1.98 , 0.17 , 3.5 , inf );
+ SetCut( 6, 0.53 , 1. , 0.12, 1.75 , 0.16 , 3.0 , inf );
+//----------------------------------------------------------------
+ SetCut( 7, 0.59 , 0. , 0. , 1.18 , 0.125 , 2.7 , inf );
+ SetCut( 8, 0.65 , 0. , 0. , 1.18 , 0.125 , 2.5 , inf );
+ SetCut( 9, 0.73 , 0. , 0. , 0. , 0.125 , 2.0 , inf );
+//----------------------------------------------------------------
+ SetCut( 10, 0.83 , 0. , 0. , 1.25 , 0.13 , 2.14 , 0.20 );
+ SetCut( 11, 0.93 , 0. , 0. , 1.18 , 0.125 , 1.88 , 0.18 );
+ SetCut( 12, 1.03 , 0. , 0. , 1.13 , 0.12 , 1.68 , 0.155);
+}
+//-----------------------------------------------------------
+
--- /dev/null
+#ifndef ALIITSPID_H
+#define ALIITSPID_H
+#include <TObject.h>
+#include <TClonesArray.h>
+#include <TVector.h>
+#include <assert.h>
+//#include <stl.h>
+//___________________________________________________________________________
+class AliITSPid :
+ public TObject {
+
+public:
+ AliITSPid(Int_t ntrs=1000);
+ virtual ~AliITSPid(){}
+ void SetEdep(Int_t track,Float_t Edep);
+ void SetPmom(Int_t track,Float_t Pmom);
+ void SetPcod(Int_t track,Int_t Pcod);
+ void Print(Int_t track);
+ void Tab(void);
+ void Reset(void);
+ void SetVec(Int_t track,TVector info);
+ TVector* GetVec(Int_t track);
+ Int_t GetPcode(TClonesArray*,Float_t);
+ Int_t GetPcode(Float_t,Float_t);
+ void SetCut(Int_t,Float_t,Float_t,Float_t,
+ Float_t,Float_t,Float_t,Float_t);
+ Float_t GetWpi(){return fWpi;}
+ Float_t GetWk(){return fWk;}
+ Float_t GetWp(){return fWp;}
+ Int_t GetPid(){return fPcode;};
+protected:
+public:
+ Float_t cut[13][7];
+ Int_t mxtrs;
+ TClonesArray *trs;
+ Float_t qtot;
+ Float_t fWpi,fWk,fWp;
+ Float_t fRpik,fRppi,fRpka,fRp;
+ 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);
+ Int_t wpik(Int_t,Float_t);
+ Int_t wpikp(Int_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;}
+public:
+ ClassDef(AliITSPid,1) // Class for ITS PID
+};
+
+#endif
+
+
+
+
fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44;
Int_t n=GetNumberOfClusters();
- for (Int_t i=0; i<n; i++) fIndex[i]=t.fIndex[i];
+ //for (Int_t i=0; i<n; i++) fIndex[i]=t.fIndex[i];
+ //b.b.
+ for (Int_t i=0; i<n; i++) {
+ fIndex[i]=t.fIndex[i];
+ fdEdxSample[i]=t.fdEdxSample[i];
+ }
+
}
/*
//_____________________________________________________________________________
//if (TMath::Abs(fP1)>11.5)
//if (fP1*fP4<0) {
- // if (n>kWARN) cout<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
+ // if (n>kWARN) cerr<<"fP1*fP4="<<fP1*fP4<<' '<<fP1<<endl; return 0;}
- if (TMath::Abs(fP2)>=1) {if (n>kWARN) cout<<"fP2="<<fP2<<endl; return 0;}
+ if (TMath::Abs(fP2)>=1) {if (n>kWARN) cerr<<"fP2="<<fP2<<endl; return 0;}
- if (fC00<=0) {if (n>kWARN) cout<<"fC00="<<fC00<<endl; return 0;}
- if (fC11<=0) {if (n>kWARN) cout<<"fC11="<<fC11<<endl; return 0;}
- if (fC22<=0) {if (n>kWARN) cout<<"fC22="<<fC22<<endl; return 0;}
- if (fC33<=0) {if (n>kWARN) cout<<"fC33="<<fC33<<endl; return 0;}
- if (fC44<=0) {if (n>kWARN) cout<<"fC44="<<fC44<<endl; return 0;}
+ if (fC00<=0) {if (n>kWARN) cerr<<"fC00="<<fC00<<endl; return 0;}
+ if (fC11<=0) {if (n>kWARN) cerr<<"fC11="<<fC11<<endl; return 0;}
+ if (fC22<=0) {if (n>kWARN) cerr<<"fC22="<<fC22<<endl; return 0;}
+ if (fC33<=0) {if (n>kWARN) cerr<<"fC33="<<fC33<<endl; return 0;}
+ if (fC44<=0) {if (n>kWARN) cerr<<"fC44="<<fC44<<endl; return 0;}
/*
TMatrixD m(5,5);
m(0,0)=fC00;
Double_t det=m.Determinant();
if (det <= 0) {
- if (n>kWARN) { cout<<" bad determinant "<<det<<endl; m.Print(); }
+ if (n>kWARN) { cerr<<" bad determinant "<<det<<endl; m.Print(); }
return 0;
}
*/
fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
}
+
+void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
+ //-----------------------------------------------------------------
+ // This funtion calculates dE/dX within the "low" and "up" cuts.
+ //-----------------------------------------------------------------
+ Int_t i;
+ Int_t nc=GetNumberOfClusters();
+ if(nc != 6) cout<<"!!!Warning: ncl isn't 6, ="<<nc<<endl;
+
+ // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
+ // Take only SSD and SDD
+ nc=4;
+ Int_t swap;//stupid sorting
+
+ cout<<"Start sorting (low,up,nc)..."<<low<<" "<<up<<" "<<nc<<endl;
+
+ /*
+ for (i=0; i<6; i++) { // b.b.
+ cout<<"! cl befor sort: cl,dEdx ="<<i<<","<<fdEdxSample[i]<<endl;
+ }
+ */
+
+ do {
+ swap=0;
+ for (i=0; i<nc-1; i++) {
+ if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
+ Float_t tmp=fdEdxSample[i];
+ fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
+ swap++;
+ }
+ } while (swap);
+
+ for (i=0; i<nc; i++) { // b.b.
+ cout<<" i, sorted dEdx ="<<i<<","<<fdEdxSample[i]<<endl;
+ }
+ Int_t nl=Int_t(low*nc), nu=Int_t(up*(nc-1)); //b.b. to take two lowest dEdX
+ // values from four ones choose
+ // nu=2
+
+ cout<<" Cook: nl,nu, dEdX samples ="<<nl<<" "<<nu<<" ";
+ Float_t dedx=0;
+ for (i=nl; i<nu; i++){
+ dedx += fdEdxSample[i]; cout<<" "<<fdEdxSample[i]<<" ";}
+ cout<<endl;
+ dedx /= float(nu);
+
+ cout<<"! CookdEdx end: dedx ="<<dedx<<endl;
+ SetdEdx(dedx);
+}
+
+
Int_t Update(const AliCluster* cl,Double_t chi2,UInt_t i);
Int_t Improve(Double_t x0,Double_t yv,Double_t zv);
void SetdEdx(Double_t dedx) {fdEdx=dedx;}
- void CookdEdx(Double_t low=0., Double_t up=1.) {}
+ void SetSampledEdx(Float_t q, Int_t i); // b.b.
+ //void CookdEdx(Double_t low=0., Double_t up=1.) {}
+ void CookdEdx(Double_t low=0., Double_t up=0.8); // b.b.
void SetDetectorIndex(Int_t i) {SetLabel(i);}
void ResetCovariance();
void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
UInt_t fIndex[kMaxLayer]; // indices of associated clusters
+ Float_t fdEdxSample[8]; // array of dE/dx samples b.b.
+
ClassDef(AliITStrackV2,1) //ITS reconstructed track
};
xr=fX;
x[0]=GetY(); x[1]=GetZ(); x[2]=GetSnp(); x[3]=GetTgl(); x[4]=Get1Pt();
}
+//b.b.
+inline
+void AliITStrackV2::SetSampledEdx(Float_t q, Int_t i) {
+ //----------------------------------------------------------------------
+ //
+ //----------------------------------------------------------------------
+ Double_t s=GetSnp(), t=GetTgl();
+ //cout<<"before corr!! i,q="<<i<<" "<<q<<" "<<endl;
+ q *= TMath::Sqrt((1-s*s)/(1+t*t));
+ fdEdxSample[i]=q;
+ //cout<<"after corr!! i,q="<<i<<" "<<q<<" "<<endl;
+}
+
+
#endif
--- /dev/null
+#ifndef ALIITSTRACKV2PID_H
+#define ALIITSTRACKV2PID_H
+
+#include <TMath.h>
+#include <iostream.h>
+#include "AliITStrackV2.h"
+
+//_____________________________________________________________________________
+class AliITStrackV2Pid : public TObject {
+public:
+ AliITStrackV2Pid();
+ virtual ~AliITStrackV2Pid(){}
+public:
+ Float_t fSignal,fMom;
+ Int_t fPcode;
+ Float_t fWpi,fWk,fWp;
+
+ ClassDef(AliITStrackV2Pid,1) // ITS trackV2 PID
+};
+
+#endif
+
+
//cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
return 0;
}
+ // b.b. change for the dEdX
+ fTrackToFollow.SetSampledEdx(c->GetQ(),
+ fTrackToFollow.GetNumberOfClusters()-1);
+
if (constraint) fTrackToFollow.Improve(d/21.82*2.33,GetY(),GetZ());
#ifdef DEBUG
--- /dev/null
+enum gentype_t {hijing, hijingParam, gun, box, pythia,
+ param1, param2, param3, param4,
+ cocktail, fluka, halo, ntuple, scan, doublescan};
+
+//gentype_t gentype=hijingParam;
+//gentype_t gentype=hijing;
+ gentype_t gentype=
+ //box;
+ cocktail;
+ //hijingParam;
+ //hijing;
+
+void Config()
+{
+ // 7-DEC-2000 09:00
+ // Switch on Transition Radiation simulation. 6/12/00 18:00
+ // iZDC=1 7/12/00 09:00
+ // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
+
+ // Set Random Number seed
+ // gRandom->SetSeed(12345);
+
+ new AliGeant3("C++ Interface to Geant3");
+
+ if (!gSystem->Getenv("CONFIG_FILE"))
+ {
+ TFile *rootfile = new TFile("galice.root", "recreate");
+ //TFile *rootfile = new TFile("galice.root", "recreate");
+ rootfile->SetCompressionLevel(2);
+ }
+
+ TGeant3 *geant3 = (TGeant3 *) gMC;
+
+ //
+ // Set External decayer
+ AliDecayer *decayer = new AliDecayerPythia();
+
+ decayer->SetForceDecay(kAll);
+ decayer->Init();
+ gMC->SetExternalDecayer(decayer);
+ //
+ //
+ //=======================================================================
+ // ******* GEANT STEERING parameters FOR ALICE SIMULATION *******
+ geant3->SetTRIG(1); //Number of events to be processed
+ geant3->SetSWIT(4, 10);
+ geant3->SetDEBU(0, 0, 1);
+ //geant3->SetSWIT(2,2);
+ geant3->SetDCAY(1);
+ geant3->SetPAIR(1);
+ geant3->SetCOMP(1);
+ geant3->SetPHOT(1);
+ geant3->SetPFIS(0);
+ geant3->SetDRAY(0);
+ geant3->SetANNI(1);
+ geant3->SetBREM(1);
+ geant3->SetMUNU(1);
+ geant3->SetCKOV(1);
+ geant3->SetHADR(1); //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3)
+ geant3->SetLOSS(2);
+ geant3->SetMULS(1);
+ geant3->SetRAYL(1);
+ geant3->SetAUTO(1); //Select automatic STMIN etc... calc. (AUTO 1) or manual (AUTO 0)
+ geant3->SetABAN(0); //Restore 3.16 behaviour for abandoned tracks
+ geant3->SetOPTI(2); //Select optimisation level for GEANT geometry searches (0,1,2)
+ geant3->SetERAN(5.e-7);
+
+ Float_t cut = 1.e-3; // 1MeV cut by default
+ Float_t tofmax = 1.e10;
+
+ // GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA TOFMAX
+ geant3->SetCUTS(cut, cut, cut, cut, cut, cut, cut, cut, cut, cut,
+ tofmax);
+ //
+ //=======================================================================
+ // ************* STEERING parameters FOR ALICE SIMULATION **************
+ // --- Specify event type to be tracked through the ALICE setup
+ // --- All positions are in cm, angles in degrees, and P and E in GeV
+ if (gSystem->Getenv("CONFIG_NPARTICLES"))
+ {
+ int nParticles = atoi(gSystem->Getenv("CONFIG_NPARTICLES"));
+ } else
+ {
+ int nParticles = 50;
+ }
+//-----------------------------------------------------------
+switch(gentype)
+{
+ case hijing:
+//-----------------------------------------HIJING------------
+AliGenHijing *gener = new AliGenHijing(-1);
+ gener->SetEnergyCMS(5500);
+ gener->SetReferenceFrame("CMS ");
+ gener->SetProjectile("A ", 208, 82);
+ gener->SetTarget ("A ", 208, 82);
+ gener->SetImpactParameterRange(0, 3.);
+ gener->SetEvaluate(1);
+ gener->KeepFullEvent();
+ gener->SetJetQuenching(1);
+ gener->SetShadowing(1);
+ gener->SetDecaysOff(1);
+ gener->SetTrigger(0);
+ gener->SetSelectAll(0);
+ gener->SetMomentumRange( 0.05 , 1.3 );
+ gener->SetPhiRange( 0.0 , 1.0 );
+// gener->SetPhiRange(-180.,180.)
+ gener->SetThetaRange(45.0,135.);
+ gener->SetFlavor(0); // 0 - all, 4 - charm and beaty
+ gener->SetOrigin(0., 0.0 ,0);
+ gener->SetSigma(0,0,5.3);
+ gener->SetVertexSmear(kPerEvent);
+// gener->SetTrackingFlag(0);
+ gener->Init();
+ break;
+//-----------------------------------------HijingPARA----------
+ case hijingParam:
+
+ AliGenHIJINGpara *gener = new AliGenHIJINGpara(20);
+
+ gener->SetMomentumRange(0.1, .7);
+ gener->SetThetaRange(45.,135.);
+ gener->SetPhiRange(0, 360);
+ // gener->SetThetaRange(0.28,179.72);
+ // gener->SetThetaRange(45., 135.);
+ gener->SetOrigin(0, 0, 0); //vertex position
+ gener->SetSigma(0, 0, 0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener->Init();
+ break;
+//------------------------------------------ Box --------------
+ case box:
+ AliGenBox *gener = new AliGenBox(4000); //ntracks=100
+ gener->SetMomentumRange(0.2,1.0);
+ gener->SetPhiRange(0,360);
+ gener->SetThetaRange(45, 135. );
+ gener->SetOrigin(0,0,0);
+ //vertex position
+ gener->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener->SetPart(321); // K+
+ //gener->SetPart(211); // Pi+
+ break;
+
+ //------------------------------------------------------------
+ case cocktail:
+ AliGenCocktail *gener = new AliGenCocktail();
+ int Nmax=50;
+ AliGenBox *gener1 = new AliGenBox(2*Nmax); //ntracks=100
+ gener1->SetMomentumRange(.050,1.00);
+ //gener1->SetMomentumRange(.050,0.80);
+ gener1->SetPhiRange(0,360);
+ gener1->SetThetaRange(45, 135. );
+ gener1->SetOrigin(0,0,0);
+ gener1->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener1->SetPart(321); // K+
+
+ AliGenBox *gener2 = new AliGenBox(2*Nmax); //ntracks=100
+ gener2->SetMomentumRange(.050,1.00);
+ //gener1->SetMomentumRange(.050,0.80);
+ gener2->SetPhiRange(0,360);
+ gener2->SetThetaRange(45, 135. );
+ gener2->SetOrigin(0,0,0);
+ gener2->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener2->SetPart(2212); // Protons
+
+ AliGenBox *gener3 = new AliGenBox(Nmax); //ntracks=100
+ gener3->SetMomentumRange(.050,1.00);
+ //gener1->SetMomentumRange(.050,0.80);
+ gener3->SetPhiRange(0,360);
+ gener3->SetThetaRange(45, 135. );
+ gener3->SetOrigin(0,0,0);
+ gener3->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener3->SetPart(11); // e+
+
+
+ AliGenBox *gener4 = new AliGenBox(Nmax); //ntracks=100
+ gener4->SetMomentumRange(.050,1.00);
+ //gener1->SetMomentumRange(.050,0.80);
+ gener4->SetPhiRange(0,360);
+ gener4->SetThetaRange(45, 135. );
+ gener4->SetOrigin(0,0,0);
+ gener4->SetSigma(0,0,0); //Sigma in (X,Y,Z) (cm) on IP position
+ gener4->SetPart(211); // pi+
+
+ gener->SetMomentumRange(.050,1.00);
+ gener->SetPhiRange(0,360);
+ gener->SetThetaRange(45., 135. );
+ gener->SetOrigin(0,0,0);
+ gener->SetSigma(0,0,0);
+ gener->AddGenerator(gener1,"box1",0.5);
+ gener->AddGenerator(gener2,"box2",0.5);
+ gener->AddGenerator(gener3,"box3",0.5);
+ gener->AddGenerator(gener4,"box4",0.5);
+ gener->Init();
+ break;
+
+ default:
+ break;
+}
+ //
+ // Activate this line if you want the vertex smearing to happen
+ // track by track
+ //
+ //gener->SetVertexSmear(perTrack);
+
+ gAlice->SetField(-999, 2); //Specify maximum magnetic field in Tesla (neg. ==> default field)
+
+ Int_t iABSO = 0;
+ Int_t iCASTOR = 0;
+ Int_t iDIPO = 0;
+ Int_t iFMD = 0;
+ Int_t iFRAME = 0;
+ Int_t iHALL = 0;
+ Int_t iITS = 1;
+ Int_t iMAG = 0;
+ Int_t iMUON = 0;
+ Int_t iPHOS = 0;
+ Int_t iPIPE = 1;
+ Int_t iPMD = 0;
+ Int_t iRICH = 0;
+ Int_t iSHIL = 0;
+ Int_t iSTART = 0;
+ Int_t iTOF = 0;
+ Int_t iTPC = 1;
+ Int_t iTRD = 0;
+ Int_t iZDC = 0;
+
+ //=================== Alice BODY parameters =============================
+ AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
+
+
+ if (iMAG)
+ {
+ //=================== MAG parameters ============================
+ // --- Start with Magnet since detector layouts may be depending ---
+ // --- on the selected Magnet dimensions ---
+ AliMAG *MAG = new AliMAG("MAG", "Magnet");
+ }
+
+
+ if (iABSO)
+ {
+ //=================== ABSO parameters ============================
+ AliABSO *ABSO = new AliABSOv0("ABSO", "Muon Absorber");
+ }
+
+ if (iDIPO)
+ {
+ //=================== DIPO parameters ============================
+
+ AliDIPO *DIPO = new AliDIPOv2("DIPO", "Dipole version 2");
+ }
+
+ if (iHALL)
+ {
+ //=================== HALL parameters ============================
+
+ AliHALL *HALL = new AliHALL("HALL", "Alice Hall");
+ }
+
+
+ if (iFRAME)
+ {
+ //=================== FRAME parameters ============================
+
+ AliFRAME *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
+
+ }
+
+ if (iSHIL)
+ {
+ //=================== SHIL parameters ============================
+
+ AliSHIL *SHIL = new AliSHILv0("SHIL", "Shielding");
+ }
+
+
+ if (iPIPE)
+ {
+ //=================== PIPE parameters ============================
+
+ AliPIPE *PIPE = new AliPIPEv0("PIPE", "Beam Pipe");
+ }
+
+ if(iITS) {
+
+//=================== ITS parameters ============================
+ //
+ // As the innermost detector in ALICE, the Inner Tracking System "impacts" on
+ // almost all other detectors. This involves the fact that the ITS geometry
+ // still has several options to be followed in parallel in order to determine
+ // the best set-up which minimizes the induced background. All the geometries
+ // available to date are described in the following. Read carefully the comments
+ // and use the default version (the only one uncommented) unless you are making
+ // comparisons and you know what you are doing. In this case just uncomment the
+ // ITS geometry you want to use and run Aliroot.
+ //
+ // Detailed geometries:
+ //
+ //
+ //AliITS *ITS = new AliITSv5symm("ITS","Updated ITS TDR detailed version with symmetric services");
+ //
+ //AliITS *ITS = new AliITSv5asymm("ITS","Updates ITS TDR detailed version with asymmetric services");
+ //
+ AliITSvPPRasymm *ITS = new AliITSvPPRasymm("ITS","New ITS PPR detailed version with asymmetric services");
+ ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRasymm2.det"); // don't touch this parameter if you're not an ITS developer
+ ITS->SetThicknessDet1(300.); // detector thickness on layer 1 must be in the range [100,300]
+ ITS->SetThicknessDet2(300.); // detector thickness on layer 2 must be in the range [100,300]
+ ITS->SetThicknessChip1(300.); // chip thickness on layer 1 must be in the range [150,300]
+ ITS->SetThicknessChip2(300.); // chip thickness on layer 2 must be in the range [150,300]
+ ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ //
+ //AliITSvPPRsymm *ITS = new AliITSvPPRsymm("ITS","New ITS PPR detailed version with symmetric services");
+ //ITS->SetMinorVersion(2); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetReadDet(kFALSE); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetWriteDet("$ALICE_ROOT/ITS/ITSgeometry_vPPRsymm2.det"); // don't touch this parameter if you're not an ITS developer
+ //ITS->SetThicknessDet1(300.); // detector thickness on layer 1 must be in the range [100,300]
+ //ITS->SetThicknessDet2(300.); // detector thickness on layer 2 must be in the range [100,300]
+ //ITS->SetThicknessChip1(300.); // chip thickness on layer 1 must be in the range [150,300]
+ //ITS->SetThicknessChip2(300.); // chip thickness on layer 2 must be in the range [150,300]
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetCoolingFluid(1); // 1 --> water ; 0 --> freon
+ //
+ //
+ // Coarse geometries (warning: no hits are produced with these coarse geometries and they unuseful
+ // for reconstruction !):
+ //
+ //
+ //AliITSvPPRcoarseasymm *ITS = new AliITSvPPRcoarseasymm("ITS","New ITS PPR coarse version with asymmetric services");
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+ //
+ //AliITS *ITS = new AliITSvPPRcoarsesymm("ITS","New ITS PPR coarse version with symmetric services");
+ //ITS->SetRails(1); // 1 --> rails in ; 0 --> rails out
+ //ITS->SetSupportMaterial(0); // 0 --> Copper ; 1 --> Aluminum ; 2 --> Carbon
+ //
+ //
+ //
+ // Geant3 <-> EUCLID conversion
+ // ============================
+ //
+ // SetEUCLID is a flag to output (=1) or not to output (=0) both geometry and
+ // media to two ASCII files (called by default ITSgeometry.euc and
+ // ITSgeometry.tme) in a format understandable to the CAD system EUCLID.
+ // The default (=0) means that you dont want to use this facility.
+ //
+ ITS->SetEUCLID(0);
+ }
+
+
+ if (iTPC)
+ {
+ //============================ TPC parameters ================================
+ // --- This allows the user to specify sectors for the SLOW (TPC geometry 2)
+ // --- Simulator. SecAL (SecAU) <0 means that ALL lower (upper)
+ // --- sectors are specified, any value other than that requires at least one
+ // --- sector (lower or upper)to be specified!
+ // --- Reminder: sectors 1-24 are lower sectors (1-12 -> z>0, 13-24 -> z<0)
+ // --- sectors 25-72 are the upper ones (25-48 -> z>0, 49-72 -> z<0)
+ // --- SecLows - number of lower sectors specified (up to 6)
+ // --- SecUps - number of upper sectors specified (up to 12)
+ // --- Sens - sensitive strips for the Slow Simulator !!!
+ // --- This does NOT work if all S or L-sectors are specified, i.e.
+ // --- if SecAL or SecAU < 0
+ //
+ //
+ //-----------------------------------------------------------------------------
+
+ // gROOT->LoadMacro("SetTPCParam.C");
+ // AliTPCParam *param = SetTPCParam();
+ AliTPC *TPC = new AliTPCv2("TPC", "Default");
+
+ // All sectors included
+ TPC->SetSecAL(-1);
+ TPC->SetSecAU(-1);
+
+ }
+
+ if (iTOF)
+ {
+ //=================== TOF parameters ============================
+ AliTOF *TOF = new AliTOFv2("TOF", "normal TOF");
+ }
+
+ if (iRICH)
+ {
+ //=================== RICH parameters ===========================
+ AliRICH *RICH = new AliRICHv1("RICH", "normal RICH");
+
+ }
+
+
+ if (iZDC)
+ {
+ //=================== ZDC parameters ============================
+
+ AliZDC *ZDC = new AliZDCv1("ZDC", "normal ZDC");
+ }
+
+ if (iCASTOR)
+ {
+ //=================== CASTOR parameters ============================
+
+ AliCASTOR *CASTOR = new AliCASTORv1("CASTOR", "normal CASTOR");
+ }
+
+ if (iTRD)
+ {
+ //=================== TRD parameters ============================
+
+ AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
+
+ // Select the gas mixture (0: 97% Xe + 3% isobutane, 1: 90% Xe + 10% CO2)
+ TRD->SetGasMix(1);
+
+ // With hole in front of PHOS
+ TRD->SetPHOShole();
+ // With hole in front of RICH
+ TRD->SetRICHhole();
+ // Switch on TR
+ AliTRDsim *TRDsim = TRD->CreateTR();
+ }
+
+ if (iFMD)
+ {
+ //=================== FMD parameters ============================
+
+ AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
+ }
+
+ if (iMUON)
+ {
+ //=================== MUON parameters ===========================
+
+ AliMUON *MUON = new AliMUONv1("MUON", "default");
+ }
+ //=================== PHOS parameters ===========================
+
+ if (iPHOS)
+ {
+ AliPHOS *PHOS = new AliPHOSv1("PHOS", "GPS2");
+ }
+
+
+ if (iPMD)
+ {
+ //=================== PMD parameters ============================
+
+ AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
+
+ PMD->SetPAR(1., 1., 0.8, 0.02);
+ PMD->SetIN(6., 18., -580., 27., 27.);
+ PMD->SetGEO(0.0, 0.2, 4.);
+ PMD->SetPadSize(0.8, 1.0, 1.0, 1.5);
+
+ }
+
+ if (iSTART)
+ {
+ //=================== START parameters ============================
+ AliSTART *START = new AliSTARTv1("START", "START Detector");
+ }
+
+
+}
--- /dev/null
+{
+ AliKalmanTrack::SetConvConst(100/.299792458/.2);
+ cout<<" Starting..."<<endl;
+if(gClassTable->GetID("AliRun")<0){
+ gROOT->LoadMacro("$ALICE_ROOT/macros/loadlibs.C");
+ loadlibs();
+}
+
+if(gROOT->IsBatch()){
+ cout<<"Start tracking..."<<endl;
+ gSystem->Exec("make all");
+ cout<<"AliITSv2PID.root written."<<endl;
+}else{
+
+ fpid = new TFile("pidhit.root","RECREATE");
+
+ gROOT->Macro("$ALICE_ROOT/ITS/load_particles.C");
+ AliITSPid *pid = new AliITSPid(npart);
+ gROOT->LoadMacro("$ALICE_ROOT/ITS/dEdXgeant.C");
+ gROOT->LoadMacro("$ALICE_ROOT/ITS/dedxanal.C");
+//----------------------------------------------------
+NStat=pid.trs->GetEntries();
+//----------------------------------------------------
+ TControlBar menu("vertical","PID menu",920,5);
+
+ menu.AddButton("dEdX.C","pid->Reset();totpid=0;dEdXyy(0,0,pid);pid->Tab();"," Create new PID table ");
+ menu.AddButton("Save TAB",
+ "pid->Tab();fpid->cd();pid->Write();fpid->Close();"," ");
+ menu.AddButton("Load TAB","loadpid();"," ");
+
+ menu.AddButton("EFFALL","effall(); "," Efficiency if PID ");
+ menu.AddButton("dEdX spectra","qhisall(); "," dEdX for PI,K and P ");
+ menu.AddButton("dEdX-P plot","dedxhis(0); "," dEdX-P plot for PI,K,P ");
+ menu.AddButton("dEdX-P pions","dedxhis(211); "," dEdX-P plot for PI ");
+ menu.AddButton("dEdX-P kaons","dedxhis(321); "," dEdX-P plot for K ");
+ menu.AddButton("dEdX-P elect","dedxhis(11); "," dEdX-P plot for e+ ");
+ menu.AddButton("dEdX-P prot ","dedxhis(2212); "," dEdX-P plot for P ");
+
+ menu.AddButton("Fit Kaons","fitkall(); "," Gaus Fit for Kaons ");
+ menu.AddButton("Fit Pions","fitpiall(); "," Gaus Fit for Kaons ");
+ menu.AddButton("Fit Protons","fitpall(); "," Gaus Fit for Protons ");
+ menu.AddButton("New cuts","newcuts(); "," Corrected cuts for PID object ");
+ menu.AddButton("pcode","pcode(); "," ... ");
+ menu.AddButton("signal (mip)","signal(); "," ... ");
+ menu.AddButton("pmom (MeV)","pmom(); "," ... ");
+ menu.AddButton("tracks","tracks(); "," Track number histogram ");
+ menu.AddButton("1 track","track(); "," Print next track ");
+ menu.AddButton("test module","dEdXxx(0,0,pid,1); "," ... ");
+ menu.AddButton("fill tab test","filltab(); ","Fill track table with test data ");
+ menu.AddButton("fill tab_tr","filltab_tracks(); ","Fill track table with reconstr. tracks ");
+
+ menu.AddButton("Config.C","gSystem->Exec(\"make conf\");","Edit Config.C");
+ menu.AddButton("Do tracking",
+ "gSystem->Exec(\"make all\");pid->Reset();totpid=0;filltab_tracks();dedxhis(0)",
+ "Start tracking");
+ menu.AddButton("Exit","quit();","Quit");
+ menu.Show();
+
+}//if batch
+
+}
+
+
+
+
+
+
+
+
--- /dev/null
+#!/bin/sh
+# delete eventual old files from the last run
+./DeleteOldV2Files
+# run the hit generation
+aliroot -q -b "$ALICE_ROOT/macros/grun.C"
+# digitize TPC and prepare TPC tracks for matching with the ITS
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCtest.C"
+# digitize ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHitsToDigitsDefault.C"
+# create reconstructed points for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/ITSDigitsToClusters.C"
+#prepare slow or fast recpoints for the tracking
+aliroot -q -b "$ALICE_ROOT/ITS/"AliITSFindClustersV2.C\(\'$2\'\)
+# ITS tracking
+aliroot -q -b "$ALICE_ROOT/ITS/SetConvConst.C" "$ALICE_ROOT/ITS/AliITSFindTracksV2.C"
+# do the PID
+aliroot -q -b "$ALICE_ROOT/ITS/SetConvConst.C" "$ALICE_ROOT/ITS/save_pid.C"
+aliroot "$ALICE_ROOT/ITS/ITSPIDtest.C"
+#
+# After all of the above the PID menu will appear. To check the PID package
+# put the buttons:
+#
+#1. "fill tab_tr" and "dEdX-P plot" (or dEdX-P pions, ... kaons, ... elect,
+# ... prot) and the dEdX versus momentum taken from the tracking for all
+# particle (pions, kaons, electrons, protons) will be drown.
+#
+#2. "dEdX.C" and the same other ones as for the point 1. In this case the
+# same plots but fot the generated particle momentum will be drown.
+#
+# The output file AliITStrackV2Pid.root will be written 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.
+
+
+
+
+
+
+
+
AliITSRad.cxx AliITSgeoinfo.cxx AliITSTrackerV1.cxx\
AliITSvtest.cxx \
AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \
+ AliITSPid.cxx AliITStrackV2Pid.cxx \
AliITSVertex.cxx \
- AliV0vertex.cxx AliV0vertexer.cxx \
- AliITSDigitizer.cxx
+ AliV0vertex.cxx AliV0vertexer.cxx
# AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \
# Fortran sources
--- /dev/null
+{
+AliKalmanTrack::SetConvConst(100/.299792458/.2);
+}
--- /dev/null
+
+Int_t max_modules=2269;
+
+void dEdXyy (Int_t evNumber1=0,Int_t evNumber2=0,AliITSPid* pid=0,int mtest=0)
+{
+ if (gClassTable->GetID("AliRun") < 0) {
+ gROOT->LoadMacro("loadlibs.C");
+ loadlibs();
+ } else {
+ delete gAlice;
+ gAlice=0;
+ }
+
+// Connect the Root Galice file containing Geometry, Kine and Hits
+
+ TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+ printf("file %p\n",file);
+ if (file) file->Close();
+ file = new TFile("galice.root","UPDATE");
+ file->ls();
+
+ printf ("I'm after Map \n");
+
+// Get AliRun object from file or create it if not on file
+
+ if (!gAlice) {
+ gAlice = (AliRun*)file->Get("gAlice");
+ if (gAlice) printf("AliRun object found on file\n");
+ if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
+ }
+ printf ("I'm after gAlice \n");
+
+ // Create Histogramms -------------------------------------------------
+
+ TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,400.);
+ //TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,400.);
+ TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,3000.);
+ //TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,6.);
+ TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,30.);
+ //TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,2000.);
+ TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,15000.);
+ //TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,6.);
+ TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,30.);
+ TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,400.);
+ //TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,400.);
+ TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,3000.);
+ //TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,6.);
+ TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,30.);
+ //TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,200.);
+ TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.);
+ //TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,6.);
+ TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,30.);
+
+
+
+ /*
+ TH1F *dedxSDD = new TH1F("dedxSDD","Particle energy (KeV) for layer 3 and 4",100,0.,4000.);
+ TH1F *dedxSDDnorm = new TH1F("dedxSDDnorm","Particle energy (KeV) for layer 3 and 4",100,0.,4000.);
+ TH1F *dedxSDDmip = new TH1F("dedxSDDmip","Particle energy (mips) for layer 3 and 4",100,0.,40.);
+ TH1F *signalSDD = new TH1F("signalSDD","SDD signal (ADC) for layer 3 and 4",100,0.,20000.);
+ TH1F *signalSDDmip = new TH1F("signalSDDmip","SDD signal (mips) for layer 3 and 4",100,0.,40.);
+ TH1F *dedxSSD = new TH1F("dedxSSD","Particle energy (KeV) for layer 5 and 6",100,0.,4000.);
+ TH1F *dedxSSDnorm = new TH1F("dedxSSDnorm","Particle energy (KeV) for layer 5 and 6",100,0.,4000.);
+ TH1F *dedxSSDmip = new TH1F("dedxSSDmip","Particle energy (mips) for layer 5 and 6",100,0.,40.);
+ TH1F *signalSSD = new TH1F("signalSSD","SSD signal (ADC) for layer 5 and 6",100,0.,2000.);
+ TH1F *signalSSDmip = new TH1F("signalSSDmip","SSD signal (mips) for layer 5 and 6",100,0.,40.);
+ */
+
+ TH1F *pathSDD = new TH1F("pathSDD","Path length in the SDD",100,0.,1000.);
+ TH1F *pathSSD = new TH1F("pathSSD","Path length in the SSD",100,0.,1000.);
+
+ // -------------------------------------------------------------
+ AliITS *ITS = (AliITS*) gAlice->GetModule("ITS");
+ if (!ITS) return;
+ AliITSgeom *geom = ITS->GetITSgeom();
+
+ // Set the models for cluster finding
+
+ // SPD
+
+ ITS->MakeTreeC();
+ Int_t nparticles=gAlice->GetEvent(0);
+
+ AliITSDetType *iDetType=ITS->DetType(0);
+ AliITSsegmentationSPD *seg0=(AliITSsegmentationSPD*)iDetType->GetSegmentationModel();
+ TClonesArray *dig0 = ITS->DigitsAddress(0);
+ TClonesArray *recp0 = ITS->ClustersAddress(0);
+ AliITSClusterFinderSPD *rec0=new AliITSClusterFinderSPD(seg0,dig0,recp0);
+ ITS->SetReconstructionModel(0,rec0);
+
+ // SDD
+
+ AliITSDetType *iDetType=ITS->DetType(1);
+ AliITSgeom *geom = ITS->GetITSgeom();
+
+ AliITSsegmentationSDD *seg1=(AliITSsegmentationSDD*)iDetType->GetSegmentationModel();
+ if (!seg1) seg1 = new AliITSsegmentationSDD(geom);
+ AliITSresponseSDD *res1 = (AliITSresponseSDD*)iDetType->GetResponseModel();
+ if (!res1) res1=new AliITSresponseSDD();
+
+ Float_t baseline,noise;
+ res1->GetNoiseParam(noise,baseline);
+ Float_t noise_after_el = res1->GetNoiseAfterElectronics();
+ Float_t thres = baseline;
+ thres += (4.*noise_after_el); // TB // (4.*noise_after_el);
+ printf("thres %f\n",thres);
+ //res1->Print();
+
+ TClonesArray *dig1 = ITS->DigitsAddress(1);
+ TClonesArray *recp1 = ITS->ClustersAddress(1);
+ AliITSClusterFinderSDD *rec1=new AliITSClusterFinderSDD(seg1,res1,dig1,recp1);
+ rec1->SetCutAmplitude((int)thres);
+ ITS->SetReconstructionModel(1,rec1);
+ //rec1->Print();
+
+ // SSD
+
+ AliITSDetType *iDetType=ITS->DetType(2);
+ AliITSsegmentationSSD *seg2=(AliITSsegmentationSSD*)iDetType->GetSegmentationModel();
+ TClonesArray *dig2 = ITS->DigitsAddress(2);
+ AliITSClusterFinderSSD *rec2=new AliITSClusterFinderSSD(seg2,dig2);
+ ITS->SetReconstructionModel(2,rec2);
+
+//
+// Loop over events
+//
+ Int_t Nh=0;
+ Int_t Nh1=0;
+ for (int nev=0; nev<= evNumber2; nev++) {
+ Int_t nparticles = 0;
+ nparticles = gAlice->GetEvent(nev);
+ cout << "nev " << nev <<endl;
+ cout << "nparticles " << nparticles <<endl;
+ if (nev < evNumber1) continue;
+ if (nparticles <= 0) return;
+
+ AliITShit *itsHit;
+ AliITShit *itsHitPrev;
+ AliITSRecPoint *itsPnt = 0;
+ AliITSRawClusterSDD *itsClu = 0;
+
+ // Get Hit, Cluster & Recpoints Tree Pointers
+
+ TTree *TH = gAlice->TreeH();
+ Int_t nenthit=TH->GetEntries();
+ printf("Found %d entries in the Hit tree (must be one per track per event!)\n",nenthit);
+
+ ITS->GetTreeC(nev);
+ TTree *TC=ITS->TreeC();
+ Int_t nentclu=TC->GetEntries();
+ printf("Found %d entries in the Cluster tree (must be one per module per event!)\n",nentclu);
+
+ TTree *TR = gAlice->TreeR();
+ Int_t nentrec=TR->GetEntries();
+ printf("Found %d entries in the RecPoints tree\n",nentrec);
+
+ // Get Pointers to Clusters & Recpoints TClonesArrays
+
+ TClonesArray *ITSclu = ITS->ClustersAddress(1);
+ printf ("ITSclu %p \n",ITSclu);
+ // Int_t ncl = ITSclu->GetEntries();
+ // cout<<"ncluster ="<<ncl<<endl;
+ TClonesArray *ITSrec = ITS->RecPoints();
+ printf ("ITSrec %p \n",ITSrec);
+
+ // check recpoints
+
+ //Int_t nbytes = 0;
+ Int_t totpoints = 0;
+ Int_t totclust = 0;
+
+ // check hits
+
+ Int_t nmodules=0;
+
+ ITS->InitModules(-1,nmodules);
+ ITS->FillModules(nev,0,nmodules,"","");
+
+ TObjArray *fITSmodules = ITS->GetModules();
+
+ Int_t first0 = geom->GetStartDet(0); // SPD
+ Int_t last0 = geom->GetLastDet(0); // SPD
+ Int_t first1 = geom->GetStartDet(1); // SDD
+ Int_t last1 = geom->GetLastDet(1); // SDD
+ Int_t first2 = geom->GetStartDet(2); // SSD
+ Int_t last2 = geom->GetLastDet(2); // SSD
+
+ // For the SPD: first0 = 0, last0 = 239 (240 modules);
+ // for the SDD: first1 = 240, last1 = 499 (260 modules);
+ // for the SSD: first2 = 500, last2 = 2269 (1770 modules).
+
+ Int_t mod;
+ printf("det type %d first0, last0 %d %d \n",0,first0,last0);
+ printf("det type %d first1, last1 %d %d \n",1,first1,last1);
+ printf("det type %d first2, last2 %d %d \n",2,first2,last2);
+ Int_t negtrSDD = 0;
+ Int_t negtrSSD = 0;
+
+ for (mod=0; mod<last2+1; mod++) {
+ if(mod>max_modules)continue;
+ if(mod < first1) continue; // for the SDD/SSD only
+
+ AliITSmodule *Mod = (AliITSmodule *)fITSmodules->At(mod);
+ Int_t nhits = Mod->GetNhits();
+ Int_t layer = Mod->GetLayer();
+ Int_t ladder = Mod->GetLadder();
+ Int_t detector = Mod->GetDet();
+ Int_t hit;
+ Int_t parent;
+ Int_t dray;
+ Int_t hitstat;
+ Int_t partcode;
+ Int_t hitprim;
+ Float_t epart;
+ Float_t pathInSi=300;
+ Float_t xhit;
+ Float_t yhit;
+ Float_t zhit;
+ Float_t px;
+ Float_t py;
+ Float_t pz;
+ Float_t pmod;
+ Float_t xhit0 = 1e+7;
+ Float_t yhit0 = 1e+7;
+ Float_t zhit0 = 1e+7;
+
+ TTree *TR = gAlice->TreeR();
+ ITS->ResetClusters();
+ //cout << "CLUSTERS: get" << endl;
+ TC->GetEvent(mod);
+ //cout << "RECPOINTS: reset" << endl;
+ ITS->ResetRecPoints();
+ //cout << "RECPOINTS: get" << endl;
+ //TR->GetEvent(mod+1);
+ TR->GetEvent(mod);
+
+ Int_t nrecp = ITSrec->GetEntries();
+ totpoints += nrecp;
+ if (!nrecp) continue;
+
+ Int_t trackRecp[3];
+ Int_t itr;
+ Int_t TrackRecPoint;
+ Float_t PmodRecP;
+ Float_t signalRP;
+
+// cout <<" module,layer,ladder,detector,nrecp,nhits ="<<
+// mod<<","<<layer<<","<<ladder<<","<<detector<<","<<nrecp<<","<<nhits<< endl;
+ cout<<".";
+
+ // ---------- RecPoint signal analysis for the SDD/SSD ---------
+
+ for (Int_t pnt=0;pnt<nrecp;pnt++) { // recpoint loop
+
+ itsPnt = (AliITSRecPoint*)ITSrec->At(pnt);
+
+ Int_t RecPointPrim = 0;
+ if(!itsPnt) continue;
+
+ signalRP = itsPnt->GetQ();
+ trackRecp[0] = itsPnt->GetLabel(0);
+ trackRecp[1] = itsPnt->GetLabel(1);
+ trackRecp[2] = itsPnt->GetLabel(2);
+
+// cout<<"New Recp: pnt,signal,tr0,tr1,tr2 ="<<pnt
+// <<","<<signalRP<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl;
+
+ /*
+ if(trackRecp[0]<0&&trackRecp[1]<0&&trackRecp[2]<0) {
+cout<<"pnt,tr0,1,2 ="<<pnt<<","<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<endl;
+ if(mod<first2) negtrSDD += 1;
+ if(mod>=first2) negtrSSD += 1;
+ }
+ */
+
+ xhit0 = 1e+7;
+ yhit0 = 1e+7;
+ zhit0 = 1e+7;
+
+ // Hit loop
+ for (hit=0;hit<nhits;hit++) {
+
+ itsHit = (AliITShit*)Mod->GetHit(hit);
+
+ hitprim = 0;
+ dray = 0;
+ Int_t hitRecp = 0;
+ Int_t trackHit = itsHit->GetTrack();
+ hitstat = itsHit->GetTrackStatus();
+ zhit = 10000*itsHit->GetZL();
+ xhit = 10000*itsHit->GetXL();
+ yhit = 10000*itsHit->GetYL();
+ Float_t dEn = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
+
+ // cout<<" New hit: hit,track,hitstat,yhit,ehit ="<<hit<<","<<trackHit<<","<<hitstat<<","<<yhit<<","<<dEn<<endl;
+
+ // check the primary particle
+ partcode = itsHit->GetParticle()->GetPdgCode();
+ parent = itsHit->GetParticle()->GetFirstMother();
+ if(parent < 0) hitprim = 1; // primery particle
+
+
+
+ // select the primery hits with the track number equaled to
+ // the one of the RecPoint
+ for(itr=0;itr<3;itr++) {
+ if(trackRecp[itr] == trackHit) hitRecp = 1;
+ if(trackRecp[itr] == trackHit && hitprim == 1) {
+ RecPointPrim = 1;
+ TrackRecPoint = trackRecp[itr];
+ }
+ if(trackRecp[itr] == trackHit && hitprim == 0) trackRecp[itr]=-100;
+ }
+
+ if(hitRecp != 1) continue; // hit doesn't correspond to the recpoint
+
+ //cout<<"Hit Ok: trackRecp0,1,2,hitRecp,RecPointPrim,TrackRecPoint ="<<trackRecp[0]<<","<<trackRecp[1]<<","<<trackRecp[2]<<","<<hitRecp<<","<<RecPointPrim<<","<<TrackRecPoint<<endl;
+
+ // if(hitstat == 66 && yhit < -136.) { // entering hit
+ if(hitstat == 66 && hitprim == 1) { // entering hit
+ xhit0 = xhit;
+ yhit0 = yhit;
+ zhit0 = zhit;
+ }
+ // cout<<" xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl;
+
+ if(hitstat == 66) continue; // Take only the hits with the not
+ // zero energy
+
+ if(xhit0 > 9e+6 || zhit0 > 9e+6) {
+ // cout<<"default xhit0,zhit0 ="<<xhit0<<","<<zhit0<<","<<endl;
+ continue;
+ }
+
+ // check the particle code (type)
+ // Int_t parent = itsHit->GetParticle()->GetFirstMother();
+ //partcode = itsHit->GetParticle()->GetPdgCode();
+
+ // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
+ // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
+
+ /*
+ px = itsHit->GetPXL(); // momentum for the hit
+ py = itsHit->GetPYL();
+ pz = itsHit->GetPZL();
+ pmod = 1000*sqrt(px*px+py*py+pz*pz);
+ */
+
+ pmod = itsHit->GetParticle()->P(); // momentum at the vertex
+ pmod *= 1.0e+3;
+
+ // cout<<"track,partcode,pmod,hitprim ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<endl;
+
+ if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
+ // at p < 6 MeV/c
+
+ // find the primery particle path length in Si and pmod (MeV/c)
+ if((hitstat == 68 || hitstat == 33) && hitprim == 1) {
+ pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
+ PmodRecP = itsHit->GetParticle()->P();
+ PmodRecP *= 1.0e+3;
+ // cout<<"path in Si="<<pathInSi<<endl;
+
+ Float_t Signal = signalRP*(300/pathInSi)/38;
+ if(Signal<0.5) {
+ //cout<<"track,partcode,pmod,hitprim,hitstat,Signal,dEn ="<<trackHit<<","<<partcode<<","<<pmod<<","<<hitprim<<","<<hitstat<<","<<Signal<<","<<dEn<<endl;
+ }
+ }
+ } // hit loop
+
+ if(RecPointPrim == 1) {
+
+ //cout<<" SDD/SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
+
+ Int_t xpartcode=gAlice->Particle(TrackRecPoint)->GetPdgCode();
+
+ if(mod < first2) { // SDD
+ signalRP *= (280./pathInSi); // 280 microns of Si thickness
+ signalSDD->Fill(signalRP);
+ signalSDDmip->Fill(signalRP/280.); // ADC/280 -> mips
+ InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/280. ,
+ (Float_t)PmodRecP,(Int_t)xpartcode);
+ //cout<<" SDD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
+ if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SDD ="<<pathInSi<<endl;
+
+ //if(signalRP/280 < 5) cout<<" small signalSDDmip, path ="<<signalRP/280<<","<<pathInSDD<<endl;
+ }else{ // SSD
+ signalRP *= (300/pathInSi); // 300 microns of Si thickness
+ //cout<<" SSD RecPoints: lay,lad,det,track,pmod,signal,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackRecPoint<<","<<PmodRecP<<","<<signalRP<<","<<pathInSi<<endl;
+ signalSSD->Fill(signalRP);
+ signalSSDmip->Fill(signalRP/38.); // ADC/38 -> mips
+ InPID(pid,(Int_t)nev,(Int_t)TrackRecPoint, signalRP/38. ,
+ (Float_t)PmodRecP,(Int_t)xpartcode);
+ if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathInSi in SSD ="<<pathInSi<<endl;
+
+ if(signalRP/38 < 0.5) cout<<" small signalSSD (mip), path, module ="<<signalRP/38<<","<<pathInSi<<mod<<endl;
+ }
+ } // primery particle
+ } // pnt, recpoint loop
+
+ // dEdX hit analysis for the SDD/SSD
+
+ Int_t track = -100;
+ Int_t trackPrev = -100;
+ parent = -100;
+ Int_t parentPrev = -100;
+ Int_t flagprim = 0;
+ Int_t TrackPart;
+ epart = 0;
+ pathInSi = 1.0e+7;
+ Float_t PmodPart;
+ zhit0 = 1.0e+7;
+ xhit0 = 1.0e+7;
+ yhit0 = 1.0e+7;
+
+ /*
+ if(mod <= 259) {
+ cout<<" SDD hits: nhits ="<<nhits<<endl;
+ }else{
+ cout<<" SSD hits: nhits ="<<nhits<<endl;
+ }
+ */
+
+ for (hit=0;hit<nhits;hit++) {
+
+ itsHit = (AliITShit*)Mod->GetHit(hit);
+ if(hit>0) itsHitPrev = (AliITShit*)Mod->GetHit(hit-1);
+ hitprim = 0;
+ Int_t hitprimPrev = 0;
+ track = itsHit->GetTrack();
+ if(hit > 0) Int_t trackPrev = itsHitPrev->GetTrack();
+ dray = 0;
+ hitstat = itsHit->GetTrackStatus();
+
+ zhit = 10000*itsHit->GetZL();
+ xhit = 10000*itsHit->GetXL();
+ yhit = 10000*itsHit->GetYL();
+ Float_t ehit = 1.0e+6*itsHit->GetIonization(); // hit energy, KeV
+
+ // cout<<"New hit,lay,hitstat,yhit,ehit ="<<hit<<","<<layer<<","<<hitstat<<","<<yhit<<","<<ehit<<endl;
+ // cout<<"befor: track,trackPrev ="<<track<<","<<trackPrev<<endl;
+
+ parent = itsHit->GetParticle()->GetFirstMother();
+ if(hit>0) Int_t parentPrev = itsHitPrev->GetParticle()->GetFirstMother();
+ partcode = itsHit->GetParticle()->GetPdgCode();
+
+ // partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - i+
+ // 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
+
+ px = itsHit->GetPXL();
+ py = itsHit->GetPYL();
+ pz = itsHit->GetPZL();
+ pmod = 1000*sqrt(px*px+py*py+pz*pz);
+
+ // cout<<"partcode,pmod,parent,parentPrev,hitprim,flagprim ="<<partcode<<","<<pmod<<","<<parent<<","<<parentPrev<<","<<hitprim<<","<<flagprim<<endl;
+
+
+ if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
+ // at p < 6 MeV/c
+
+
+ if(parent < 0 && parent > -90) hitprim += 1;
+ if(parentPrev < 0 && parentPrev >-90) hitprimPrev += 1;
+ // hitprim=1 for the primery particles
+ // if(hit==0&&hitprim==1&&hitstat==66) {
+ if(hitprim==1&&hitstat==66) {
+ zhit0 = zhit;
+ xhit0 = xhit;
+ yhit0 = yhit;
+ }
+
+ if((hitstat == 68 || hitstat == 33) && hitprim == 1) {
+ pathInSi = TMath::Sqrt((xhit0-xhit)*(xhit0-xhit)+(yhit0-yhit)*(yhit0-yhit)+(zhit0-zhit)*(zhit0-zhit));
+ TrackPart = track;
+ PmodPart = itsHit->GetParticle()->P();
+ PmodPart *= 1.0e+3;
+ // cout<<"xhit0,yhit0,zhit0,pathInSi ="<<xhit0<<","<<yhit0<<","<<zhit0<<","<<pathInSi<<endl;
+ }
+
+ if(hitprim == 0) track = parent;
+ if(hitprimPrev == 0) trackPrev = parentPrev;
+ // cout<<"after: track,trackPrev ="<<track<<","<<trackPrev<<endl;
+
+ if(hit > 0 && track == trackPrev) epart += ehit;
+ if((hit > 0 && track == trackPrev) && (hitprim==1)) flagprim +=1;
+
+ // cout<<"new hitprim, flagprim ="<<hitprim<<","<<flagprim<<endl;
+
+ if((hit > 0 && track != trackPrev) || (hit == nhits-1)) {
+ if(flagprim > 0) {
+ if(layer > 2 && layer < 5) {
+ if(epart < 40) cout<<" small SDD dedx ="<<epart<<endl;
+ if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSDD ="<<pathInSi<<endl;
+ if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) {
+ dedxSDD->Fill(epart);
+ epart *= (280/pathInSi);
+ dedxSDDnorm->Fill(epart);
+ dedxSDDmip->Fill(epart/75);
+ //InPID(pid,(Int_t)nev,(Int_t)track, epart/75. ,(Float_t)pmod,(Int_t)partcode);
+ pathSDD->Fill(pathInSi);
+ // cout<<"SDD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl;
+ }
+ }
+ if(layer > 4) {
+ if(epart < 40) cout<<" small SSD dedx ="<<epart<<endl;
+ if(pathInSi < 100 || pathInSi > 2000) cout<<" No normal pathSSD ="<<pathInSi<<endl;
+ if(epart > 40 && pathInSi > 100 && pathInSi < 1.0e+5) {
+ dedxSSD->Fill(epart);
+ epart *= (280/pathInSi);
+ dedxSSDnorm->Fill(epart);
+ dedxSSDmip->Fill(epart/80);
+ //InPID(pid,(Int_t)nev,(Int_t)track, epart/80. ,(Float_t)pmod,(Int_t)partcode);
+ pathSSD->Fill(pathInSi);
+ //cout<<"SSD hit: lay,lad,det,track,pmod,epart,path ="<<layer<<","<<ladder<<","<<detector<<","<<TrackPart<<","<<PmodPart<<","<<epart<<","<<pathInSi<<endl;
+ }
+ }
+ flagprim=0;
+ epart = 0;
+ }else{
+ //cout<<"No!, epart for the secondary"<<endl;
+ epart = 0;
+ }
+ }
+ } // SDD/SSD hit loop
+ } // module loop
+ //cout<<"negtrSDD, negtrSSD ="<<negtrSDD<<","<<negtrSSD<<endl;
+ } // event loop
+
+ cout<<endl;
+
+ TFile fhistos("dEdX_his.root","RECREATE");
+
+ dedxSDD->Write();
+ dedxSDDnorm->Write();
+ dedxSDDmip->Write();
+ dedxSSD->Write();
+ dedxSSDnorm->Write();
+ dedxSSDmip->Write();
+ signalSDD->Write();
+ signalSDDmip->Write();
+ signalSSD->Write();
+ signalSSDmip->Write();
+ pathSDD->Write();
+ pathSSD->Write();
+
+ fhistos.Close();
+ cout<<"!!! Histogramms and ntuples were written"<<endl;
+ // ------------------------------------------------------------
+
+ TCanvas *c1 = new TCanvas("c1","ITS clusters",400,10,600,700);
+ c1->Divide(2,2);
+
+
+// c1->cd(1);
+// gPad->SetFillColor(33);
+// signalSDD->SetFillColor(42);
+// signalSDD->Draw();
+
+// c1->cd(2);
+// gPad->SetFillColor(33);
+// signalSDDmip->SetFillColor(42);
+// signalSDDmip->Draw();
+
+// c1->cd(3);
+// gPad->SetFillColor(33);
+// dedxSDDnorm->SetFillColor(46);
+// dedxSDDnorm->Draw();
+
+// c1->cd(4);
+// gPad->SetFillColor(33);
+// dedxSDDmip->SetFillColor(46);
+// dedxSDDmip->Draw();
+
+ c1->cd(1);
+ gPad->SetFillColor(33);
+ signalSSD->SetFillColor(42);
+ signalSSD->Draw();
+
+ c1->cd(2);
+ gPad->SetFillColor(33);
+ signalSSDmip->SetFillColor(42);
+ signalSSDmip->Draw();
+
+ c1->cd(3);
+ gPad->SetFillColor(33);
+ dedxSSDnorm->SetFillColor(46);
+ dedxSSDnorm->Draw();
+
+ c1->cd(4);
+ gPad->SetFillColor(33);
+ dedxSSDmip->SetFillColor(46);
+ dedxSSDmip->Draw();
+
+
+ /*
+ c1->cd(1);
+ gPad->SetFillColor(33);
+ dedxSDD->SetFillColor(42);
+ dedxSDD->Draw();
+
+ c1->cd(2);
+ gPad->SetFillColor(33);
+ pathSDD->SetFillColor(42);
+ pathSDD->Draw();
+
+ c1->cd(3);
+ gPad->SetFillColor(33);
+ dedxSSD->SetFillColor(46);
+ dedxSSD->Draw();
+
+ c1->cd(4);
+ gPad->SetFillColor(33);
+ pathSSD->SetFillColor(46);
+ pathSSD->Draw();
+ */
+
+// if(!pid)pid->Tab();
+ cout<<"END test for clusters and hits "<<endl;
+
+}
+//============================ InPID() =========================================
+int totpid;
+void InPID (AliITSPid *pid=0,int nev,
+ Int_t track,
+ Float_t signal,
+ Float_t pmom,
+ Int_t pcode)
+{
+// Includes recp in PID table.
+ if(pid){
+ //cout<<" In pid: track,sig,pmod,pcode="<<track<<","<<signal<<","<<pmom<<","<<pcode<<endl;
+ if( abs(pcode)==321 || abs(pcode)==211 ||abs(pcode)==2212 ||abs(pcode)==11 )
+ {
+ pid->SetEdep(10000*nev+track,signal);
+ pid->SetPmom(10000*nev+track,pmom);
+ pid->SetPcod(10000*nev+track,abs(pcode));
+ }
+ }
+ totpid++;
+}
+//=======================
+
+
+
+
+
+
--- /dev/null
+//=========================================
+int HMax=80,PMax=6.;
+int NStat=15000;
+Float_t pmin,pmax;
+
+char tit[]="--------------------------";
+TH1F pions=TH1F("Qpi",tit,100,0.5,PMax);
+TH1F *qpi =&pions;
+
+TH2F qplotP, qplotKa, qplotPi, qplotE;
+TH1F *q1plot;
+TH2F *qplot;
+void dedxhis(Int_t pcod=0,Float_t pmin=0,Float_t pmax=1300)
+{
+ if(!q1plot)q1plot = new TH1F("Qhis","Particle charge",100,0.5,PMax);
+ if(!qplot){ qplot = new TH2F("Qtrm","Qtrm vs Pmom",100,0,1300,100,0,13);
+ // TH2F qplotP, qplotKa, qplotPi, qplotE;
+ qplot->Copy(qplotP);
+ qplotP.SetMarkerStyle(8); qplotP.SetMarkerColor(kBlack); qplotP.SetMarkerSize(.2);
+ qplot->Copy(qplotKa);
+ qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.2);
+ qplot->Copy(qplotPi);
+ qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); qplotPi.SetMarkerSize(.2);
+ qplot->Copy(qplotE);
+ qplotE.SetMarkerStyle(8); qplotE.SetMarkerColor(kGreen); qplotE.SetMarkerSize(.2);
+ }else{
+ qplotP.Reset();
+ qplotPi.Reset(); qplotKa.Reset(); qplotE.Reset();
+ }
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ if(pcod==0)
+ { if(tabv(0)>=1 ){
+ //qplot->Fill( Pmom,Qtr );
+ if(tabv(11)==2212)qplotP.Fill( Pmom,Qtr );
+ if(tabv(11)== 321)qplotKa.Fill( Pmom,Qtr );
+ if(tabv(11)== 211)qplotPi.Fill( Pmom,Qtr );
+ if(tabv(11)== 11)qplotE.Fill( Pmom,Qtr );
+ }
+ }
+ else { if(tabv(11)==pcod && tabv(0)>=1 )qplot->Fill( Pmom,Qtr );}
+
+ //if( tabv(0)>=2 && tabv(11)==pcod )qplot->Fill( Pmom,Qtr );
+ //if(tabv(0)>=1)qplot->Fill( Pmom,Qtr );
+ if(Pmom>pmin&&Pmom<pmax)q1plot->Fill(Qtr);
+ }
+// ------------------------------------------------------------
+// c1=new TCanvas("c1"," ",10,10,700,500);
+// pad1 =new TPad("p1"," " ,0 , 0,0,1,21);
+// pad1 =new TPad("p1"," " ,0 , 0,.3,1,21);
+// pad2 =new TPad("p2"," " ,.35, 0, 1,1,0);
+// pad2->Draw(); pad2->cd();
+//.....................
+ //if(pcod==0){ qplotP.Draw(); qplotKa.Draw("same"); qplotPi.Draw("same"); qplotE.Draw("same"); }
+
+ if(pcod==0){
+
+ //qplotP->GetXaxis()->SetTitleSize(0.05);
+ //qplotP->GetYaxis()->SetTitleSize(0.05);
+ gStyle->SetOptStat(0);
+ qplotP->SetXTitle("(Mev/c)");
+ qplotP->SetYTitle("(mips)");
+ qplotP.Draw(); qplotKa.Draw("same"); qplotPi.Draw("same");
+ qplotE.Draw("same");
+ TText *text = new TText(800.,11.,"Protons");
+ text->SetTextSize(0.05);
+ text->SetTextColor(kBlack);
+ text->Draw();
+ TText *text = new TText(800.,10.,"Kaons");
+ text->SetTextSize(0.05);
+ text->SetTextColor(kRed);
+ text->Draw();
+ TText *text = new TText(800.,9.,"Pions");
+ text->SetTextSize(0.05);
+ text->SetTextColor(kBlue);
+ text->Draw();
+ TText *text = new TText(800.,8.,"Electrons");
+ text->SetTextSize(0.05);
+ text->SetTextColor(kGreen);
+ text->Draw();
+
+ }
+ else{
+ qplot->Draw();}
+ 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();
+ }
+ //Draw kaons-protons cuts.
+ gStyle->SetLineWidth(1.);
+}
+
+void qhispi(Float_t pmin=0,Float_t pmax=1300)
+{
+ char tit[]="--------------------------";
+ sprintf(tit,"%s%.1f %.1f","Pmom ",pmin,pmax);
+// TH1F *qpi = new TH1F("Qpi",tit,100,0.5,PMax);
+// TH1F *nhit = new TH1F("Nhit",tit,100,0,6);
+ TH2F *nhit = new TH2F("Nhit",tit,100,0,6,100,0,1300);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ if(Pmom>pmin&&Pmom<pmax&&abs(tabv(11))==211)qpi->Fill(Qtr);
+// if(tabv(11)==2212)nhit->Fill(tabv(0));
+ if(tabv(11)==211)nhit->Fill(tabv(0),Pmom);
+ }
+// ------------------------------------------------------------
+ int pistat=qpi->GetEntries();
+ //qpi->SetMaximum(pistat+1);
+ qpi->SetMaximum(HMax);
+ qpi->SetFillColor(42);
+ qpi->Fit("gaus");
+ TF1 *fun=qpi->GetFunction("gaus");
+ fun->SetLineWidth(.1);
+ qpi->SetLineWidth(.1);
+ qpi->Draw();
+//---------------------------
+ if(0){
+ gaus1=new TF1("g1","gaus",.5,2);
+ qpi->SetMaximum(-1111);
+ gaus1->SetParameter(0, qpi->GetMaximum() );
+ qpi->SetMaximum(HMax);
+ gaus1->SetParameter(1,1);
+ gaus1->SetParameter(2,.12);
+ gaus1->SetLineWidth(.1);
+ gaus1->Draw("same");
+ }
+}
+
+void qhiska(Float_t pmin=0,Float_t pmax=1300)
+{
+ TH1F *qka = new TH1F("Qka","Particle charge",100,0.5,PMax);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+// qplot->Fill( Pmom,Qtr );
+ if(Pmom>pmin&&Pmom<pmax&&tabv(11)==321)qka->Fill(Qtr);
+ }
+// ------------------------------------------------------------
+ qka->SetFillColor(0);
+ qka->SetLineColor(kRed);
+//...................
+ qka->Fit("gaus","0");
+ TF1 *fun=qka->GetFunction("gaus");
+ fun->SetLineWidth(.1);
+ qka->SetLineWidth(.1);
+ qka->Draw("same");
+//...................
+// qka->Draw("same");
+// gaus_k(pmin,pmax,qka);
+
+}
+
+void qhisp(Float_t pmin=0,Float_t pmax=1300)
+{
+ TH1F *qpr = new TH1F("Qpr","Particle charge",100,0.5,PMax);
+// -------------------------------------------------------------
+ qpr->Clear();
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+// qplot->Fill( Pmom,Qtr );
+ if(Pmom>pmin&&Pmom<pmax&&tabv(11)==2212)qpr->Fill(Qtr);
+ }
+// ------------------------------------------------------------
+ qpr->SetFillColor(16);
+//...................
+ qpr->Fit("gaus","0");
+ TF1 *fun=qpr->GetFunction("gaus");
+ fun->SetLineWidth(.1);
+ qpr->SetLineWidth(.1);
+ qpr->Draw("same");
+//...................
+// qpr->Draw("same");
+// gaus_p(pmin,pmax,qpr);
+}
+//---------------------------------------
+
+void gaus_k(Float_t pmin=450,Float_t pmax=470,TH1F *qka){
+ Float_t xmean,xsig;
+ if(pmin==410&&pmax==470){xmean=pid->cut[5][3]; xsig=pid->cut[5][4];}
+ else
+ if(pmin==470&&pmax==530){xmean=pid->cut[6][3]; xsig=pid->cut[6][4];}
+ else
+ if(pmin==730&&pmax==830){xmean=pid->cut[10][3]; xsig=pid->cut[10][4];}
+ else
+ if(pmin==830&&pmax==930){xmean=pid->cut[11][3]; xsig=pid->cut[11][4];}
+ else
+ if(pmin==930&&pmax==1030){xmean=pid->cut[12][3]; xsig=pid->cut[12][4];}
+ else{return;}
+ gaus1=new TF1("g1","gaus",xmean-3*xsig,xmean+3*xsig);
+ qka->SetMaximum(-1111);
+ gaus1->SetParameter(0, qka->GetMaximum() );
+ qka->SetMaximum(HMax);
+ gaus1->SetParameter(1,xmean);
+ gaus1->SetParameter(2,xsig);
+ gaus1->SetLineWidth(.1);
+ gaus1->Draw("same");
+}
+//---------------------------------------
+void gaus_p(Float_t pmin=730,Float_t pmax=830,TH1F *qp){
+ Float_t xmean,xsig;
+ if(pmin==730&&pmax==830){xmean=pid->cut[10][5]; xsig=pid->cut[10][6];}
+ else
+ if(pmin==830&&pmax==930){xmean=pid->cut[11][5]; xsig=pid->cut[11][6];}
+ else
+ if(pmin==930&&pmax==1030){xmean=pid->cut[12][5]; xsig=pid->cut[12][6];}
+ else{return;}
+ gaus1=new TF1("g1","gaus",xmean-3*xsig,xmean+3*xsig);
+ qp->SetMaximum(-1111);
+ gaus1->SetParameter(0, qp->GetMaximum() );
+ qp->SetMaximum(HMax);
+ gaus1->SetParameter(1,xmean);
+ gaus1->SetParameter(2,xsig);
+ gaus1->SetLineWidth(.1);
+ gaus1->Draw("same");
+}
+
+//----b.b.---------------------------------
+void fitpi(Float_t pmin=0,Float_t pmax=1300,char *tit)
+{
+ cout<<"fitpi: NStat, PMax, pmin, pmax ="<<NStat<<","<<PMax<<","<<pmin<<","<<pmax<<endl;
+ TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ if(Pmom>pmin&&Pmom<pmax&&tabv(11)==211)q1fit->Fill(Qtr);
+ }
+ q1fit->SetFillColor(0);
+ q1fit->SetLineColor(kRed);
+ q1fit->Draw();
+ q1fit->Fit("gaus");
+}
+
+//---------------------------------------
+void fitka(Float_t pmin=0,Float_t pmax=1300,char *tit)
+{
+ //TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
+ TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,10.);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ if(Pmom>pmin&&Pmom<pmax&&tabv(11)==321)q1fit->Fill(Qtr);
+ }
+ q1fit->SetFillColor(0);
+ q1fit->SetLineColor(kRed);
+ q1fit->Draw();
+ q1fit->Fit("gaus");
+}
+//---------------------------------------
+void fitp(Float_t pmin=0,Float_t pmax=1300,char *tit)
+{
+ TH1F *q1fit = new TH1F("Qfit",tit,100,0.5,PMax);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ if(Pmom>pmin&&Pmom<pmax&&tabv(11)==2212)q1fit->Fill(Qtr);
+ }
+ q1fit->SetFillColor(0);
+ q1fit->Draw();
+ q1fit->Fit("gaus");
+}
+//---------------------------------------
+void effall(){
+
+ eff(211,211,"",3);
+ eff(321,321,"same",2);
+ eff(2212,2212,"same",1);
+}
+//---------------------------------------
+void eff(int pc=211,int pc2=211,char *opt="",int color=2)
+{
+ const Int_t HSize=100;
+ TH1F *efpi = new TH1F("Effpi","Eff of PID",HSize,0.025,1400);
+ TH1F *pmom = new TH1F("Pmom" ,"Eff of PID",HSize,0.025,1400);
+ TH1F *heff = new TH1F("Eff%" ,"Eff of PID",HSize,0.025,1400);
+// -------------------------------------------------------------
+ Float_t Qtr,Pmom;
+ Int_t Pcode;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Qtr=tabv(6);Pmom=tabv(10);
+ Pcode= ( pid->GetPcode(Qtr,Pmom/1000.) );
+ if(tabv(0)>=4){
+ if(Pcode==pc&&tabv(11)==pc2) efpi->Fill(Pmom);
+ if(tabv(11)==pc2) pmom->Fill(Pmom);
+ }
+ }
+// ------------------------------------------------------------
+ for(int i=1;i<=HSize;i++){
+ if( pmom->fArray[i] > 0 )
+ heff->fArray[i] = color * ( efpi->fArray[i] )/( pmom->fArray[i] );
+ }
+ heff->SetLineColor(color);
+ if(color==1)heff->SetLineWidth(2)else heff->SetLineWidth(.5);
+ heff->Draw(opt);
+return;
+ pmom->SetFillColor(16); pmom->Draw();
+ efpi->SetFillColor(0); //42
+ efpi->Draw("same");
+}
+//---------------------------------------
+Int_t NRange=0;
+ Float_t xpmin[]={410,470,730,830,930};
+ Float_t xpmax[]={470,530,830,930,1030};
+void all(){
+ PMax=3; HMax=500; NStat=99000;
+ Float_t pmin,pmax;
+ pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
+ qhispi(pmin,pmax); qhisp(pmin,pmax); qhiska(pmin,pmax);
+}
+//---------------------------------------
+void fitkall(){
+ PMax=3;
+ Float_t pmin,pmax;
+ //pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
+ pmin=200;pmax=300; // b.b.
+ char str[]=" ";
+ sprintf(str,"Kaons %d - %d MeV/c",pmin,pmax);
+ gStyle->SetOptFit();
+ fitka(pmin,pmax,str);
+}
+
+//--b.b.-------------------------------------
+void fitpiall(){
+ PMax=3;
+ NRange=3; // b.b.
+ Float_t pmin,pmax;
+ cout<<"fitpiall: NRange ="<<NRange<<endl;
+ //pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=0;
+ pmin=100;pmax=200; // b.b.
+ char str[]=" ";
+ sprintf(str,"Pions %d - %d MeV/c",pmin,pmax);
+ gStyle->SetOptFit();
+ cout<<"fitpiall: pmin, pmax ="<<pmin<<","<<pmax<<endl;
+ fitpi(pmin,pmax,str);
+}
+
+//---------------------------------------
+void fitpall(){
+ PMax=3;
+ Float_t pmin,pmax;
+ if(NRange==0)NRange=2;
+ pmin=xpmin[NRange];pmax=xpmax[NRange];NRange++;if(NRange==5)NRange=2;
+ char str[]=" ";
+ sprintf(str,"Protons %d - %d MeV/c",pmin,pmax);
+ gStyle->SetOptFit();
+ fitp(pmin,pmax,str);
+}
+//-------------
+void newcuts(){
+//
+
+ pid->SetCut(3,.3, 0 , 2.5, 2.5 , 9, 9, 10 ); //200-300
+
+ pid->SetCut(5,.47, 1 , 0.12, 1.98 , 1.17 , 2.5, 10 );//410-470
+ pid->SetCut(6,.53, 1 , 0.12, 1.73 , 0.15 , 2.5, 10 );//470-530
+ pid->SetCut(7,.59, 0 , 0, 1.18 , 1.125 , 2.5, 10 );//530-590
+
+ pid->SetCut(8,.65, 0 , 0, 1.18, 1.125 , 2.3, 10 );//590-650
+}
+//---------------------------------------
+void qhisall(Float_t pmin=0.25,Float_t pmax=.700)
+//void qhisall()
+{
+ qhispi(pmin,pmax);
+ qhisp(pmin,pmax);
+ qhiska(pmin,pmax);
+}
+//--------------------------------------
+void pcode(){
+ TH1F *his = new TH1F("pcode","tit",100,10,2300);
+Float_t Pcod;
+for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ Pcod=tabv(11); if(Pcod>0) his->Fill(Pcod);
+}
+his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
+}
+
+void signal(){
+ TH1F *hisR = new TH1F("Rec signal","tit",100,0,13);
+ TH1F *hisH = new TH1F("Hit signal","tit",100,0,13);
+ Float_t xx;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ xx=tabv(1); if(tabv(0)>0) hisR->Fill(xx);
+ }
+ hisR->SetFillColor(0);hisR->SetLineColor(kRed);hisR->Draw();
+ hisH->SetFillColor(0);hisH->SetLineColor(kBlue);hisH->Draw("same");
+}
+//-------------------------------------
+
+void pmom(){
+ TH1F *his = new TH1F("pmom","tit",100,0,1000);
+ Float_t xx;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ xx=tabv(10); if(tabv(0)>0) his->Fill(xx);
+ }
+ his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
+}
+
+// Print next track
+int jtrack=0;
+void track(){
+ if(jtrack==NStat)jtrack=0;
+ for(Int_t i=jtrack;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ if(tabv(0)>0)
+ { jtrack=i;break; }
+ }
+cout<<"Track No "<<jtrack<<endl;
+pid->Print(jtrack++);
+}
+
+// Fill histogram for track number
+//--------------------------------
+void tracks(){
+ TH1F *his = new TH1F("tracks","tit",100,0,15000);
+ Float_t xx;
+ for(Int_t i=0;i<NStat;i++){
+ TVector tabv(*( pid->GetVec(i) ));
+ if(tabv(0)>0)
+ { his->Fill(i); }
+ }
+ his->SetFillColor(0);his->SetLineColor(kRed);his->Draw();
+}
+// Fill pid table with reconstructed tracks
+
+
+#ifndef __CINT__
+ #include <iostream.h>
+ #include <fstream.h>
+
+ #include "AliRun.h"
+ #include "AliITS.h"
+ #include "AliITSgeom.h"
+ #include "AliITStrackerV2.h"
+ #include "AliITStrackV2.h"
+ #include "AliITSclusterV2.h"
+
+ #include "TFile.h"
+ #include "TTree.h"
+ #include "TH1.h"
+ #include "TObjArray.h"
+ #include "TStyle.h"
+ #include "TCanvas.h"
+ #include "TLine.h"
+ #include "TText.h"
+#endif
+
+struct GoodTrackITS {
+ Int_t lab;
+ Int_t code;
+ Float_t px,py,pz;
+ Float_t x,y,z;
+};
+
+
+
+void filltab_tracks(){
+
+ cerr<<"Filling of track table...\n";
+
+ Int_t event=0;
+Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event);
+
+ const Int_t MAX=15000;
+ Int_t nentr=0; TObjArray tarray(2000);
+ {/* Load tracks */
+ TFile *tf=TFile::Open("AliITStracksV2.root");
+ if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
+ char tname[100]; sprintf(tname,"TreeT_ITS_%d",event);
+ TTree *tracktree=(TTree*)tf->Get(tname);
+ if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ nentr=(Int_t)tracktree->GetEntries();
+ for (Int_t i=0; i<nentr; i++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+ tarray.AddLast(iotrack);
+ }
+ delete tracktree; //Thanks to Mariana Bondila
+ tf->Close();
+ }
+
+ /* Generate a list of "good" tracks */
+ GoodTrackITS gt[MAX];
+ Int_t ngood=0;
+ Float_t ConvRadAng=180./TMath::Pi();
+ ifstream in("good_tracks_its");
+ if (in) {
+ cerr<<"Reading good tracks...\n";
+ while (in>>gt[ngood].lab>>gt[ngood].code>>
+ gt[ngood].px>>gt[ngood].py>>gt[ngood].pz>>
+ gt[ngood].x >>gt[ngood].y >>gt[ngood].z) {
+ ngood++;
+ cerr<<ngood<<'\r';
+ if (ngood==MAX) {
+ cerr<<"Too many good tracks !\n";
+ break;
+ }
+ }
+ if (!in.eof()) cerr<<"Read error (good_tracks_its) !\n";
+ } else {
+ cerr<<"Marking good tracks (this will take a while)...\n";
+ ngood=good_tracks_its(gt,MAX,event);
+ ofstream out("good_tracks_its");
+ if (out) {
+ for (Int_t ngd=0; ngd<ngood; ngd++)
+ out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '
+ <<gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '
+ <<gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
+ } else cerr<<"Can not open file (good_tracks_its) !\n";
+ out.close();
+ }
+ cerr<<"Number of good tracks : "<<ngood<<endl;
+ TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
+ TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
+ TH1F *hpt=new TH1F("hpt","Relative Pt resolution",30,-10.,10.);
+ hpt->SetFillColor(2);
+ TH1F *hmpt=new TH1F("hmpt","Transverse impact parameter",30,-300,300);
+ hmpt->SetFillColor(6);
+ TH1F *hz=new TH1F("hz","Longitudinal impact parameter",30,-300,300);
+ //hmpt->SetFillColor(6);
+
+ AliITStrackV2 *trk=(AliITStrackV2*)tarray.UncheckedAt(0);
+ Double_t pmin=0.1*(100/0.299792458/0.2/trk->GetConvConst());
+ Double_t pmax=6.0+pmin;
+
+ TH1F *hgood=new TH1F("hgood","Good tracks",30,pmin,pmax);
+ TH1F *hfound=new TH1F("hfound","Found tracks",30,pmin,pmax);
+ TH1F *hfake=new TH1F("hfake","Fake tracks",30,pmin,pmax);
+ TH1F *hg=new TH1F("hg","",30,pmin,pmax); //efficiency for good tracks
+ hg->SetLineColor(4); hg->SetLineWidth(2);
+ TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,pmin,pmax);
+ hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
+
+ TH1F *hptw=new TH1F("hptw","Weghted pt",30,pmax,pmin);
+
+ while (ngood--) {
+ Int_t lab=gt[ngood].lab, tlab=-1;
+ Double_t pxg=gt[ngood].px, pyg=gt[ngood].py, pzg=gt[ngood].pz;
+ Double_t ptg=TMath::Sqrt(pxg*pxg+pyg*pyg);
+ Double_t pg=1000.*TMath::Sqrt(pxg*pxg+pyg*pyg+pzg*pzg); // b.b.
+
+
+ if (ptg<pmin) continue;
+
+ hgood->Fill(ptg);
+
+ AliITStrackV2 *track=0;
+ Int_t j;
+ for (j=0; j<nentr; j++) {
+ track=(AliITStrackV2*)tarray.UncheckedAt(j);
+ tlab=track->GetLabel();
+ if (lab==TMath::Abs(tlab)) break;
+ }
+ if (j==nentr) {
+ cerr<<"Track "<<lab<<" was not found !\n";
+ continue;
+ }
+ track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ track->PropagateToVertex();
+
+ if (lab==tlab) hfound->Fill(ptg);
+ else { hfake->Fill(ptg); cerr<<lab<<" fake\n";}
+
+ Double_t xv,par[5]; track->GetExternalParameters(xv,par);
+ Float_t phi=TMath::ASin(par[2]) + track->GetAlpha();
+ if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+ if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+ Float_t lam=TMath::ATan(par[3]);
+ Float_t pt_1=TMath::Abs(par[4]);
+ Double_t phig=TMath::ATan2(pyg,pxg);
+ //hp->Fill((phi - phig)*1000.);
+
+ Double_t lamg=TMath::ATan2(pzg,ptg);
+ //hl->Fill((lam - lamg)*1000.);
+
+ Double_t d=10000*track->GetD();
+ //hmpt->Fill(d);
+
+ //hptw->Fill(ptg,TMath::Abs(d));
+
+ Double_t z=10000*track->GetZ();
+ //hz->Fill(z);
+
+
+
+// ------------------------------------------------ b.b. -----
+ Int_t nev=0;
+ Float_t mom=1000.*(1./(pt_1*TMath::Cos(lam)));
+ Float_t dedx=track->GetdEdx();
+ //hep->Fill(mom,dedx,1.);
+ Int_t pcode=gt[ngood].code;
+ pid->SetEdep(10000*nev+ngood,dedx);
+ pid->SetPmom(10000*nev+ngood,mom);
+ pid->SetPcod(10000*nev+ngood,abs(pcode));
+ //cout<<"!!!!! pcode, pmod, dedx ="<<pcode<<","<<mom<<","<<dedx<<endl;
+
+ }
+
+ pid->Tab();
+
+ Stat_t ng=hgood->GetEntries();
+ //cerr<<"Good tracks "<<ng<<endl;
+ Stat_t nf=hfound->GetEntries();
+ Stat_t nfak=hfake->GetEntries();
+ if (ng!=0)
+ cerr<<"Integral efficiency is about "<<nf/ng*100.<<" %\n";
+ //delete gAlice; //b.b.
+ cout<<"Done with nfound(good): "<<nf<<" + nfake: "<<nfak<<" = "<<nf+nfak<<" tracks."<<endl;
+// -------------------- b.b. -------------------
+
+ //return 0;
+}
+
+Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event) {
+ if (gAlice) {delete gAlice; gAlice=0;}
+
+ TFile *file=TFile::Open("galice.root");
+ if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);}
+ if (!(gAlice=(AliRun*)file->Get("gAlice"))) {
+ cerr<<"gAlice have not been found on galice.root !\n";
+ exit(5);
+ }
+
+ Int_t np=gAlice->GetEvent(event);
+
+ Int_t *good=new Int_t[np];
+ Int_t k;
+ for (k=0; k<np; k++) good[k]=0;
+
+ AliITS *ITS=(AliITS*)gAlice->GetDetector("ITS");
+ if (!ITS) {
+ cerr<<"can't get ITS !\n"; exit(8);
+ }
+ AliITSgeom *geom=ITS->GetITSgeom();
+ if (!geom) {
+ cerr<<"can't get ITS geometry !\n"; exit(9);
+ }
+
+ TFile *cf=TFile::Open("AliITSclustersV2.root");
+ if (!cf->IsOpen()){
+ cerr<<"Can't open AliITSclustersV2.root !\n"; exit(6);
+ }
+ char cname[100]; sprintf(cname,"TreeC_ITS_%d",event);
+ TTree *cTree=(TTree*)cf->Get(cname);
+ if (!cTree) {
+ cerr<<"Can't get cTree !\n"; exit(7);
+ }
+ TBranch *branch=cTree->GetBranch("Clusters");
+ if (!branch) {
+ cerr<<"Can't get clusters branch !\n"; exit(8);
+ }
+ TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+ branch->SetAddress(&clusters);
+
+ Int_t entr=(Int_t)cTree->GetEntries();
+ for (k=0; k<entr; k++) {
+ cTree->GetEvent(k);
+ Int_t ncl=clusters->GetEntriesFast(); if (ncl==0) continue;
+ Int_t lay,lad,det; geom->GetModuleId(k,lay,lad,det);
+ if (lay<1 || lay>6) {
+ cerr<<"wrong layer !\n"; exit(10);
+ }
+ while (ncl--) {
+ AliITSclusterV2 *pnt=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
+ Int_t l0=pnt->GetLabel(0);
+ Int_t l1=pnt->GetLabel(1);
+ Int_t l2=pnt->GetLabel(2);
+ Int_t mask=1<<(lay-1);
+ if (l0>=0) good[l0]|=mask;
+ if (l1>=0) good[l1]|=mask;
+ if (l2>=0) good[l2]|=mask;
+ }
+ }
+ clusters->Delete(); delete clusters;
+ delete cTree; //Thanks to Mariana Bondila
+ cf->Close();
+
+ ifstream in("good_tracks_tpc");
+ if (!in) {
+ cerr<<"can't get good_tracks_tpc !\n"; exit(11);
+ }
+ Int_t nt=0;
+ Double_t px,py,pz,x,y,z;
+ Int_t code,lab;
+ while (in>>lab>>code>>px>>py>>pz>>x>>y>>z) {
+ if (good[lab] != 0x3F) continue;
+ TParticle *p = (TParticle*)gAlice->Particle(lab);
+ gt[nt].lab=lab;
+ gt[nt].code=p->GetPdgCode();
+ //**** px py pz - in global coordinate system
+ gt[nt].px=p->Px(); gt[nt].py=p->Py(); gt[nt].pz=p->Pz();
+ gt[nt].x=gt[nt].y=gt[nt].z=0.;
+ nt++;
+ if (nt==max) {cerr<<"Too many good tracks !\n"; break;}
+ }
+
+ delete[] good;
+
+ delete gAlice; gAlice=0;
+ file->Close();
+
+ return nt;
+}
+
+
+//----------------------------------------
+void filltab(){
+ cout<<"Fill tab..";
+ int nev=0;
+ int track,pcode;
+ float signal,pmom;
+ for(int t=0;t<10000;t++){
+ track=t; signal=5.; pmom=1200.; pcode=321;
+ pid->SetEdep(10000*nev+track,signal);
+ pid->SetPmom(10000*nev+track,pmom);
+ pid->SetPcod(10000*nev+track,abs(pcode));
+ }
+ pid->Tab();
+ cout<<"Done."<<endl;
+}
+//
+void loadpid(){
+ if(pid==0){
+ TFile *f=new TFile("pidhit.root");
+ AliITSPid *pid=(AliITSPid*)f->Get("AliITSPid");
+ if(pid)
+ {cout<<"Load PID object from PIDHIT.ROOT file"<<endl;}
+ else{cout<<"ERROR load PID object from PIDHIT.ROOT file"<<endl;}
+ }
+ pid->Print(0);
+}
+
+void quit(){
+ gROOT->Reset();
+ gROOT->ProcessLine(".q");
+}
+//---------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+{
+ TFile *file = new TFile("galice.root");
+ file->ls();
+ delete gAlice;
+ gAlice = (AliRun*)(file->Get("gAlice"));
+
+ Int_t npart = gAlice->GetEvent(0);
+ assert(npart);
+ TObjArray* parray = gAlice->Particles();
+ cout<<" load_particles: N part="<<npart<<endl;
+
+}
+
+
+
+
+
+
+
+
+
+
+
--- /dev/null
+{
+ TObjArray tarray(2000);
+ TObjArray parray(2000);
+//{
+ TFile *cf=TFile::Open("AliITSclustersV2.root"); assert(cf);
+ AliITSgeom *geom=(AliITSgeom*)cf->Get("AliITSgeom"); assert(geom);
+ AliITStrackerV2 tracker(geom);
+ cf->Close();
+ TFile *tf=TFile::Open("AliITStracksV2.root");
+ if (!tf->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 3;}
+// TObjArray tarray(2000);
+ TTree *tracktree=(TTree*)tf->Get("TreeT_ITS_0;1");
+ if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
+ TBranch *tbranch=tracktree->GetBranch("tracks");
+ Int_t nentr=(Int_t)tracktree->GetEntries(),i;
+ for (i=0; i<nentr; i++) {
+ AliITStrackV2 *iotrack=new AliITStrackV2;
+ tbranch->SetAddress(&iotrack);
+ tracktree->GetEvent(i);
+
+ Int_t tpcLabel=iotrack->GetLabel();
+ // tracker.CookLabel(iotrack,0.);
+ Int_t itsLabel=iotrack->GetLabel();
+ if (itsLabel != tpcLabel) iotrack->SetLabel(-TMath::Abs(itsLabel));
+ if (tpcLabel < 0) iotrack->SetLabel(-TMath::Abs(itsLabel));
+ tarray.AddLast(iotrack);
+ }
+ cout<<" N rec. tracks="<<nentr<<endl;
+ tf->Close();
+//----------------------------------------
+TFile *fpid = new TFile("AliITStracksV2Pid.root","recreate");
+AliITStrackV2Pid pidtmp;
+AliITSPid* pid = new AliITSPid(1000);
+ TTree itsTree("ITSf","Tree with PID");
+ AliITStrackV2Pid *outpid=&pidtmp;
+ itsTree.Branch("pids","AliITStrackV2Pid",&outpid,32000,0);
+Float_t signal,pmom;
+Double_t xv,par[5];
+AliITStrackV2* track;
+Float_t lam,pt_1;
+Int_t pcode;
+for(Int_t ii=0;ii<nentr;ii++)
+{
+ track = (AliITStrackV2*)tarray[ii];
+ track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+ track->PropagateToVertex();
+ signal=track->GetdEdx();
+ track->GetExternalParameters(xv,par);
+ lam=TMath::ATan(par[3]);
+ pt_1=TMath::Abs(par[4]);
+ //cout<<"lam,pt_1="<<lam<<" " <<pt_1<<endl;
+ if(TMath::Abs(pt_1)>0)pmom=1.*(1./(pt_1*TMath::Cos(lam)));
+ pidtmp.fSignal=signal;
+ pcode=pid->GetPcode(signal,pmom);
+ pidtmp.fPcode=pcode;
+ pidtmp.fWpi=pidtmp.fWk=pidtmp.fWp=-1;
+ if(pcode==211)pidtmp.fWpi=1;
+ if(pcode==321)pidtmp.fWk=1;
+ if(pcode==2212)pidtmp.fWp=1;
+
+ pidtmp.fMom=pmom;
+ itsTree.Fill();
+}
+itsTree.Write();
+fpid->Close();
+cout<<"File AliITStracksV2Pid.root written"<<endl;
+//----------------------------------------
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root");
+fpid->ls();
+fpid->Close();
+//-------------------------------------------------
+ TFile *tfpid=TFile::Open("AliITStracksV2Pid.root");
+ if (!tfpid->IsOpen()) {cerr<<"Can't open AliITStracksV2Pid.root !\n"; return 3;}
+// TObjArray parray(2000);
+ TTree *pidtree=(TTree*)tfpid->Get("ITSf");
+ if (!pidtree) {cerr<<"Can't get a tree with ITS pid !\n"; return 4;}
+ TBranch *tbranch=pidtree->GetBranch("pids");
+ Int_t nentr=(Int_t)pidtree->GetEntries(),i;
+ for (i=0; i<nentr; i++) {
+ AliITStrackV2Pid *iopid=new AliITStrackV2Pid;
+ tbranch->SetAddress(&iopid);
+ pidtree->GetEvent(i);
+ parray.AddLast(iopid);
+ //cout<<" fWpi,k,p, Signal = "<<iopid->fWpi<<","<<iopid->fWk<<","<<iopid->fWp<<","<<iopid->fSignal <<endl;
+ //cout<<" Pcode,pmom="<<iopid->fPcode<<" "<<iopid->fMom<<endl;
+ }
+ cout<<" N rec. pids="<<nentr<<endl;
+ tfpid->Close();
+
+//-------------------------------------------------
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff