426757ed78a6a8ef4c05f9e58a8632ca945c8f92
[u/mrichter/AliRoot.git] / HLT / hough / GetGoodParticles.cxx
1 #ifndef __CINT__
2 #include "alles.h"
3 #include "GetGoodParticles.h"
4 #include "AliL3Transform.h"
5 #endif
6
7
8 void GetGoodParticles(Int_t slice,char *eventfile,char *digitfile)
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);
69       if(sl < slice) continue;
70       if(sl != slice) break;
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 }