3 #include "GetGoodParticles.h"
4 #include "AliL3Transform.h"
8 void GetGoodParticles(Int_t minslice,Int_t maxslice,char *eventfile,char *digitfile,Int_t event)
11 Int_t good_number = 10; //Minimum number of points on a good track
13 struct GoodTrack goodtracks[15000];
15 char *fname = "good_tracks";
20 printf("Delete you old good_tracks file\n");
24 TFile *evfile = TFile::Open(eventfile);
27 printf("No eventfile\n");
30 gAlice = (AliRun*)evfile->Get("gAlice");
31 gAlice->GetEvent(event);
33 AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC");
34 AliTPCParam *param = (AliTPCParam*)evfile->Get("75x40_100x60");
38 Int_t ver = TPC->IsVersion();
41 printf("Wrong TPC version\n");
45 Int_t zero=TPC->GetParam()->GetZeroSup();
47 Int_t np = gAlice->GetNtrack();
48 Int_t *good = new Int_t[np];
49 for(Int_t ii=0; ii<np; ii++)
52 AliL3Transform *transform = new AliL3Transform();
54 TFile *digfile = TFile::Open(digitfile);
56 sprintf(dname,"TreeD_75x40_100x60_%d",event);
57 TTree *tree=(TTree*)digfile->Get(dname);
58 AliSimDigits da, *digits=&da;
59 tree->GetBranch("Segment")->SetAddress(&digits);
61 Int_t *count = new Int_t[np]; //np number of particles.
63 for (i=0; i<np; i++) count[i]=0;
64 for(i=0; i<tree->GetEntries(); i++)
66 if (!tree->GetEvent(i)) continue;
68 param->AdjustSectorRow(digits->GetID(),sec,row);
70 transform->Sector2Slice(sl,padrow,sec,row);
71 if(sl < minslice) continue;
72 if(sl > maxslice) break;
75 Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
76 Short_t dig = digits->GetDigit(it,ip);
77 Int_t idx0=digits->GetTrackID(it,ip,0);
78 Int_t idx1=digits->GetTrackID(it,ip,1);
79 Int_t idx2=digits->GetTrackID(it,ip,2);
81 if (idx0>=0 && dig>=zero) count[idx0]+=1;
82 if (idx1>=0 && dig>=zero) count[idx1]+=1;
83 if (idx2>=0 && dig>=zero) count[idx2]+=1;
84 } while (digits->Next());
86 for (Int_t j=0; j<np; j++)
88 if (count[j]>1) //at least two digits at this padrow
98 //The following code has been taken from offlinemacro->AliTPCComparison.C
100 TTree *TH=gAlice->TreeH();
101 Int_t npart=(Int_t)TH->GetEntries();
102 cout<<"Found "<<npart<<" particles in this event"<<endl;
110 AliTPChit * hit = (AliTPChit*) TPC->FirstHit(-1);
112 if (hit->fQ==0.) break;
113 hit = (AliTPChit*) TPC->NextHit();
116 hit0 = new AliTPChit(*hit); //Make copy of hit
120 AliTPChit *hit1=(AliTPChit*)TPC->NextHit();
121 if (hit1==0) continue;
122 if (hit1->fQ != 0.) continue;
123 Int_t i=hit->Track();
124 TParticle *p = (TParticle*)gAlice->Particle(i);
126 if (p->GetFirstMother()>=0) continue; //secondary particle
127 if (good[i] < good_number) continue;
128 if (p->Pt() < 0.1) continue;
129 if (TMath::Abs(p->Pz()/p->Pt())>0.999) continue;
130 cout<<"Checking particle with "<<good[i]<<" hits, and pt: "<<p->Pt()<<endl;
131 goodtracks[nt].eta = p->Eta();
132 goodtracks[nt].label=i;
133 goodtracks[nt].code=p->GetPdgCode();
134 goodtracks[nt].px=hit->X(); goodtracks[nt].py=hit->Y(); goodtracks[nt].pz=hit->Z();
135 goodtracks[nt].pt = p->Pt();
136 goodtracks[nt].nhits = good[i];
139 if (hit0) delete hit0;
141 cout<<"Found "<<nt<<" good tracks in this event"<<endl;
145 for (Int_t ngd=0; ngd<nt; ngd++)
146 out<<goodtracks[ngd].label<<' '<<goodtracks[ngd].code<<' '<<
147 goodtracks[ngd].px<<' '<<goodtracks[ngd].py<<' '<<goodtracks[ngd].pz<<' '<<
148 goodtracks[ngd].pt<<' '<<goodtracks[ngd].eta<<' '<<goodtracks[ngd].nhits<<endl;
152 printf("cannot open file\n");