New PID classes and macros for Dubna group - PID weights tuned on HIJING
authormasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 May 2003 13:47:13 +0000 (13:47 +0000)
committermasera <masera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 16 May 2003 13:47:13 +0000 (13:47 +0000)
ITS/AliITSPid.cxx
ITS/AliITSPid.h
ITS/AliITSSavePIDV2.C
ITS/AliITSScanPIDV2.C
ITS/AliITStestPIDV2.sh
ITS/AliITStrackV2Pid.cxx
ITS/AliITStrackV2Pid.h

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