3 #include "AliTRDtracker.h"
5 #include "AliTRDcluster.h"
8 #include "AliTRDgeometry.h"
9 #include "AliTRDparameter.h"
11 #include "AliTRDmcTrack.h"
13 #include "AliTRDtrack.h"
16 #include "TParticle.h"
17 #include "TStopwatch.h"
20 Int_t Find(Int_t prtcl, TObjArray *mctarray) {
22 // Returns index of the mct which corresponds to particle <prtcl>
24 Int_t b=0, e=mctarray->GetEntriesFast(), m=(b+e)/2;
25 AliTRDmcTrack *mct = 0;
26 while ((b<e)&&(e!=0)) {
28 mct = (AliTRDmcTrack*) mctarray->UncheckedAt(m);
29 if (prtcl > mct->GetTrackIndex()) b=m+1;
33 if(mct->GetTrackIndex() == prtcl) return m;
35 if((m+1) < mctarray->GetEntriesFast()) {
36 mct = (AliTRDmcTrack*) mctarray->UncheckedAt(m+1);
37 if(mct->GetTrackIndex() == prtcl) return m+1;
43 void AliTRDsaveMCtracks() {
45 TObjArray mctracks(2000);
46 TObjArray *TRDmcTracks = &mctracks;
48 Char_t *alifile = "galice.root";
50 Float_t low_pt_cut = 0.05;
51 Float_t low_p_cut = 0.02;
52 Int_t min_no_of_clusters = 10;
54 // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
55 TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
57 cout << "Open the ALIROOT-file " << alifile << endl;
58 gafl = new TFile(alifile);
61 cout << alifile << " is already open" << endl;
64 // Get AliRun object from file or create it if not on file
65 gAlice = (AliRun*) gafl->Get("gAlice");
67 cout << "AliRun object found on file" << endl;
69 gAlice = new AliRun("gAlice","Alice test program");
71 AliTRDv1 *fTRD = (AliTRDv1*) gAlice->GetDetector("TRD");
73 // Import the Trees for the event nEvent in the file
74 const Int_t nparticles = gAlice->GetEvent(nEvent);
76 printf("found %d particles in event %d \n", nparticles, nEvent);
78 // Create TRDmcTracks for each particle
79 Int_t label = -1, charge = 0;
86 for (Int_t ii=0; ii<nparticles; ii++) {
88 printf("\r particle %d out of %d",ii,nparticles);
90 TParticle *p = gAlice->Particle(ii);
91 if(p->P() < low_p_cut) continue;
94 if (p->GetFirstMother()>=0) primary = kFALSE;
96 pdg_code = (Int_t) p->GetPdgCode();
98 if ((pdg_code == 10010020) ||
99 (pdg_code == 10010030) ||
100 (pdg_code == 50000050) ||
101 (pdg_code == 50000051) ||
102 (pdg_code == 10020040)) {
108 TParticlePDG *pdg = p->GetPDG();
109 charge = (Int_t) pdg->Charge();
112 if(charge == 0) continue;
114 AliTRDmcTrack *t = new AliTRDmcTrack(ii, primary, mass, charge, pdg_code);
115 TRDmcTracks->AddLast(t);
119 // Loop through TRD clusters and assign indexes to MC tracks
121 TFile *geofile =TFile::Open("AliTRDclusters.root");
122 AliTRDtracker *Tracker = new AliTRDtracker(geofile);
123 Tracker->SetEventNumber(nEvent);
125 AliTRDgeometry *fGeo = (AliTRDgeometry*) geofile->Get("TRDgeometry");
126 AliTRDparameter *fPar = (AliTRDparameter*) geofile->Get("TRDparameter");
129 Char_t *clusterfile = "AliTRDclusters.root";
130 TObjArray carray(2000);
131 TObjArray *ClustersArray = &carray;
132 Tracker->ReadClusters(ClustersArray,clusterfile);
134 Int_t nClusters = carray.GetEntriesFast();
137 AliTRDmcTrack *tpr = NULL;
139 for (Int_t i = 0; i < nClusters; i++) {
141 printf("\r assigning cluster %d out of %d", i, nClusters);
142 AliTRDcluster *cl = (AliTRDcluster *) ClustersArray->UncheckedAt(i);
144 for(Int_t j=0; j<3; j++) {
145 track = cl->GetLabel(j);
146 if(track < 0) continue;
147 if(track >= nparticles)
148 printf("Track index %d is larger than total number of particles %d\n"
151 index = Find(track, TRDmcTracks);
153 if(index > 0) tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
156 label = tpr->GetTrackIndex();
157 if(label != track) printf("Got label %d while expected %d !\n",
159 TParticle *p = gAlice->Particle(track);
161 if (p->Pt() > low_pt_cut) tpr->Update(i);
167 // Loop through the TRD hits and set XYZ and Pin and Pout
169 TTree *hitTree = gAlice->TreeH();
171 Int_t preamp, amp, det, plane = 0, nBytes = 0;
172 Int_t nTrack = (Int_t) hitTree->GetEntries();
173 Float_t pos[3], rot[3];
174 Float_t prepos[3], prerot[3];
176 Bool_t entrance = kFALSE;
177 Bool_t exit = kFALSE;
180 Double_t xin[6], xout[6];
181 for(Int_t j = 0; j < 6; j++) {
182 xin[j] = Tracker->GetX(0,j,14);
183 xout[j] = Tracker->GetX(0,j,0);
188 Int_t current_track = -1;
194 nBytes += hitTree->GetEvent(nTrack);
198 for(Int_t j = 0; j < 3; j++) { pos[j] = 1111.; rot[j] = 1111.;}
200 for(hit = (AliTRDhit *) fTRD->FirstHit(-1);
202 hit = (AliTRDhit *) fTRD->NextHit()) {
205 for(Int_t j = 0; j < 3; j++) { prepos[j] = pos[j]; prerot[j] = rot[j];}
207 amp = hit->GetCharge();
208 if(amp != 0) continue;
209 track = hit->GetTrack();
211 TParticle *p = gAlice->Particle(track);
212 if (p->Pt() <= low_pt_cut) continue;
214 if(track != current_track) {
215 current_track = track;
219 det = hit->GetDetector();
220 plane = fGeo->GetPlane(det);
224 fGeo->Rotate(det,pos,rot);
226 if(TMath::Abs(rot[0] - xin[plane]) < 1.2) entrance = kTRUE;
227 else entrance = kFALSE;
228 if(TMath::Abs(rot[0] - xout[plane]) < 1.2) exit = kTRUE;
231 if(entrance && (preamp == 0) && (prerot[0] < 200)) {
232 index = Find(track, TRDmcTracks);
235 tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
237 label = tpr->GetTrackIndex();
238 if(label != track) printf("Got label %d while expected %d !\n",
242 printf("\n momentum at plane %d entrance, Pxyz: %f, %f, %f\n",
243 plane, prepos[0], prepos[1], prepos[2]);
244 printf(" XYZ at plane %d entrance: %f, %f, %f\n",
245 plane, rot[0], rot[1], rot[2]);
246 printf("expected x range: %f .. %f\n",xin[plane]-1.2,xin[plane]+1.2);
248 tpr->SetPin(plane, prepos[0], prepos[1], prepos[2]);
249 tpr->SetXYZin(plane, rot[0], rot[1], rot[2]);
256 if(exit && (preamp == 0) && (prerot[0] < 200)) {
257 index = Find(track, TRDmcTracks);
260 tpr = (AliTRDmcTrack *) TRDmcTracks->UncheckedAt(index);
262 label = tpr->GetTrackIndex();
263 if(label != track) printf("Got label %d while expected %d !\n",
267 printf("\n momentum at plane %d exit, Pxyz: %f, %f, %f\n",
268 plane, prepos[0], prepos[1], prepos[2]);
269 printf(" XYZ at plane %d exit: %f, %f, %f\n",
270 plane, rot[0], rot[1], rot[2]);
271 printf("expected x range: %f .. %f\n",xout[plane]-1.2,xout[plane]+1.2);
273 tpr->SetPout(plane, prepos[0], prepos[1], prepos[2]);
274 tpr->SetXYZout(plane, rot[0], rot[1], rot[2]);
283 // Write TRDmcTracks in output file
285 TDirectory *savedir=gDirectory;
286 Char_t *filename = "AliTRDmcTracks.root";
287 TFile *out = new TFile(filename,"RECREATE");
289 TTree *tracktree = new TTree("MCtracks","TRD MC tracks");
291 AliTRDmcTrack *iotrack=0;
292 tracktree->Branch("MCtracks","AliTRDmcTrack",&iotrack,32000,0);
294 Int_t ntracks = TRDmcTracks->GetEntriesFast();
296 for (Int_t i=0; i<ntracks; i++) {
297 AliTRDmcTrack *pt=(AliTRDmcTrack*)TRDmcTracks->UncheckedAt(i);
299 Int_t n = pt->GetNumberOfClusters();
300 if(n > min_no_of_clusters) {
303 printf("Put track with label %d and %d clusters in the output tree \n",
304 pt->GetTrackIndex(),n);