Possible to specify a number of slices, instead of only one
[u/mrichter/AliRoot.git] / HLT / hough / GetGoodParticles.cxx
CommitLineData
85da5382 1#ifndef __CINT__
2#include "alles.h"
3#include "GetGoodParticles.h"
4#include "AliL3Transform.h"
5#endif
6
7
722130ee 8void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitfile)
85da5382 9{
10
11 Int_t good_number = 70;
12
13 struct GoodTrack goodtracks[15000];
14 Int_t nt=0;
15 char *fname = "good_tracks";
16
17 ifstream in(fname);
18 if(in)
19 {
20 printf("Delete you old good_tracks file\n");
21 return;
22 }
23
24 TFile *evfile = TFile::Open(eventfile);
25 if(!evfile)
26 {
27 printf("No eventfile\n");
28 return;
29 }
30 gAlice = (AliRun*)evfile->Get("gAlice");
31 gAlice->GetEvent(0);
32
33 AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
34 AliTPCParam *param = (AliTPCParam*)evfile->Get("75x40_100x60");
35
36 TPC->SetParam(param);
37
38 Int_t ver = TPC->IsVersion();
39 if(ver!=2)
40 {
41 printf("Wrong TPC version\n");
42 return;
43 }
44
45 Int_t zero=TPC->GetParam()->GetZeroSup();
46
47 Int_t np = gAlice->GetNtrack();
48 Int_t *good = new Int_t[np];
49 for(Int_t ii=0; ii<np; ii++)
50 good[ii] = 0;
51
52 AliL3Transform *transform = new AliL3Transform();
53
54 TFile *digfile = TFile::Open(digitfile);
55 TTree *tree = (TTree*)digfile->Get("TreeD_75x40_100x60_0");
56 AliSimDigits da, *digits=&da;
57 tree->GetBranch("Segment")->SetAddress(&digits);
58
59 Int_t *count = new Int_t[np]; //np number of particles.
60 Int_t i;
61 for (i=0; i<np; i++) count[i]=0;
62 for(i=0; i<tree->GetEntries(); i++)
63 {
64 if (!tree->GetEvent(i)) continue;
65 Int_t sec,row;
66 param->AdjustSectorRow(digits->GetID(),sec,row);
67 Int_t sl,padrow;
68 transform->Sector2Slice(sl,padrow,sec,row);
722130ee 69 if(sl < minslice) continue;
70 if(sl > maxslice) break;
85da5382 71 digits->First();
72 do {
73 Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
74 Short_t dig = digits->GetDigit(it,ip);
75 Int_t idx0=digits->GetTrackID(it,ip,0);
76 Int_t idx1=digits->GetTrackID(it,ip,1);
77 Int_t idx2=digits->GetTrackID(it,ip,2);
78
79 if (idx0>=0 && dig>=zero) count[idx0]+=1;
80 if (idx1>=0 && dig>=zero) count[idx1]+=1;
81 if (idx2>=0 && dig>=zero) count[idx2]+=1;
82 } while (digits->Next());
83
84 for (Int_t j=0; j<np; j++)
85 {
86 if (count[j]>1) //at least two digits at this padrow
87 good[j]++;
88
89 count[j]=0;
90 }
91 }
92
93 delete[] count;
94
95
96 //The following code has been taken from offlinemacro->AliTPCComparison.C
97
98 TTree *TH=gAlice->TreeH();
99 Int_t npart=(Int_t)TH->GetEntries();
100
101
102
103 while (npart--) {
104 AliTPChit *hit0=0;
105
106 TPC->ResetHits();
107 TH->GetEvent(npart);
108 AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
109 while (hit){
110 if (hit->fQ==0.) break;
111 hit = (AliTPChit*) TPC->NextHit();
112 }
113 if (hit) {
114 hit0 = new AliTPChit(*hit); //Make copy of hit
115 hit = hit0;
116 }
117 else continue;
118 AliTPChit *hit1=(AliTPChit*)TPC->NextHit();
119 if (hit1==0) continue;
120 if (hit1->fQ != 0.) continue;
121 Int_t i=hit->Track();
122 TParticle *p = (TParticle*)gAlice->Particle(i);
123
124 if (p->GetFirstMother()>=0) continue; //secondary particle
125 if (good[i] < good_number) continue;
126 if (p->Pt() < 0.1) continue;
127 if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
128 printf("checking paricle %d with %d hits \n",i,good[i]);
129
130 goodtracks[nt].eta = p->Eta();
131 goodtracks[nt].label=i;
132 goodtracks[nt].code=p->GetPdgCode();
133 goodtracks[nt].px=hit->X(); goodtracks[nt].py=hit->Y(); goodtracks[nt].pz=hit->Z();
134 goodtracks[nt].pt = p->Pt();
135 goodtracks[nt].nhits = good[i];
136 nt++;
137
138 if (hit0) delete hit0;
139 }
140
141 ofstream out(fname);
142 if(out)
143 {
144 for (Int_t ngd=0; ngd<nt; ngd++)
145 out<<goodtracks[ngd].label<<' '<<goodtracks[ngd].code<<' '<<
146 goodtracks[ngd].px<<' '<<goodtracks[ngd].py<<' '<<goodtracks[ngd].pz<<' '<<
147 goodtracks[ngd].pt<<' '<<goodtracks[ngd].eta<<' '<<goodtracks[ngd].nhits<<endl;
148
149 }
150 else
151 printf("cannot open file\n");
152
153 out.close();
154
155 delete [] good;
156
157 evfile->Close();
158 digfile->Close();
159
160 delete transform;
161
162}