Example how to read reconstructed tracks, and load generated particles in pp-events
[u/mrichter/AliRoot.git] / HLT / exa / trigger_pp.C
CommitLineData
b7556f61 1void trigger()
2{
3
4 //Get the tracks:
5 AliL3TrackArray *tracks = new AliL3TrackArray();
6 AliL3FileHandler *file = new AliL3FileHandler();
7 file->SetBinaryInput("tracks.raw");
8 file->Binary2TrackArray(tracks);
9 file->CloseBinaryInput();
10 delete file;
11
12 Int_t ntracks=0;
13 Double_t xc,yc,zc;
14 Double_t impact;
15 AliL3Vertex vertex;
16 AliL3Fitter *fitter = new AliL3Fitter(&vertex);
17 fitter->LoadClusters("./");
18 //fitter->NoVertex();
19 for(int i=0; i<tracks->GetNTracks(); i++)
20 {
21 track = (AliL3Track*)tracks->GetCheckedTrack(i);
22 if(!track) continue;
23 if(fabs(track->GetPseudoRapidity())>0.9) continue;
24 if(track->GetNHits() < 50) continue;
25 if(track->GetPt()<0.1) continue;
26 track->CalculateHelix();
27 //cout<<"Pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<" Nhits "<<track->GetNHits()<<endl;
28 //cout<<"Before second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
29 fitter->FitHelix(track);
30 track->CalculateHelix();
31
32 track->GetClosestPoint(&vertex,xc,yc,zc);
33 impact = sqrt(xc*xc+yc*yc+zc*zc);
34 if(impact>3) continue;
35 cout<<"Number of hits "<<track->GetNHits()<<endl;
36 cout<<"Transversal impact "<<sqrt(xc*xc+yc*yc)<<endl;
37 cout<<"Longitudinal impact "<<zc<<endl;
38 cout<<"Total impact "<<sqrt(xc*xc+yc*yc+zc*zc)<<endl;
39 cout<<"xc "<<xc<<" yc "<<yc<<" zc "<<zc<<endl;
40
41 //cout<<"After second fit; pt "<<track->GetPt()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
42 cout<<"pt "<<track->GetPt()<<" eta "<<track->GetPseudoRapidity()<<endl;
43
44 cout<<"-------------------------------------"<<endl;
45 ntracks++;
46 }
47 cout<<endl<<"There was "<<ntracks<<" found tracks"<<endl;
48 delete tracks;
49 delete fitter;
50}
51
52enum tagprimary {kPrimaryCharged=0x4000};
53void LoadEvent(Int_t event=0)
54{
55 //Load the generated particles
56
57 gROOT->LoadMacro("$(ALICE_ROOT)/macros/loadlibs.C");
58 loadlibs();
59 TFile *rootfile = TFile::Open("/prog/alice/data/pro/25event_pp.root","READ");
60 gAlice = (AliRun*)rootfile->Get("gAlice");
61 Int_t nparticles = gAlice->GetEvent(event);
62 Int_t nprim = FindPrimaries(nparticles,0.,0.,0.);
63 cout<<"Number of primaries "<<nprim<<endl;
64 int co=0;
65 for(Int_t i=0; i<nparticles; i++)
66 {
67 TParticle *part = gAlice->Particle(i);
68 if(!part->TestBit(kPrimaryCharged)) continue;
69 if(fabs(part->Eta())>0.9) continue;
70 cout<<part->GetName()<<" pt "<<part->Pt()<<" eta "<<part->Eta()<<" xvert "<<part->Vx()<<" yvert "<<part->Vy()<<" zvert "<<part->Vz()<<endl;
71 co++;
72 }
73 cout<<endl<<"Number of primary tracks in the detector: "<<co<<endl;
74 gAlice=0;
75 rootfile->Close();
76
77}
78
79Int_t FindPrimaries(Int_t nparticles, Double_t xori, Double_t yori, Double_t zori)
80{
81 //Define primary particles in a pp-event. Code taken from offline.
82
83
84 // cuts:
85 Double_t vertcut = 0.001; // cut on the vertex position
86 Double_t decacut = 3.; // cut if the part. decays close to the vert.
87 Double_t timecut = 0.;
88
89 TList *listprim = new TList();
90 listprim->SetOwner(kFALSE);
91
92 Int_t nprch1=0;
93 Int_t nprch2=0;
94 for(Int_t iprim = 0; iprim<nparticles; iprim++){ //loop on tracks
95
96 TParticle * part = gAlice->Particle(iprim);
97 char *xxx=strstr(part->GetName(),"XXX");
98 if(xxx)continue;
99
100 TParticlePDG *ppdg = part->GetPDG();
101 if(TMath::Abs(ppdg->Charge())!=3)continue; // only charged (no quarks)
102
103 Double_t dist=TMath::Sqrt((part->Vx()-xori)*(part->Vx()-xori)+(part->Vy()-yori)*(part->Vy()-yori)+(part->Vz()-zori)*(part->Vz()-zori));
104 if(dist>vertcut)continue; // cut on the vertex
105
106 if(part->T()>timecut)continue;
107
108 Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
109 if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
110
111 Bool_t prmch = kTRUE; // candidate primary track
112 Int_t fidau=part->GetFirstDaughter(); // cut on daughters
113 Int_t lasdau=0;
114 Int_t ndau=0;
115 if(fidau>=0){
116 lasdau=part->GetLastDaughter();
117 ndau=lasdau-fidau+1;
118 }
119 if(ndau>0){
120 for(Int_t j=fidau;j<=lasdau;j++){
121 TParticle *dau=gAlice->Particle(j);
122 Double_t distd=TMath::Sqrt((dau->Vx()-xori)*(dau->Vx()-xori)+(dau->Vy()-yori)*(dau->Vy()-yori)+(dau->Vz()-zori)*(dau->Vz()-zori));
123 if(distd<decacut)prmch=kFALSE; // eliminate if the decay is near the vertex
124 }
125 }
126
127 if(prmch){
128 nprch1++;
129 part->SetBit(kPrimaryCharged);
130 listprim->Add(part); // list of primary particles (before cleanup)
131 }
132 }
133
134
135 nprch2=0;
136 for(Int_t iprim = 0; iprim<nparticles; iprim++){ // cleanup loop
137 TParticle * part = gAlice->Particle(iprim);
138 if(part->TestBit(kPrimaryCharged)){
139 Int_t mothind=part->GetFirstMother();
140 if(mothind<0)continue;
141 TParticle *moth=gAlice->Particle(mothind);
142 TParticle *mothb=(TParticle*)listprim->FindObject(moth);
143 if(mothb){
144 listprim->Remove(moth);
145 moth->ResetBit(kPrimaryCharged);
146 nprch2++;
147 }
148 }
149 }
150
151 listprim->Clear("nodelete");
152 delete listprim;
153 return nprch1-nprch2;
154
155}