ITS/AliITSSecondaryVertices.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 //-- --- standard headers-------------
   3 #include <Riostream.h>
   4 //--------Root headers ---------------
   5 #include <TSystem.h>
   6 #include <TFile.h>
   7 #include <TStopwatch.h>
   8 #include <TObject.h>
   9 #include <TTree.h>
  10 #include <TChain.h>
  11 #include <TF1.h>
  12 #include <TLegend.h>
  13 #include <TLegendEntry.h>
  14 #include <TH1F.h>
  15 #include <TH2F.h>
  16 #include <TCut.h>
  17 #include <TParticle.h>
  18 #include <TCanvas.h>
  19 #include <TStyle.h>
  20 #include <TGraph.h>
  21 //----- AliRoot headers ---------------
  22 //#include "alles.h"
  23 #include "AliRun.h"
  24 #include "AliMagF.h"
  25 #include "AliKalmanTrack.h"
  26 #include "AliESDVertex.h"
  27 #include "AliITSVertexer.h"
  28 #include "AliITSVertexerTracks.h"
  29 #include <AliHeader.h>
  30 #include <AliGenEventHeader.h>
  31 #include "AliRun.h"
  32 #include "AliRunLoader.h"
  33 #include "AliITStrackV2.h"
  34 #include "AliITSSimpleVertex.h"
  35 #include "AliITSStrLine.h"
  36 #include "AliITSLoader.h"
  37 #include <AliESD.h>
  38 #include <AliStack.h>
  39 //-------------------------------------
  40 #endif
  41
  42
  43 void AliITSSecondaryVertices(Int_t ntracks, Int_t *pos) {
  44
  45   // example of usige AliITSVertexerTracks to get a (secodnary) vertex from a set of tracks.
  46
  47   Double_t field = 0.4;
  48
  49   AliRunLoader *rl = AliRunLoader::Open("galice.root");
  50   if (!rl) {
  51     cerr<<"Can't start session !\n";
  52   }
  53
  54   rl->LoadgAlice();
  55   if (rl->GetAliRun()){
  56     AliMagF * magf = gAlice->Field();
  57     AliKalmanTrack:: SetFieldMap(magf);
  58     AliKalmanTrack::SetUniformFieldTracking();
  59   }else {
  60     cerr<<"Can't get AliRun !\n";
  61
  62   }
  63
  64   field=rl->GetAliRun()->Field()->SolenoidField()/10.;
  65   rl->LoadHeader();
  66
  67
  68   TFile *ef=TFile::Open("AliESDs.root");
  69   if (!ef || !ef->IsOpen()) {cerr<<"Can't AliESDs.root !n";}
  70   AliESD* event = new AliESD;
  71   TTree* tree = (TTree*) ef->Get("esdTree");
  72   if (!tree) {cerr<<"no ESD tree foundn";};
  73   tree->SetBranchAddress("ESD", &event);
  74
  75
  76   Int_t e=0;
  77   tree->GetEvent(e);
  78   Int_t ntrk=event->GetNumberOfTracks();
  79   cout<<"Number of ESD tracks : "<<ntrk<<endl;
  80
  81   // Open input and output files
  82   TFile *outFile = TFile::Open("AliITSVertices.root","recreate");
  83
  84
  85   TTree *trkTree = new TTree("TreeT_ITS","its tracks");
  86   AliITStrackV2 *itstrack = 0;
  87   trkTree->Branch("tracks","AliITStrackV2",&itstrack,ntrk,0);
  88   Int_t pipe=0;
  89   for(Int_t i=0; i<ntrk; i++) {
  90     AliESDtrack *esdTrack = (AliESDtrack*)event->GetTrack(i);
  91     UInt_t status=esdTrack->GetStatus();
  92     UInt_t flags=AliESDtrack::kTPCin|AliESDtrack::kITSin;
  93     if ((status&AliESDtrack::kTPCin)==0)continue;
  94     if ((status&AliESDtrack::kITSin)==0)continue;
  95     if ((status&AliESDtrack::kITSrefit)==0) continue;
  96     if ((status&flags)==status) pipe=1;
  97
  98     itstrack = new AliITStrackV2(*esdTrack);
  99     if (pipe!=0) {
 100       itstrack->PropagateTo(3.,0.0028,65.19);
 101       itstrack->PropagateToVertex();  // This is not "exactly" the vertex
 102     }
 103     Int_t nclus=itstrack->GetNumberOfClusters();
 104
 105     if(nclus<6)continue;
 106     trkTree->Fill();
 107     itstrack = 0;
 108     // delete esdTrack;
 109   }
 110   // Primary vertex
 111   Double_t vtx[3];
 112   event->GetVertex()->GetXYZ(vtx);
 113   cout<<"Primary Vertex:"<<endl;
 114   event->GetVertex()->PrintStatus();
 115   cout<<"\n\n\n"<<endl;
 116
 117
 118   // Create vertexer
 119   AliITSVertexerTracks *vertexer = new AliITSVertexerTracks();
 120   vertexer->SetDCAcut(0);
 121   vertexer->SetDebug(0);
 122
 123
 124   Double_t vertF1[3],vertF2[3],vertF3[3],vertF4[3],vertF5[3];
 125   Int_t ncombi1,ncombi2,ncombi3,ncombi4,ncombi5;
 126   Double_t sig1,sig2,sig3,sig4,sig5;
 127   AliITSSimpleVertex *vert=0;
 128
 129
 130   // Calculate vertex from tracks in pos
 131   // BEWARE: the method VertexForSelectedTracks returns the address of the data member fSimpVert
 132   // which is always the same for all calls to VertexForSelectedTracks.
 133
 134   cout<<"\nStraight Line Min Dist finder with weights"<<endl;
 135   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,1);
 136   vert->Print();
 137   vert->GetXYZ(vertF1);
 138   ncombi1=vert->GetNContributors();
 139   sig1=vert->GetDispersion();
 140
 141
 142   cout<<"Straight Line Min Dist finder"<<endl;
 143   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,2);
 144   vert->Print();
 145   vert->GetXYZ(vertF2);
 146   ncombi2=vert->GetNContributors();
 147   sig2=vert->GetDispersion();
 148
 149   cout<<"Helix finder"<<endl;
 150   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,3);
 151   vert->Print();
 152   vert->GetXYZ(vertF3);
 153   ncombi3=vert->GetNContributors();
 154   sig3=vert->GetDispersion();
 155
 156   cout<<"Straight Line finder with weights"<<endl;
 157   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,4);
 158   vert->Print();
 159   vert->GetXYZ(vertF4);
 160   ncombi4=vert->GetNContributors();
 161   sig4=vert->GetDispersion();
 162
 163   cout<<"Straight Line finder with weights"<<endl;
 164   vert=(AliITSSimpleVertex*) vertexer->VertexForSelectedTracks(event,ntracks,pos,5);
 165   vert->Print();
 166   vert->GetXYZ(vertF5);
 167   ncombi5=vert->GetNContributors();
 168   sig5=vert->GetDispersion();
 169
 170
 171   delete vertexer;
 172
 173   outFile->Close();
 174   return;
 175 }