1 // This is an example macro which loops on all ESD events in a file set
2 // and stores the selected information in form of mini ESD tree.
3 // It also extracts the corresponding MC information in a parallel MC tree.
5 // Input: galice.root,Kinematics.root,AliESDs.root (read only)
6 // Output: miniesd.root (recreate)
8 // Usage: faster if compiled
9 // gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT/TPC");
10 // gROOT->LoadMacro("mcminiesd.C+");
13 // Author: Peter Hristov
14 // e-mail: Peter.Hristov@cern.ch
16 #if !defined(__CINT__) || defined(__MAKECINT__)
19 #include <Riostream.h>
23 #include <TStopwatch.h>
25 #include <TParticle.h>
26 #include <TLorentzVector.h>
29 #include <TDatabasePDG.h>
30 #include <TParticlePDG.h>
32 // AliRoot include files
34 #include "AliESDVertex.h"
35 #include "AliRunLoader.h"
38 #include "AliHeader.h"
39 #include "AliGenEventHeader.h"
40 #include "AliTPCtrack.h"
41 #include "AliTracker.h"
45 void selectRecTracks(AliESD* esdOut, Int_t * select, Int_t &nkeep0) {
46 // Select also all the reconstructed tracks in case it was not done before
47 Int_t nrec = esdOut->GetNumberOfTracks();
48 for(Int_t irec=0; irec<nrec; irec++) {
49 AliESDtrack * track = esdOut->GetTrack(irec);
50 UInt_t label = TMath::Abs(track->GetLabel());
51 if (label>=10000000) continue; // Underlying event. 10000000 is the
52 // value of fkMASKSTEP in AliRunDigitizer
60 void copyGeneralESDInfo(AliESD* esdIn, AliESD* esdOut) {
61 // Event and run number
62 esdOut->SetEventNumber(esdIn->GetEventNumber());
63 esdOut->SetRunNumber(esdIn->GetRunNumber());
66 esdOut->SetTriggerMask(esdIn->GetTriggerMask());
69 esdOut->SetMagneticField(esdIn->GetMagneticField());
72 const AliESDVertex * vtxIn = esdIn->GetVertex();
73 AliESDVertex * vtxOut = 0x0;
78 vtxIn->GetCovMatrix(cov);
80 vtxOut = new AliESDVertex(pos,cov,
82 vtxIn->GetNContributors());
84 vtxIn->GetTruePos(tp);
85 vtxOut->SetTruePos(tp);
88 vtxOut = new AliESDVertex();
90 esdOut->SetVertex(vtxOut);
93 void selectMiniESD(AliESD* esdIn, AliESD* &esdOut) {
94 // This function copies the general ESD information
95 // and selects the reconstructed tracks of interest
96 // The second argument is a reference to pointer since we
97 // want to operate not with the content, but with the pointer itself
99 printf("--------------------\n");
100 printf("Selecting data mini ESD\n");
104 // Create the new output ESD
105 esdOut = new AliESD();
107 // Copy the general information
108 copyGeneralESDInfo(esdIn, esdOut);
111 Int_t ntrk = esdIn->GetNumberOfTracks();
114 for (Int_t itrk=0; itrk<ntrk; itrk++) {
116 AliESDtrack * trackIn = esdIn->GetTrack(itrk);
117 UInt_t status=trackIn->GetStatus();
119 //select only tracks with TPC or ITS refit
120 if ((status&AliESDtrack::kTPCrefit)==0
121 && (status&AliESDtrack::kITSrefit)==0
124 //select only tracks with "combined" PID
125 if ((status&AliESDtrack::kESDpid)==0) continue;
127 AliESDtrack * trackOut = new AliESDtrack(*trackIn);
129 // Reset most of the redundant information
130 trackOut->MakeMiniESDtrack();
132 esdOut->AddTrack(trackOut);
134 // Do not delete trackIn, it is a second pointer to an
135 // object belonging to the TClonesArray
139 printf("%d tracks selected\n",esdOut->GetNumberOfTracks());
143 printf("--------------------\n");
146 void fixMotherDaughter(TClonesArray& particles, Int_t * map) {
147 // Fix mother-daughter relationship
148 // using the map of old to new indexes
149 Int_t nkeep = particles.GetEntriesFast();
150 for (Int_t i=0; i<nkeep; i++) {
151 TParticle * part = (TParticle*)particles[i];
154 Int_t mum = part->GetFirstMother();
155 if (mum>-1 && i>map[mum])
156 part->SetFirstMother(map[mum]);
159 Int_t fd = part->GetFirstDaughter();
160 Int_t ld = part->GetLastDaughter();
162 // invalidate daughters
163 part->SetFirstDaughter(-1);
164 part->SetLastDaughter(-1);
166 for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
168 // this daughter was selected
169 if(part->GetFirstDaughter() < 0) part->SetFirstDaughter(map[id]);
170 part->SetLastDaughter(map[id]);
176 void checkConsistency(TClonesArray& particles) {
178 // each mother is before the daughters,
179 // each daughter from the mother's list
180 // points back to its mother
181 Int_t nkeep = particles.GetEntriesFast();
182 for (Int_t i=0; i<nkeep; i++) {
183 TParticle * part = (TParticle*)particles[i];
185 Int_t mum = part->GetFirstMother();
187 if (mum>-1 && i<mum) {
188 cout << "Problem: mother " << mum << " after daughter " << i << endl;
191 Int_t fd = part->GetFirstDaughter();
192 Int_t ld = part->GetLastDaughter();
194 if (fd > ld ) cout << "Problem " << fd << " > " << ld << endl;
195 if (fd > -1 && fd < i )
196 cout << "Problem: mother " << i << " after daughter " << ld << endl;
198 for (Int_t id=TMath::Max(fd,0); id<=ld; id++) {
199 TParticle * daughter = (TParticle*)particles[id];
200 if (daughter->GetFirstMother() != i) {
201 cout << "Problem "<< i << " not "
202 << daughter->GetFirstMother() << endl;
211 void kinetree(AliRunLoader* runLoader, AliESD* &esdOut, TClonesArray* &ministack) {
213 printf("--------------------\n");
214 printf("Selecting MC mini ESD\n");
218 // Particle properties
219 TDatabasePDG * pdgdb = TDatabasePDG::Instance();
222 AliStack * stack = runLoader->Stack();
224 Int_t npart = stack->GetNtrack();
226 // Particle momentum and vertex. Will be extracted from TParticle
227 TLorentzVector momentum;
229 runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
232 // Counter for selected particles
235 Int_t * select = new Int_t[npart];
237 // Loop on particles: physics selection
238 for (Int_t ipart=0; ipart<npart; ipart++) {
242 TParticle * part = stack->Particle(ipart);
245 cout << "Empty particle " << ipart << endl;
249 // Particle momentum and vertex of origin
250 part->Momentum(momentum);
251 dvertex.SetXYZ(part->Vx()-vertex[0],
252 part->Vy()-vertex[1],
253 part->Vz()-vertex[2]);
255 // Now select only the "interesting" particles
257 // Select all particles from the primary vertex:
258 // resonances and products of resonance decays
259 if(dvertex.Mag()<0.0001) {
265 // Reject particles born outside ITS
266 if (dvertex.Perp()>100) continue;
268 // Reject particles with too low momentum, they don't rich TPC
269 if (part->Pt()<0.1) continue;
271 // Reject "final state" neutral particles
272 // This has to be redone after this loop
273 Int_t pdgcode = part->GetPdgCode();
274 TParticlePDG * pdgpart = pdgdb->GetParticle(pdgcode);
276 if (pdgpart) charge = (Int_t) pdgpart->Charge();
279 if (part->GetFirstDaughter()<0) continue;
282 // Select the rest of particles except the case
283 // particle -> particle + X
284 // for example bremstrahlung and delta electrons
285 Int_t mumid = part->GetFirstMother();
286 TParticle * mother = stack->Particle(mumid);
288 Int_t mumpdg = mother->GetPdgCode();
289 Bool_t skip = kFALSE;
290 for (Int_t id=mother->GetFirstDaughter();
291 id<=mother->GetLastDaughter();
293 TParticle * daughter=stack->Particle(id);
294 if (!daughter) continue;
295 if (mumpdg == daughter->GetPdgCode()) {
303 // All the criteria are OK, this particle is selected
307 } // end loop over particles in the current event
309 selectRecTracks(esdOut, select, nkeep0);
311 // Container for the selected particles
312 TClonesArray * pparticles = new TClonesArray("TParticle",nkeep0);
313 TClonesArray &particles = *pparticles;
315 // Map of the indexes in the stack and the new ones in the TClonesArray
316 Int_t * map = new Int_t[npart];
318 // Counter for selected particles
321 for (Int_t ipart=0; ipart<npart; ipart++) {
323 map[ipart] = -99; // Reset the map
327 TParticle * part = stack->Particle(ipart);
329 new(particles[nkeep++]) TParticle(*part);
334 if (nkeep0 != nkeep) printf("Strange %d is not %d\n",nkeep0,nkeep);
338 // Fix mother-daughter relationship
339 fixMotherDaughter(particles,map);
342 checkConsistency(particles);
344 // Now remove the "final state" neutral particles
346 TClonesArray * particles1 = new TClonesArray("TParticle",nkeep);
347 Int_t * map1 = new Int_t[nkeep];
351 for (Int_t ipart=0; ipart<nkeep; ipart++) {
353 map1[ipart] = -99; // Reset the map
355 TParticle * part = (TParticle*)particles[ipart];
357 // Reject "final state" neutral particles
358 // This has to be redone after this loop
359 Int_t pdgcode = part->GetPdgCode();
360 TParticlePDG * pdgpart = pdgdb->GetParticle(pdgcode);
362 if (pdgpart) charge = (Int_t) pdgpart->Charge();
365 if (part->GetFirstDaughter()<0) continue;
368 // Select the particle
369 map1[ipart] = nkeep1;
370 TClonesArray &ar = *particles1;
371 new(ar[nkeep1++]) TParticle(*part);
376 cout << nkeep1 << " particles finally selected" << endl;
378 fixMotherDaughter(*particles1,map1);
379 checkConsistency(*particles1);
381 // Remap the labels of reconstructed tracks
383 AliESD * esdNew = new AliESD();
385 copyGeneralESDInfo(esdOut,esdNew);
388 Int_t nrec = esdOut->GetNumberOfTracks();
389 for(Int_t irec=0; irec<nrec; irec++) {
390 AliESDtrack * track = esdOut->GetTrack(irec);
391 UInt_t label = TMath::Abs(track->GetLabel());
392 if (label<10000000) { // Signal event. 10000000 is the
393 // value of fkMASKSTEP in AliRunDigitizer
394 track->SetLabel(TMath::Sign(map1[map[label]],track->GetLabel()));
396 esdNew->AddTrack(track);
402 ministack = particles1;
411 printf("--------------------\n");
416 AliTPCtrack * particle2track(TParticle * part) {
418 // Converts TParticle to AliTPCtrack
426 // Calculate alpha: the rotation angle of the corresponding TPC sector
427 alpha = part->Phi()*180./TMath::Pi();
428 if (alpha<0) alpha+= 360.;
429 if (alpha>360) alpha -= 360.;
431 Int_t sector = (Int_t)(alpha/20.);
432 alpha = 10. + 20.*sector;
434 alpha *= TMath::Pi();
437 // Reset the covariance matrix
438 for (Int_t i=0; i<15; i++) cc[i]=0.;
442 index = part->GetUniqueID();
444 // Get the vertex of origin and the momentum
445 TVector3 ver(part->Vx(),part->Vy(),part->Vz());
446 TVector3 mom(part->Px(),part->Py(),part->Pz());
449 // Rotate to the local (sector) coordinate system
453 // X of the referense plane
457 // fX = xref X-coordinate of this track (reference plane)
458 // fAlpha = Alpha Rotation angle the local (TPC sector)
459 // fP0 = YL Y-coordinate of a track
460 // fP1 = ZG Z-coordinate of a track
461 // fP2 = C*x0 x0 is center x in rotated frame
462 // fP3 = Tgl tangent of the track momentum dip angle
463 // fP4 = C track curvature
466 TVector3 field(0.,0.,1/AliKalmanTrack::GetConvConst());
467 // radius [cm] of track projection in (x,y)
469 mom.Pt()*AliKalmanTrack::GetConvConst();
472 TDatabasePDG::Instance()->GetParticle(part->GetPdgCode())->Charge();
475 if (TMath::Abs(charge) < 0.9) charge=1.e-9;
477 TVector3 vrho = mom.Cross(field);
482 Double_t x0 = vrho.X();
486 xx[3] = mom.Pz()/mom.Pt();
489 // if (TMath::Abs(charge) < 0.9) xx[4] = 0;
491 AliTPCtrack * tr = new AliTPCtrack(index, xx, cc, xref, alpha);
499 printf("====================\n");
500 printf("Main program\n");
504 // Input "data" file, tree, and branch
505 TFile * fIn = TFile::Open("AliESDs.root");
506 TTree * tIn = (TTree*)fIn->Get("esdTree");
507 TBranch * bIn = tIn->GetBranch("ESD");
509 AliESD * esdIn = 0; // The input ESD object is put here
510 bIn->SetAddress(&esdIn);
512 // Output file, trees, and branches.
513 // Both "data" and MC are stored in the same file
514 TFile * fOut = TFile::Open("miniesd.root","recreate");
515 TTree * tOut = new TTree("esdTree","esdTree");
516 TTree * tMC = new TTree("mcStackTree","mcStackTree");
518 AliESD * esdOut = 0; // The newly created ESD object
519 TBranch * bOut = tOut->Branch("ESD","AliESD",&esdOut);
520 bOut->SetCompressionLevel(9);
522 // TParticles are stored in TClonesArray
523 // The corresponding branch is created using "dummy" TClonesArray
524 TClonesArray * ministack = new TClonesArray("TParticle");
525 TBranch * bMC = tMC->Branch("Stack",&ministack);
528 bMC->SetCompressionLevel(9);
531 AliHeader * headerIn = 0; //
532 TBranch * bHeader = tMC->Branch("Header","AliHeader",&headerIn);
533 bHeader->SetCompressionLevel(9);
536 AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
539 runLoader->LoadgAlice();
540 gAlice = runLoader->GetAliRun();
543 AliTracker::SetFieldMap(gAlice->Field(),1); // 1 means uniform magnetic field
545 // Now load kinematics and event header
546 runLoader->LoadKinematics();
547 runLoader->LoadHeader();
549 // Loop on events: check that MC and data contain the same number of events
550 Int_t nevESD = bIn->GetEntries();
551 Int_t nevMC = (Int_t)runLoader->GetNumberOfEvents();
554 cout << "Inconsistent number of ESD("<<nevESD<<") and MC("<<nevMC<<") events" << endl;
557 Int_t nev=TMath::Min(nevMC,nevESD);
559 cout << nev << " events to process" << endl;
561 for (Int_t iev=0; iev<nev; iev++) {
562 cout << "Event " << iev << endl;
567 selectMiniESD(esdIn,esdOut);
570 runLoader->GetEvent(iev);
573 kinetree(runLoader,esdOut,ministack);
574 bMC->SetAddress(&ministack); // needed because ministack has changed
577 headerIn = runLoader->GetHeader();
591 fOut = tOut->GetCurrentFile();
599 printf("====================\n");