Part of new PID code from Boris Batyunya
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2002 09:56:44 +0000 (09:56 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Jan 2002 09:56:44 +0000 (09:56 +0000)
15 files changed:
ITS/AliITSPid.cxx [new file with mode: 0644]
ITS/AliITSPid.h [new file with mode: 0644]
ITS/AliITStrackV2.cxx
ITS/AliITStrackV2.h
ITS/AliITStrackV2Pid.h [new file with mode: 0644]
ITS/AliITStrackerV2.cxx
ITS/Config_pid.C [new file with mode: 0644]
ITS/ITSPIDtest.C [new file with mode: 0644]
ITS/ITSPidV2.sh [new file with mode: 0644]
ITS/Makefile
ITS/SetConvConst.C [new file with mode: 0644]
ITS/dEdXgeant.C [new file with mode: 0644]
ITS/dedxanal.C [new file with mode: 0644]
ITS/load_particles.C [new file with mode: 0644]
ITS/save_pid.C [new file with mode: 0644]

diff --git a/ITS/AliITSPid.cxx b/ITS/AliITSPid.cxx
new file mode 100644 (file)
index 0000000..884651d
--- /dev/null
@@ -0,0 +1,289 @@
+#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);
+}
+//-----------------------------------------------------------
+
diff --git a/ITS/AliITSPid.h b/ITS/AliITSPid.h
new file mode 100644 (file)
index 0000000..96f76d7
--- /dev/null
@@ -0,0 +1,58 @@
+#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 
+
+
+
+
index 3dfb99f..4c36090 100644 (file)
@@ -88,7 +88,13 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
   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];
+  }
+
 }
 /*
 //_____________________________________________________________________________
@@ -512,15 +518,15 @@ Int_t AliITStrackV2::Invariant() const {
   
   //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; 
@@ -537,7 +543,7 @@ Int_t AliITStrackV2::Invariant() const {
   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;
   }
   */
@@ -833,3 +839,54 @@ void AliITStrackV2::ResetCovariance() {
   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);
+}
+  
+
index aab0a25..012c961 100644 (file)
@@ -43,7 +43,9 @@ public:
   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); }
@@ -93,6 +95,8 @@ private:
 
   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
 };
 
@@ -104,6 +108,20 @@ void AliITStrackV2::GetExternalParameters(Double_t& xr, Double_t x[5]) const {
      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
 
diff --git a/ITS/AliITStrackV2Pid.h b/ITS/AliITStrackV2Pid.h
new file mode 100644 (file)
index 0000000..584928c
--- /dev/null
@@ -0,0 +1,23 @@
+#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
+
+
index 5564549..d161efd 100644 (file)
@@ -630,6 +630,10 @@ cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
      //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
diff --git a/ITS/Config_pid.C b/ITS/Config_pid.C
new file mode 100644 (file)
index 0000000..7c1fa74
--- /dev/null
@@ -0,0 +1,466 @@
+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");
+    }
+
+
+}
diff --git a/ITS/ITSPIDtest.C b/ITS/ITSPIDtest.C
new file mode 100644 (file)
index 0000000..46aebca
--- /dev/null
@@ -0,0 +1,69 @@
+{
+ 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
+
+}
+
+
+
+
+
+
+
+
diff --git a/ITS/ITSPidV2.sh b/ITS/ITSPidV2.sh
new file mode 100644 (file)
index 0000000..5a9b71d
--- /dev/null
@@ -0,0 +1,42 @@
+#!/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.
+
+
+
+
+
+
+
+
index 028ffeb..2eeefa3 100644 (file)
@@ -42,9 +42,9 @@ SRCS          = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \
                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
diff --git a/ITS/SetConvConst.C b/ITS/SetConvConst.C
new file mode 100644 (file)
index 0000000..feee514
--- /dev/null
@@ -0,0 +1,3 @@
+{
+AliKalmanTrack::SetConvConst(100/.299792458/.2);
+}
diff --git a/ITS/dEdXgeant.C b/ITS/dEdXgeant.C
new file mode 100644 (file)
index 0000000..95869db
--- /dev/null
@@ -0,0 +1,647 @@
+
+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++;                            
+}
+//=======================
+
+
+
+
+
+
diff --git a/ITS/dedxanal.C b/ITS/dedxanal.C
new file mode 100644 (file)
index 0000000..b612aaa
--- /dev/null
@@ -0,0 +1,801 @@
+//=========================================
+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");
+}
+//---------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ITS/load_particles.C b/ITS/load_particles.C
new file mode 100644 (file)
index 0000000..ea85bc4
--- /dev/null
@@ -0,0 +1,23 @@
+{
+   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;
+
+}
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ITS/save_pid.C b/ITS/save_pid.C
new file mode 100644 (file)
index 0000000..bfa48e9
--- /dev/null
@@ -0,0 +1,115 @@
+{
+   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();
+
+//-------------------------------------------------
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+