]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCcalibMaterial.cxx
moving AliEMCALGeoUtils to AliEMCALGeometry
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibMaterial.cxx
index 795e3461612f4f809ad6dc160b531a98d8fb8260..8ecc6421562e644c68e293856ef6db9191afc11f 100644 (file)
 #include "AliDCSSensorArray.h"
 #include "AliDCSSensor.h"
 
+#include "AliKFParticle.h"
+#include "AliKFVertex.h"
+
+#include "AliESDTZERO.h"
+
 ClassImp(AliTPCcalibMaterial)
 
 AliTPCcalibMaterial::AliTPCcalibMaterial():
@@ -119,7 +124,6 @@ Long64_t AliTPCcalibMaterial::Merge(TCollection *li) {
       return -1;
     }
     AliTPCcalibMaterial* calib= (AliTPCcalibMaterial*)(cal);
-    if(!calib) return 0;
     //
     if (!fHisMaterial) fHisMaterial=MakeHisto();
     fHisMaterial->Add(calib->fHisMaterial);
@@ -140,6 +144,8 @@ void AliTPCcalibMaterial::Process(AliESDEvent *event){
   const Float_t kMinDist=5;
   const Double_t kStep=1.;
   if (!event) return;
+  ProcessPairs(event);
+  return;
   //  TTreeSRedirector * cstream =  GetDebugStreamer();
   //
   if (!fHisMaterial){
@@ -150,7 +156,6 @@ void AliTPCcalibMaterial::Process(AliESDEvent *event){
   // fill histogram of track prolongations
   Float_t dca[2];
   Int_t ntracks = event->GetNumberOfTracks();
-
   for (Int_t itrack=0; itrack<ntracks; itrack++){
     AliESDtrack *track=event->GetTrack(itrack);
     if (!track) continue;
@@ -236,3 +241,470 @@ THnSparse *AliTPCcalibMaterial::MakeHisto(){
   fHisMaterialRPhi=sparseR;
   return sparse;
 }
+
+
+
+
+void  AliTPCcalibMaterial::ProcessPairs(AliESDEvent *event){
+  //
+  // Process pairs of tracks to get a material budget map
+  //
+  //
+  if (!event) return;
+  if (event) AliKFParticle::SetField(event->GetMagneticField());  // set mean magnetic field for KF particles
+  //
+  // 1. Calculate total dEdx for all TPC tracks
+  //
+  const Int_t kMinCl=70;
+  const Double_t kEpsilon=0.000001;
+  const Float_t kMinRatio=0.7;
+  const Float_t kMinDist=1.5;
+  const Float_t kMinDistChi2=8;        // 
+  const Float_t kMaxDistZ=280;         // max distanceZ
+  const Float_t kMaxDistR=250;          // max distanceR
+  const Double_t kMaxChi2   =36;     // maximal chi2 to define the vertex
+  const Double_t kMaxDistVertexSec=3;     // maximal distance to secondary vertex
+
+  if (!event) return;
+  AliESDVertex *vertexSPD =  (AliESDVertex *)event->GetPrimaryVertexSPD();
+  AliESDVertex * spdVertex   = (AliESDVertex *)event->GetPrimaryVertexSPD();
+  AliESDVertex * trackVertex = (AliESDVertex *)event->GetPrimaryVertexTracks();
+  AliESDVertex * tpcVertex   = (AliESDVertex *)event->GetPrimaryVertexTPC();
+  AliESDTZERO  * tzero       = (AliESDTZERO  *)event->GetESDTZERO() ;
+  //
+  Double_t tpcSignalTotPrim=0; 
+  Double_t tpcSignalTotSec=0; 
+  Int_t ntracksTPC=0;
+  Int_t nTPCPrim=0; 
+  Int_t nTPCSec=0;  
+  Int_t ntracks=event->GetNumberOfTracks();
+  if ( ntracks<=2 ) return;
+  
+  //
+  Float_t dca[2]={0};
+  Float_t cov[3]={0};
+  Float_t dca0[2]={0};
+  Float_t dca1[2]={0};
+  //
+  //1. Calculate total dEdx for primary and secondary tracks
+  //   and count primaries and secondaries
+  Int_t *rejectTrack = new Int_t[ntracks];
+  for (Int_t itrack=0; itrack<ntracks; itrack++){
+    AliESDtrack *track=event->GetTrack(itrack);
+    rejectTrack[itrack]=0;
+    if (!track) continue;
+    if (track->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
+    ntracksTPC++;
+    if ((1.+track->GetTPCNcls())/(1.+track->GetTPCNclsF())<=kMinRatio) continue;
+    if (!track->GetInnerParam()) continue;  // skip not TPC tracks
+    if (track->GetKinkIndex(0)!=0)  {rejectTrack[itrack]+=16;continue;} // skip kinks
+    track->GetImpactParameters(dca[0],dca[1]);
+    if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
+    // remove too dip secondaries
+    if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist){
+      tpcSignalTotPrim+=track->GetTPCsignal();
+      nTPCPrim++;
+    }else{
+      tpcSignalTotSec+=track->GetTPCsignal();
+      nTPCSec++;
+    };
+    if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Sqrt(dca[0]*dca[0]/(TMath::Abs(cov[0])))<kMinDistChi2) rejectTrack[itrack]+=1;  // primary
+    if (cov[0]>kEpsilon &&TMath::Abs(dca[0])>kEpsilon &&TMath::Abs(dca[0])<kMinDist)  rejectTrack[itrack]+=1;  // primary
+    if (track->GetTPCsignal()<40) rejectTrack[itrack]+=16;
+    //
+    if (CheckLooper(itrack, event))   rejectTrack[itrack]+=2;   // looper
+    if (CheckV0(itrack,event))       rejectTrack[itrack]+=4;
+  }
+  
+  
+  //
+  // 2. Find secondary vertices - double loop
+  //    
+  for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
+    AliESDtrack *track0=event->GetTrack(itrack0);
+    if (!track0) continue;
+    if (track0->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
+    if ((1.+track0->GetTPCNcls())/(1.+track0->GetTPCNclsF())<=kMinRatio) continue;
+    if (!track0->GetInnerParam()) continue;  // skip not TPC tracks
+    if (track0->GetKinkIndex(0)>0) continue; // skip kinks
+    if (rejectTrack[itrack0]) continue;   // skip
+    track0->GetImpactParameters(dca[0],dca[1]);
+    track0->GetImpactParameters(dca0[0],dca0[1]);
+    if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;
+    // remove too dip secondaries
+
+    AliKFParticle part0(*track0,211);  //assuming pion mass
+    if (track0->Charge()*part0.Q()<0) part0.Q()*=-1;  // change sign if opposite
+    //
+    for (Int_t itrack1=itrack0+1; itrack1<ntracks; itrack1++){
+      AliESDtrack *track1=event->GetTrack(itrack1);
+      if (!track1) continue;
+      if (rejectTrack[itrack1]) continue;   // skip
+      if (track1->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
+      if ((1.+track1->GetTPCNcls())/(1.+track1->GetTPCNclsF())<=kMinRatio) continue;
+      if (!track1->GetInnerParam()) continue;  // skip not TPC tracks
+      if (track1->GetKinkIndex(0)!=0) continue; // skip kinks
+      track1->GetImpactParameters(dca1[0],dca1[1]);
+      track1->GetImpactParameters(dca[0],dca[1]);
+      if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
+      if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;     
+      AliKFParticle part1(*track1,211); // assuming pion mass
+      if (track1->Charge()*part1.Q()<0) part1.Q()*=-1;  // change sign if opposite
+
+      //
+      //
+      AliKFVertex vertex;
+      vertex+=part0;
+      vertex+=part1;
+      if ((vertex.GetChi2()/vertex.GetNDF())> kMaxChi2) continue;
+      if (TMath::Abs(vertex.GetX())>kMaxDistR) continue;
+      if (TMath::Abs(vertex.GetY())>kMaxDistR) continue;
+      if (TMath::Abs(vertex.GetZ())>kMaxDistZ) continue;
+      Double_t errX2=vertex.GetErrX();
+      Double_t errY2=vertex.GetErrY();
+      Double_t errZ2=vertex.GetErrZ();
+      //
+      Double_t err3D=TMath::Sqrt(errX2*errX2+errY2*errY2+errZ2*errZ2/10.);  
+      Double_t err2D=TMath::Sqrt(errX2*errX2+errY2*errY2);  
+      if (err3D>kMaxDistVertexSec) continue;
+      if (err3D*TMath::Sqrt(vertex.GetChi2()+0.00001)>kMaxDistVertexSec) continue;
+
+      Double_t dvertex=0;
+      dvertex += (vertexSPD->GetX()-vertex.GetX())*(vertexSPD->GetX()-vertex.GetX());
+      dvertex += (vertexSPD->GetY()-vertex.GetY())*(vertexSPD->GetY()-vertex.GetY());
+      dvertex += (vertexSPD->GetZ()-vertex.GetZ())*(vertexSPD->GetZ()-vertex.GetZ());
+      dvertex=TMath::Sqrt(dvertex+0.00000001);
+      if (err3D>0.2*dvertex) continue;    
+      if (err3D*TMath::Sqrt(vertex.GetChi2()+0.000001)>0.1*dvertex) continue;
+      Double_t radius = TMath::Sqrt((vertex.GetX()*vertex.GetX()+vertex.GetY()*vertex.GetY()));
+      //
+      AliKFVertex vertex2;
+      vertex2+=part0;
+      vertex2+=part1;
+      //
+
+      for (Int_t itrack2=0; itrack2<ntracks; itrack2++){  
+       if (itrack2==itrack0) continue;
+       if (itrack2==itrack1) continue;
+       if (rejectTrack[itrack2]) continue;   // skip
+       AliESDtrack *track2=event->GetTrack(itrack2);
+       if (!track2) continue;
+       if (track2->GetTPCNcls()<=kMinCl) continue;  // skip short tracks
+       if ((1.+track2->GetTPCNcls())/(1.+track2->GetTPCNclsF())<=kMinRatio) continue;
+       if (!track2->GetInnerParam()) continue;  // skip not TPC tracks
+       if (track2->GetKinkIndex(0)>0) continue; // skip kinks
+       track2->GetImpactParameters(dca[0],dca[1]);
+       if (TMath::Abs(dca[0])<kMinDist && TMath::Abs(dca[1])<kMinDist) continue;
+       if (TMath::Abs(dca[0])>kMaxDistR && TMath::Abs(dca[1])>kMaxDistZ) continue;     
+       if (TMath::Abs(track2->GetD(vertex.GetX(), vertex.GetY(),event->GetMagneticField()))>kMaxDistVertexSec) continue;
+       Double_t vtxx[3]={vertex2.GetX(),vertex2.GetY(),vertex2.GetZ()};
+       Double_t svtxx[3]={vertex.GetErrX(),vertex.GetErrY(),vertex.GetErrZ()};
+       AliESDVertex vtx(vtxx,svtxx);
+       AliExternalTrackParam param=*track2;
+       Double_t delta[2]={0,0};
+       if (!param.PropagateToDCA(&vtx,event->GetMagneticField(),kMaxDistVertexSec,delta)) continue;    
+       if (TMath::Abs(delta[0])>kMaxDistVertexSec) continue;
+       if (TMath::Abs(delta[1])>kMaxDistVertexSec) continue;      
+       if (TMath::Abs(delta[0])>6.*TMath::Sqrt(param.GetSigmaY2()+vertex2.GetErrY()*vertex2.GetErrY())+0.1) continue; 
+       if (TMath::Abs(delta[1])>6.*TMath::Sqrt(param.GetSigmaZ2()+vertex2.GetErrZ()*vertex2.GetErrZ())+0.5) continue; 
+       //
+       AliKFParticle part2(param,211); // assuming pion mass
+       if (track2->Charge()*part2.Q()<0) part2.Q()*=-1;  // change sign if opposite
+       vertex2+=part2;
+       rejectTrack[itrack0]+=10;  // do noit reuse the track
+       rejectTrack[itrack1]+=10;  // do not reuse the track      
+       rejectTrack[itrack2]+=10;
+      }
+      
+
+      TTreeSRedirector *pcstream = GetDebugStreamer();      
+      if (pcstream){
+       //
+       //
+       //
+       Float_t dedx0= track0->GetTPCsignal(); 
+       Float_t dedx1= track1->GetTPCsignal();
+       AliExternalTrackParam * p0= (AliExternalTrackParam *)track0->GetInnerParam();
+       AliExternalTrackParam * p1= (AliExternalTrackParam *)track1->GetInnerParam();
+       Double_t errX=vertex2.GetErrX();
+       Double_t errY=vertex2.GetErrY();
+       Double_t errZ=vertex2.GetErrZ();
+       Double_t vx = vertex2.GetX();
+       Double_t vy = vertex2.GetY();
+       Double_t vz = vertex2.GetZ();
+       (*pcstream)<<"mapTPC"<<
+         "run="<<fRun<<                     // run
+         "event="<<fEvent<<          //  event number
+         "time="<<fTime<<         // timeStamp
+         "trigger="<<fTrigger<<      //  trigger
+         "mag="<<fMagF<<             //  magnetic field
+         //
+         "spdV.="<<spdVertex<<            // spd vertex
+         "trackV.="<<trackVertex<<        // track vertex
+         "tpcV.="<<tpcVertex<<            // track vertex
+         "tzero.="<<tzero<<               // tzero info
+         //
+         "ntracks="<<ntracks<<
+         "ntracksTPC="<<ntracksTPC<<
+         "nPrim="<<nTPCPrim<<              // number of primaries
+         "nSec="<<nTPCSec<<                // number of secondaries
+         "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
+         "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
+         "dedx0="<<dedx0<<                 // dedx part 0
+         "dedx1="<<dedx1<<                 // dedx part 1
+         "p0.="<<p0<<                      // part 0
+         "p1.="<<p1<<                       //part 1
+         "v.="<<&vertex<<                  // KF vertex
+         "v2.="<<&vertex2<<                // KF vertex all tracks
+         "z0="<<dca0[1]<<
+         "z1="<<dca1[1]<<
+         "rphi0="<<dca0[0]<<
+         "rphi1="<<dca1[0]<<
+         "dvertex="<<dvertex<<
+         "radius="<<radius<<
+         "vx="<<vx<<
+         "vy="<<vy<<
+         "vz="<<vz<<
+         "errX="<<errX<<
+         "errY="<<errY<<
+         "errZ="<<errZ<<
+         "err2D="<<err2D<<
+         "err3D="<<err3D<<       
+         "\n";
+       //
+       if (vertex2.GetNDF()>2){
+         (*pcstream)<<"mapVertex"<<
+           "run="<<fRun<<                     // run
+           "event="<<fEvent<<          //  event number
+           "time="<<fTime<<         // timeStamp
+           "trigger="<<fTrigger<<      //  trigger
+           "mag="<<fMagF<<             //  magnetic field
+           //
+           "spdV.="<<spdVertex<<            // spd vertex
+           "trackV.="<<trackVertex<<        // track vertex
+           "tpcV.="<<tpcVertex<<            // track vertex
+           "tzero.="<<tzero<<               // tzero info
+           //
+           //
+           "ntracks="<<ntracks<<
+           "ntracksTPC="<<ntracksTPC<<
+           "nPrim="<<nTPCPrim<<              // number of primaries
+           "nSec="<<nTPCSec<<                // number of secondaries
+           "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
+           "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
+           "dedx0="<<dedx0<<                 // dedx part 0
+           "dedx1="<<dedx1<<                 // dedx part 1
+           "p0.="<<p0<<                      // part 0
+           "p1.="<<p1<<                       //part 1
+           "v.="<<&vertex<<                  // KF vertex
+           "v2.="<<&vertex2<<                // KF vertex all tracks
+           "z0="<<dca0[1]<<
+           "z1="<<dca1[1]<<
+           "rphi0="<<dca0[0]<<
+           "rphi1="<<dca1[0]<<
+           "radius="<<radius<<
+           "vx="<<vx<<
+           "vy="<<vy<<
+           "vz="<<vz<<
+           "errX="<<errX<<
+           "errY="<<errY<<
+           "errZ="<<errZ<<
+           "err2D="<<err2D<<
+           "err3D="<<err3D<<     
+           "dvertex="<<dvertex<<
+           "\n";         
+       }
+      }
+    }
+  }
+  TTreeSRedirector *pcstream = GetDebugStreamer();
+  if (pcstream){    
+    (*pcstream)<<"statTPC"<<
+      "run="<<fRun<<                     // run
+      "time="<<fTime<<         // timeStamp
+      "trigger="<<fTrigger<<      //  trigger
+      "mag="<<fMagF<<             //  magnetic field
+      "ntracks="<<ntracks<<
+      "ntracksTPC="<<ntracksTPC<<
+      //
+      "nPrim="<<nTPCPrim<<              // number of primaries
+      "nSec="<<nTPCSec<<                // number of secondaries
+      "sigPrim="<<tpcSignalTotPrim<<    // total dEdx in primaries
+      "sigSec="<<tpcSignalTotSec<<      // total dEdx in secondaries
+      //
+      "spdV.="<<spdVertex<<            // spd vertex
+      "trackV.="<<trackVertex<<        // track vertex
+      "tpcV.="<<tpcVertex<<            // track vertex
+      "tzero.="<<tzero<<               // tzero info
+      "\n";
+  }
+  delete [] rejectTrack;
+}
+
+
+Bool_t AliTPCcalibMaterial::CheckLooper(Int_t index, AliESDEvent *event){
+  //
+  // check if given track is looper candidate
+  // if looper return kTRUE
+  // 
+  Int_t ntracks=event->GetNumberOfTracks();
+  Int_t index1=-1;
+  const Double_t ktglCut=0.03;
+  const Double_t kalphaCut=0.4;
+  static Int_t counter=0;
+  //
+  AliESDtrack * track0 = event->GetTrack(index);
+  AliESDtrack * track1P = 0;
+  for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
+    if (itrack1==index) continue;
+    AliESDtrack *track1=event->GetTrack(itrack1);
+    if (!track1) continue;
+    if (TMath::Abs(TMath::Abs(track1->GetTgl())-TMath::Abs(track0->GetTgl()))>ktglCut) continue;
+    if (TMath::Abs(TMath::Abs(track1->GetAlpha())-TMath::Abs(track0->GetAlpha()))>kalphaCut) continue;
+    index1=index;
+    track1P=track1;
+  }
+  if (index1>=0){
+    TTreeSRedirector *pcstream = GetDebugStreamer();      
+    if (pcstream &&counter<100000){
+      counter++;
+      AliExternalTrackParam p0(*track0);
+      AliExternalTrackParam p1(*track1P);
+      (*pcstream)<<"checkLooper"<<
+       "p0.="<<&p0<<
+       "p1.="<<&p1<<
+       "\n";
+    }
+    return kTRUE;
+  }
+  return kFALSE;
+}
+
+Bool_t AliTPCcalibMaterial::CheckV0(Int_t index, AliESDEvent *event){
+  //
+  // check if given track is V0 candidata
+  // if looper return kTRUE
+  // 
+  return kFALSE;
+  Int_t ntracks=event->GetNumberOfTracks();
+  Int_t index1=-1;
+  const Double_t kSigmaMass=0.001;
+  const Int_t kChi2Cut=10;
+  static Int_t counter=0;
+  //
+  AliESDtrack * track0 = event->GetTrack(index);
+  AliExternalTrackParam pL(*track0);
+  AliKFParticle part0El(*track0, 11);  //assuming  mass e
+  AliKFParticle part0Pi(*track0, 211);  //assuming  mass pi0
+  AliKFParticle part0P(*track0, 2212);  //assuming  mass proton
+  if (track0->Charge()*part0El.Q()<0) {
+    part0El.Q()*=-1;  // change sign if opposite
+    part0Pi.Q()*=-1;  // change sign if opposite
+    part0P.Q()*=-1;   // change sign if opposite
+  }
+  Bool_t isGamma=0;
+  Bool_t isK0=0;
+  
+  for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
+    if (itrack1==index) continue;
+    AliESDtrack *track1=event->GetTrack(itrack1);
+    if (!track1) continue;
+    if (track1->Charge()*track0->Charge()>0) continue;
+    AliKFParticle part1El(*track1, 11);  //assuming  mass e
+    AliKFParticle part1Pi(*track1, 211);  //assuming  mass e
+    AliKFParticle part1P(*track1, 2212);  //assuming  mass e
+    if (track1->Charge()*part1El.Q()<0) {
+      part1El.Q()*=-1;  // change sign if opposite
+      part1Pi.Q()*=-1;  // change sign if opposite
+      part1P.Q()*=-1;   // change sign if opposite
+    }
+    //
+    AliKFVertex vertexG;  // gamma conversion candidate
+    vertexG+=part0El;
+    vertexG+=part1El;
+    AliKFVertex vertexGC;  // gamma conversion candidate
+    vertexGC+=part0El;
+    vertexGC+=part1El;
+    vertexGC.SetMassConstraint(0,kSigmaMass);
+    AliKFVertex vertexK0;  // gamma conversion candidate
+    vertexK0+=part0Pi;
+    vertexK0+=part1Pi;
+    AliKFVertex vertexK0C;  // gamma conversion candidate
+    vertexK0C+=part0Pi;
+    vertexK0C+=part1Pi;
+    vertexK0C.SetMassConstraint(0.497614,kSigmaMass);
+    if (vertexGC.GetChi2()<kChi2Cut && vertexG.GetMass()<0.06)      isGamma=kTRUE;
+    if (vertexK0C.GetChi2()<kChi2Cut&&TMath::Abs(vertexK0.GetMass()-0.5)<0.06)  isK0=kTRUE;
+    if (isGamma||isK0) {
+      index1=index;
+      TTreeSRedirector *pcstream = GetDebugStreamer();      
+      if (pcstream&&counter<2000){
+       counter++;
+       AliExternalTrackParam p0(*track0);
+       AliExternalTrackParam p1(*track1);
+       (*pcstream)<<"checkV0"<<
+         "p0.="<<&p0<<  //particle 0
+         "p1.="<<&p1<<  //particle 1
+         "isGamma="<<isGamma<<  // is gamma candidate
+         "isK0s="<<isK0<<      // is k0 candidate
+         "vG.="<<&vertexG<<     // is gamma candidate
+         "vGC.="<<&vertexGC<<   // is gamma candidate
+         "vK.="<<&vertexK0<<     // is K0s candidate
+         "vKC.="<<&vertexK0C<<   // is K0s candidate
+         "\n";
+      }
+      break;
+    }
+  }
+  if (index1>0) return kTRUE;
+  return kFALSE;
+}
+
+/*
+  //AliXRDPROOFtoolkit::FilterList("mater.list","* mapTPC",1);
+  AliXRDPROOFtoolkit toolkit;
+  TChain *chain   = toolkit.MakeChainRandom("mater.list.Good","mapVertex",0,4000);
+  TChain *chainTPC   = toolkit.MakeChainRandom("mater.list.Good","mapTPC",0,50000);
+  
+  TCut cutErr="sqrt(v.fChi2)*err3D<1.0";
+  TCut cutOccu="sqrt(v.fChi2*(errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.2&&sqrt((errX^2+errY^2))/min(sqrt(v.fP[0]^2+v.fP[1]^2+v.fP[2]^2/20),20)<0.3&&sqrt(v.fChi2)*err3D<1.5";
+  //
+  chainTPC->Draw(">>listTPC",cutOccu,"entryList");
+  TEntryList *elistTPC = (TEntryList*)gDirectory->Get("listTPC");
+  chainTPC->SetEntryList(elistTPC);
+
+  //
+  //
+  
+  chainTPC.Draw("v.fP[1]:v.fP[0]","sqrt(v.fP[0]^2+v.fP[1]^2)<100","",10000000);
+
+  TCut cutITS="abs(v.fP[2])-10<sqrt(v.fP[0]^2+v.fP[1]^2)";
+  chainTPC.Draw("v.fP[1]:v.fP[0]",cutITS+"(sqrt(v.fP[0]^2+v.fP[1]^2)<50)&&err2D<1&&err3D/dvertex<0.05","",500000);
+
+
+  chainTPC.Draw("v.fP[1]:v.fP[0]>>his(400,-100,100,400,-100,100)","err2D<1&&err2D/dvertex<0.2","colz",50000); 
+
+
+  chainTPC.Draw("radius:vz>>his(400,-100,100,400,0,100)","err2D<1&&err2D/dvertex<0.05","colz",5000); 
+
+  TTree * tree = chain->CloneTree();
+  TFile f("material.root","recreate");
+  tree->Write();
+  f.Close();
+  
+
+
+  TCut cutChi2="v.fChi2<2&&sqrt(errX^2+errY^2)<2."; 
+  TFile f("material.root")
+  TTree * tree = (TTree*)f.Get("mapVertex");
+
+  tree->Draw("v.fChi2",cutChi2,"",1000)
+
+  
+
+
+*/
+