Updated version of the PID code (from B. Batyunya)
authorbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2002 06:47:35 +0000 (06:47 +0000)
committerbarbera <barbera@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 29 Jan 2002 06:47:35 +0000 (06:47 +0000)
ITS/AliBarrelReconstructionV2.C [new file with mode: 0644]
ITS/AliITSPid.cxx
ITS/AliITSPid.h
ITS/Barrel.sh [new file with mode: 0644]
ITS/ITSPidV2.sh
ITS/dedxanal.C
ITS/save_pidV2.C [new file with mode: 0644]
ITS/scan_pidV2.C [new file with mode: 0644]

diff --git a/ITS/AliBarrelReconstructionV2.C b/ITS/AliBarrelReconstructionV2.C
new file mode 100644 (file)
index 0000000..f15ba9b
--- /dev/null
@@ -0,0 +1,459 @@
+/****************************************************************************
+ * This macro is supposed to do reconstruction in the barrel ALICE trackers *
+ * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
+ * It does the following steps (April 12, 2001):                            *
+ *                   1) TPC cluster finding                                 *
+ *                   2) TPC track finding                                   *
+ *                   3) ITS cluster finding V2 (via fast points !)          *
+ *                   4) ITS track finding V2                                *
+ *                   5) ITS back track propagation V2                       *
+ *                   6) TPC back track propagation                          *
+ *                (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch)          *
+ ****************************************************************************/
+
+#ifndef __CINT__
+  #include "alles.h"
+  #include "AliMagF.h"
+  #include "AliTPCtracker.h"
+
+  #include "AliITS.h"
+  #include "AliITSgeom.h"
+  #include "AliITSRecPoint.h"
+  #include "AliITSclusterV2.h"
+  #include "AliITSsimulationFastPoints.h"
+  #include "AliITStrackerV2.h"
+  #include<iostream.h>
+#endif
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n);
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n);
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname,  Int_t n);
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
+
+Int_t ITSFindClusters(const Char_t *inname,  const Char_t *outname, Int_t n);
+Int_t ITSFindTracks(const Char_t *inname, const Char_t *inname2, const Char_t *outname, Int_t n);
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
+
+
+Int_t AliBarrelReconstructionV2(Int_t n=1) {
+  //const Char_t *TPCdigName="rfio:galice.root";
+   const Char_t *TPCdigName="galice.root";
+   const Char_t *TPCclsName="AliTPCclusters.root";
+   const Char_t *TPCtrkName="AliTPCtracks.root";
+   const Char_t *TPCtrkNameS="AliTPCtracksSorted.root";
+
+
+   //const Char_t *ITSdigName="rfio:galice.root";
+   const Char_t *ITSdigName="galice.root";
+   const Char_t *ITSclsName="AliITSclustersV2.root";
+   const Char_t *ITStrkName="AliITStracksV2.root";
+
+   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
+
+    
+// ********** Find TPC clusters *********** //
+   if (TPCFindClusters(TPCdigName,TPCclsName, n)) {
+      cerr<<"Failed to get TPC clusters !\n";
+      return 1;
+   }      
+
+// ********** Find TPC tracks *********** //
+    if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
+      cerr<<"Failed to get TPC tracks !\n";
+      return 1;
+    }
+  
+// ********** Sort and label TPC tracks *********** //
+   if (TPCSortTracks(TPCclsName,TPCtrkName,TPCtrkNameS,n)) {
+      cerr<<"Failed to sort TPC tracks !\n";
+      return 1;
+    }
+
+// ********** Find ITS clusters *********** //
+   if (ITSFindClusters(ITSdigName,ITSclsName,n)) {
+      cerr<<"Failed to get ITS clusters !\n";
+      return 1;
+   }
+
+// ********** Find ITS tracks *********** //
+   {
+     //TFile *clsFile=TFile::Open(ITSclsName);
+   if (ITSFindTracks(TPCtrkNameS,ITSclsName,ITStrkName,n)) {
+      cerr<<"Failed to get ITS tracks !\n";
+      return 1;
+   }
+   //clsFile->Close();
+   }
+   return 1;
+// ********** Back propagation of the ITS tracks *********** //
+   {TFile *clsFile=TFile::Open(ITSclsName);
+   if (ITSPropagateBack(ITStrkName,TPCtrkNameS)) {
+      cerr<<"Failed to propagate back the ITS tracks !\n";
+      return 1;
+   }
+   clsFile->Close();}
+
+
+// ********** Back propagation of the TPC tracks *********** //
+   {TFile *clsFile=TFile::Open(TPCclsName);
+   if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
+      cerr<<"Failed to propagate back the TPC tracks !\n";
+      return 1;
+   }
+   clsFile->Close();}
+
+   return 0;
+}
+
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
+   Int_t rc=0;
+   const Char_t *name="TPCFindClusters";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+
+   AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+   AliTPCv2 tpc;
+   tpc.SetParam(param);
+
+   //tpc.Digits2Clusters(out); //MI change
+   cout<<"TPCFindClusters: nev ="<<n<<endl;
+   for (Int_t i=0;i<n;i++){
+     printf("Processing event %d\n",i);
+     tpc.Digits2Clusters(out,i);
+     //         AliTPCclusterer::Digits2Clusters(dig, out, i);
+   }
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname, Int_t n) {
+   Int_t rc=0;
+   const Char_t *name="TPCFindTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+   AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+   //AliTPCtracker *tracker=new AliTPCtracker(param);
+   //rc=tracker->Clusters2Tracks(0,out);
+   //delete tracker;
+   cout<<"TPCFindTracks: nev ="<<n<<endl;
+
+   for (Int_t i=0;i<n;i++){
+     printf("Processing event %d\n",i);
+     AliTPCtracker *tracker = new AliTPCtracker(param,i);
+     //Int_t rc=
+       tracker->Clusters2Tracks(0,out);
+     delete tracker;
+   }
+
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *inname2, const Char_t * outname,  Int_t eventn){
+   Int_t rc=0;
+   const Char_t *name="TPCSortTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   TFile *out =TFile::Open(outname,"recreate");
+   TFile *in1 =TFile::Open(inname);
+   TFile *in2 =TFile::Open(inname2);
+
+
+   AliTPCParam *param=(AliTPCParam *)in1->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+   cout<<"TPCSortedtracks: nev ="<<eventn<<endl;
+
+
+   // loop over events 
+   for (Int_t event=0;event<eventn; event++){
+     TObjArray tarray(10000);
+     AliTPCtrack *iotrack=0;
+     Int_t i;
+
+
+     in2->cd();
+     char   tname[100];
+     sprintf(tname,"TreeT_TPC_%d",event);
+
+     TTree *tracktree=(TTree*)in2->Get(tname);
+     TBranch *tbranch=tracktree->GetBranch("tracks");
+     Int_t nentr=(Int_t)tracktree->GetEntries();
+     for (i=0; i<nentr; i++) {
+       iotrack=new AliTPCtrack;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
+       tarray.AddLast(iotrack);
+     }   
+     tarray.Sort();
+     // in2->Close();
+     
+     //assign thacks GEANT labels
+     in1->cd();
+     AliTPCtracker *tracker = new AliTPCtracker(param,event);
+     tracker->LoadInnerSectors();
+     tracker->LoadOuterSectors();
+     for (i=0; i<nentr; i++) {
+       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+       tracker->CookLabel(iotrack,0.1);
+     }   
+     delete tracker;
+     //in->Close();
+     //end of GEANT label assignment
+     
+     tracktree=new TTree(tname,"Tree with TPC tracks");
+     tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
+     for (i=0; i<nentr; i++) {
+       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+       tracktree->Fill();
+     }
+     out->cd();
+     tracktree->Write();
+   }
+
+   out->Close();
+   in2->Close();
+   in1->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
+   Int_t rc=0;
+   const Char_t *name="ITSFindClusters";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname,"update");
+
+   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+      cerr<<"Can't get gAlice !\n";
+      return 1;
+   }
+
+   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
+   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
+   AliITSgeom *geom=ITS->GetITSgeom();
+   out->cd();   
+   geom->Write();
+     
+   cout<<"ITSFindClusters: nev ="<<n<<endl;
+   Int_t ev=0;
+   for (ev = 0; ev<n; ev++){
+     in->cd();   // !!!! MI directory must point to galice. - othervise problem with Tree -connection
+     gAlice->GetEvent(ev);
+     //gAlice->TreeR()->Reset();   //reset reconstructed tree
+
+     
+     TTree *pTree=gAlice->TreeR();
+     if (!pTree){
+       gAlice->MakeTree("R");
+       pTree = gAlice->TreeR();
+     }
+     TBranch *branch=pTree->GetBranch("ITSRecPoints");
+     /*
+     if (!branch) {
+       //if not reconstructed ITS branch do reconstruction 
+       ITS->MakeBranch("R",0);
+
+       //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
+       AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
+       for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
+       Int_t nsignal=25;
+       Int_t size=-1;
+       Int_t bgr_ev=Int_t(ev/nsignal);
+       ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
+       //////////////////////////////////////////////////////////////////////////
+       gAlice->GetEvent(ev);   //MI comment  - in HitsToFast... they reset treeR to 0 
+       //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
+       //time
+     }
+     */
+
+     
+     out->cd();
+     TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+     char   cname[100];
+     sprintf(cname,"TreeC_ITS_%d",ev);
+  
+     TTree *cTree=new TTree(cname,"ITS clusters");
+     cTree->Branch("Clusters",&clusters);
+     
+     pTree=gAlice->TreeR();
+     if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
+     branch=pTree->GetBranch("ITSRecPoints");
+     if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
+     TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
+     branch->SetAddress(&points);
+     
+     TClonesArray &cl=*clusters;
+     Int_t nclusters=0;
+     Int_t nentr=(Int_t)pTree->GetEntries();
+     Int_t lab[6]; 
+     Float_t lp[5];
+     AliITSgeom *geom=ITS->GetITSgeom();
+
+     cout<<"!!!! nentr ="<<nentr<<endl;
+     for (Int_t i=0; i<nentr; i++) {
+       if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
+       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
+       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
+       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
+       Double_t yshift = x*rot[0] + y*rot[1];
+       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
+       Int_t ncl=points->GetEntriesFast();
+       nclusters+=ncl;
+
+     Float_t kmip=1.;
+     if(lay<5&&lay>2){kmip=280.;}; // b.b.
+     if(lay<7&&lay>4){kmip=38.;};
+
+
+       for (Int_t j=0; j<ncl; j++) {
+        AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
+        lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
+        lp[1]=p->GetZ()+zshift;
+        lp[2]=p->GetSigmaX2();
+        lp[3]=p->GetSigmaZ2();
+        lp[4]=p->GetQ(); lp[4]/=kmip;    // b.b.
+        if(ncl==11)cout<<"mod,cl,Q ="<<i<<","<<j<<","<<lp[4]<<endl;
+        lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
+        lab[3]=ndet;
+        
+        Int_t label=lab[0];
+        if(label>=0) { // b.b.
+        TParticle *part=(TParticle*)gAlice->Particle(label);
+        label=-3;
+        while (part->P() < 0.005) {
+          Int_t m=part->GetFirstMother();
+          if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
+          label=m;
+          part=(TParticle*)gAlice->Particle(label);
+        }
+        }
+        if      (lab[1]<0) lab[1]=label;
+        else if (lab[2]<0) lab[2]=label;
+        else cerr<<"No empty labels !\n";
+        
+        new(cl[j]) AliITSclusterV2(lab,lp);
+       }
+       cTree->Fill(); clusters->Delete();
+       points->Delete();
+     }
+     cTree->Write();
+     cerr<<"Number of clusters: "<<nclusters<<endl;
+     delete cTree; delete clusters; delete points;
+
+   }
+
+   delete gAlice; gAlice=0;
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSFindTracks(const Char_t * inname, const Char_t *inname2, const Char_t *outname, Int_t n) {
+   Int_t rc=0;
+   const Char_t *name="ITSFindTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+   TFile *in2 =TFile::Open(inname2);
+
+   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+
+   cout<<"ITSFindtracks: nev ="<<n<<endl;
+
+   for (Int_t i=0;i<n;i++){
+     AliITStrackerV2 tracker(geom,i);
+     rc=tracker.Clusters2Tracks(in,out);
+   }
+
+   in->Close();
+   in2->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
+   
+  Int_t rc=0;
+  /*
+   const Char_t *name="ITSPropagateBack";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+   AliITStrackerV2 tracker(geom);
+
+   TFile *out=TFile::Open(outname,"update");
+   TFile *in =TFile::Open(inname);
+   rc=tracker.PropagateBack(in,out);
+   in->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+  */
+   return rc;
+}
+
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="TPCPropagateBack";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+   AliTPCtracker *tracker=new AliTPCtracker(param);
+
+   TFile *out=TFile::Open(outname,"update");
+   TFile *in =TFile::Open(inname);
+   rc=tracker->PropagateBack(in,out);
+   delete tracker;
+   in->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+
+
+
+
index 99a4fe4..1494cc0 100644 (file)
@@ -1,5 +1,6 @@
 #include "AliITSPid.h"
 #include "TMath.h"
+#include "AliITSIOTrack.h"
 #include <iostream.h>
 ClassImp(AliITSPid)
 
@@ -12,7 +13,6 @@ Float_t AliITSPid::qcorr(Float_t xc)
   if(fcorr<=0.1)fcorr=0.1;
 return qtot/fcorr;
 }
-
 //-----------------------------------------------------------
 #include "vector.h"
 #include <algorithm>
@@ -34,7 +34,28 @@ Float_t AliITSPid::qtrm(Int_t track)
     for(Int_t i=0;i<nl;i++){q(i+1)=vf[i];} this->SetVec(track,q);
     return (q(6));
 }
+//........................................................
+Float_t AliITSPid::qtrm(Float_t qarr[6],Int_t narr)
+{
+  Float_t q[6],qm;
+  Int_t nl,ml;
+  narr>6?ml=6:ml=narr;
+  nl=0; for(Int_t i=0;i<ml;i++){if(qarr[i]>0.){q[i]=qarr[i];nl++;}}
+  if(nl==0)return 0.;
+    vector<float> vf(nl);
+    switch(nl){
+      case  1:qm=q[0]; break;
+      case  2:qm=(q[0]+q[1])/2.; break;
+      default:
+       for(Int_t i=0;i<nl;i++){vf[i]=q[i];}
+        sort(vf.begin(),vf.end());
+        qm= (vf[0]+vf[1])/2.;
+        break;
+    }
+    return qm;
+}
 //-----------------------------------------------------------
+//inline
 Int_t  AliITSPid::wpik(Int_t nc,Float_t q)
 {
        return pion();   
@@ -59,6 +80,7 @@ Int_t AliITSPid::wpik(Int_t nc,Float_t q)
     return pion();    
 }
 //-----------------------------------------------------------
+//inline
 Int_t  AliITSPid::wpikp(Int_t nc,Float_t q)
 {
        return pion();   
@@ -92,6 +114,58 @@ Int_t       AliITSPid::GetPcode(TClonesArray* rps,Float_t pm)
     return 0;    
 }
 //-----------------------------------------------------------
+Int_t   AliITSPid::GetPcode(AliTPCtrack *track)
+{
+      Double_t xk,par[5]; track->GetExternalParameters(xk,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]);
+      Float_t mom=1./(pt_1*TMath::Cos(lam));
+      Float_t dedx=track->GetdEdx();
+    Int_t pcode=GetPcode(dedx/40.,mom);
+    cout<<"TPCtrack dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+    return pcode?pcode:211;
+    }
+//------------------------------------------------------------
+/*
+Int_t   AliITSPid::GetPcode(AliITSIOTrack *track)
+{
+    Double_t px,py,pz;
+    px=track->GetPx();
+    py=track->GetPy();
+    pz=track->GetPz();
+    Float_t mom=TMath::Sqrt(px*px+py*py+pz*pz);
+//???????????????????
+   // Float_t dedx=1.0;
+     Float_t dedx=track->GetdEdx();
+//???????????????????    
+    Int_t pcode=GetPcode(dedx,mom);
+    cout<<"ITSV1 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+return pcode?pcode:211;
+}
+*/
+//-----------------------------------------------------------
+Int_t   AliITSPid::GetPcode(AliITStrackV2 *track)
+{
+  if(track==0)return 0;
+  //track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+      track->PropagateTo(3.,0.0028,65.19);
+
+      track->PropagateToVertex();
+    Double_t xk,par[5]; track->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.;};
+    Float_t dedx=track->GetdEdx();
+    //cout<<"lam,pt_1,mom,dedx="<<lam<<","<<pt_1<<","<<mom<<","<<dedx<<endl;
+    Int_t pcode=GetPcode(dedx,mom);
+    //cout<<"ITS V2 dedx,mom,pcode="<<dedx<<","<<mom<<","<<pcode<<endl;
+return pcode?pcode:211;
+}
+//-----------------------------------------------------------
 Int_t  AliITSPid::GetPcode(Float_t q,Float_t pm)
 {
     fWpi=fWk=fWp=-1.;     fPcode=0;
index 96f76d7..8de43be 100644 (file)
@@ -3,6 +3,9 @@
 #include <TObject.h>
 #include <TClonesArray.h>
 #include <TVector.h>
+#include "../TPC/AliTPCtrack.h"
+#include "AliITSIOTrack.h"
+#include "AliITStrackV2.h"
 #include <assert.h>
 //#include <stl.h>
 //___________________________________________________________________________
@@ -22,6 +25,9 @@ public:
        TVector* GetVec(Int_t track);
        Int_t   GetPcode(TClonesArray*,Float_t);
        Int_t   GetPcode(Float_t,Float_t);
+       Int_t   GetPcode(AliTPCtrack*);  // For PID TPC
+//        Int_t   GetPcode(AliITSIOTrack*); // For PID ITS tracking V1
+        Int_t   GetPcode(AliITStrackV2*); // For PID ITS tracking V2
        void    SetCut(Int_t,Float_t,Float_t,Float_t,
                            Float_t,Float_t,Float_t,Float_t);
        Float_t GetWpi(){return fWpi;}
@@ -42,6 +48,7 @@ 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);
        Int_t   wpik(Int_t,Float_t);
        Int_t   wpikp(Int_t,Float_t);
        Int_t   pion(){return fWpi=1.,fPcode=211;}
diff --git a/ITS/Barrel.sh b/ITS/Barrel.sh
new file mode 100644 (file)
index 0000000..ad4c696
--- /dev/null
@@ -0,0 +1,39 @@
+#!/bin/sh
+# This script processes the following steps for n events
+# (use the parameter n):
+# - the TPC and ITS slow simulation (hits+digits+clusters),
+# - the TPC+ITS tracking,
+# - the TPC and ITS PID
+# (Dubnna version)
+# 
+# delete eventual old files from the last run
+rm -f *.root
+#
+# run the hit generation
+aliroot -q -b "$ALICE_ROOT/macros/grun.C($1)" 
+#
+# digitize TPC
+aliroot -q -b "$ALICE_ROOT/TPC/AliTPCHits2Digits.C($1)"
+#
+# digitize ITS
+aliroot -q -b "$ALICE_ROOT/ITS/AliITSHits2DigitsDefault.C" 
+#
+# create reconstructed points for the ITS
+aliroot -q -b "$ALICE_ROOT/ITS/ITSDigitsToClusters.C(0,$[$1-1])" 
+#
+# do the TPC+ITS tracking
+aliroot -q -b "$ALICE_ROOT/ITS/AliBarrelReconstructionV2.C($1)" 
+#
+# 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/save_pidV2.C(0,$[$1-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/scan_pidV2.C(0,$[$1-1])"
+
+
index 5a9b71d..63b3df7 100644 (file)
@@ -1,42 +1,22 @@
 #!/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:
+# This script processes the ITS PID for n events
+# (use the parameter n):
+# - the AliTPCtest.C and AliITStest.C should be processed befor 
+# (Dubnna version)
+# 
+# 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.
-
-
-
-
-
+#
+rm -f AliITStrackingV2Pid.root
+#
+aliroot -q -b "$ALICE_ROOT/ITS/save_pidV2.C(0,$[$1-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/scan_pidV2.C(0,$[$1-1])"
 
 
 
index b612aaa..cd12e3c 100644 (file)
@@ -60,6 +60,7 @@ void dedxhis(Int_t pcod=0,Float_t pmin=0,Float_t pmax=1300)
 
       //qplotP->GetXaxis()->SetTitleSize(0.05);
       //qplotP->GetYaxis()->SetTitleSize(0.05);
+      //gPad->SetFillColor(kWhite);    // b.b.
       gStyle->SetOptStat(0);
       qplotP->SetXTitle("(Mev/c)");
       qplotP->SetYTitle("(mips)");
@@ -603,7 +604,9 @@ Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, const Int_t event);
     cerr<<"Track "<<lab<<" was not found !\n";
     continue;
     }
-  track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+    //track->Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+    track->PropagateTo(3.,0.0028,65.19);
+
   track->PropagateToVertex();
 
   if (lab==tlab) hfound->Fill(ptg);
diff --git a/ITS/save_pidV2.C b/ITS/save_pidV2.C
new file mode 100644 (file)
index 0000000..c930630
--- /dev/null
@@ -0,0 +1,79 @@
+////////////////////////////////////////////////////////////
+// Test macro for AliITSPid class and tracking version V2 // 
+////////////////////////////////////////////////////////////
+void
+save_pidV2(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);
+
+//
+//   Loop over events 
+//
+ for (int nev=0; nev<= evNumber2; nev++) {
+   char tname[30];
+   sprintf(tname,"TreeT_ITS_%d;1",nev);
+   TTree *tracktree=(TTree*)file->Get(tname);
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+       
+   Int_t nentr=tracktree->GetEntries();
+   cout<<"Found "<<nentr<<" ITS tracks in event No "<<nev<<endl;
+
+   char tpidname[30];
+   sprintf(tpidname,"TreeT%d",nev);
+   AliITStrackV2Pid pidtmp;
+   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);
+              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;
+       pidtmp.fSignal=signal;
+       pidtmp.fMom=mom;
+       itspidTree.Fill();
+       delete iotrack;
+   }// End for i (tracks)
+ fpid->Write();
+ }// End for nev (events)
+ fpid->Close();
+ file->Close();                 
+ cout<<"File AliITStracksV2Pid.root written"<<endl;
+ return;
+///////////////////////////////////////////////////////
+}
+
+
+
+
+
+
+
+
diff --git a/ITS/scan_pidV2.C b/ITS/scan_pidV2.C
new file mode 100644 (file)
index 0000000..5e0f5b4
--- /dev/null
@@ -0,0 +1,94 @@
+///////////////////////////////////////////////////////////
+// Test macro for AliITStracksV1Pid.root file            //
+// JINR Dubna Jan 2002                                   //
+///////////////////////////////////////////////////////////
+void
+scan_pidV2(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);
+     qplotKa.SetMarkerStyle(8); qplotKa.SetMarkerColor(kRed); qplotKa.SetMarkerSize(.3);
+     qplotPi.SetMarkerStyle(8); qplotPi.SetMarkerColor(kBlue); 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.);
+
+TFile *fpid = new TFile("AliITStracksV2Pid.root","read");
+fpid->ls();
+//
+//   Loop over events 
+//
+for (int nev=0; nev<= evNumber2; nev++) {
+  char tpidname[30];
+  sprintf(tpidname,"TreeT%d",nev);
+  TTree *tracktree=(TTree*)fpid->Get(tpidname);
+  TBranch *tbranch=tracktree->GetBranch("pids");
+       
+   Int_t nentr=tracktree->GetEntries();
+   cout<<"Found PID for "<<nentr<<" ITS V1 tracks on "<<tpidname<<endl;
+
+   AliITStrackV2Pid *iopid=0;
+for(Int_t ii=0;ii<nentr;ii++)
+  {
+      AliITStrackV2Pid *iopid=new AliITStrackV2Pid;
+      tbranch->SetAddress(&iopid);
+      tracktree->GetEvent(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;
+      delete iopid;
+  }// Enf for ii (tracks)
+ }// End for nev (events)
+ fpid->Close();
+  //...................... Draw histogramms .................
+   TCanvas *c1 = new TCanvas("PID_test","Scan PID ",200,10,900,700);
+   c1->Divide(2,1);
+  //.........................................................
+   c1->cd(1); gPad->SetFillColor(33);
+   signal_mip->Draw();
+
+   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; 
+  return;
+}
+