]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliGenInfo.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New comparison macros (M.Ivanov)
[u/mrichter/AliRoot.git] / STEER / AliGenInfo.C
CommitLineData
ae7d73d2 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16
17///////////////////////////////////////////////////////////////////////////
18/*
19
20Origin: marian.ivanov@cern.ch
21
22Macro to generate comples MC information - used for Comparison later on
23How to use it?
24
25.L $ALICE_ROOT/STEER/AliGenInfo.C+
26AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
27t->Exec();
28
29*/
30
31#if !defined(__CINT__) || defined(__MAKECINT__)
32#include <stdio.h>
33#include <string.h>
34//ROOT includes
35#include "Rtypes.h"
36#include "TFile.h"
37#include "TTree.h"
38#include "TChain.h"
39#include "TCut.h"
40#include "TString.h"
41#include "TBenchmark.h"
42#include "TStopwatch.h"
43#include "TParticle.h"
44#include "TSystem.h"
45#include "TTimer.h"
46#include "TVector3.h"
47#include "TH1F.h"
48#include "TH2F.h"
49#include "TCanvas.h"
50#include "TPad.h"
51#include "TF1.h"
52
53//ALIROOT includes
54#include "AliRun.h"
55#include "AliStack.h"
56#include "AliSimDigits.h"
57#include "AliTPCParam.h"
58#include "AliTPC.h"
59#include "AliTPCLoader.h"
60#include "AliDetector.h"
61#include "AliTrackReference.h"
62#include "AliTPCParamSR.h"
63#include "AliTracker.h"
64#include "AliMagF.h"
65#endif
66#include "AliGenInfo.h"
67//
68//
69
70AliTPCParam * GetTPCParam(){
71 AliTPCParamSR * par = new AliTPCParamSR;
72 par->Update();
73 return par;
74}
75
76
77//_____________________________________________________________________________
78Float_t TPCBetheBloch(Float_t bg)
79{
80 //
81 // Bethe-Bloch energy loss formula
82 //
83 const Double_t kp1=0.76176e-1;
84 const Double_t kp2=10.632;
85 const Double_t kp3=0.13279e-4;
86 const Double_t kp4=1.8631;
87 const Double_t kp5=1.9479;
88
89 Double_t dbg = (Double_t) bg;
90
91 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
92
93 Double_t aa = TMath::Power(beta,kp4);
94 Double_t bb = TMath::Power(1./dbg,kp5);
95
96 bb=TMath::Log(kp3+bb);
97
98 return ((Float_t)((kp2-aa-bb)*kp1/aa));
99}
100
101
102
103
104
105////////////////////////////////////////////////////////////////////////
106AliMCInfo::AliMCInfo()
107{
108 fTPCReferences = new TClonesArray("AliTrackReference",10);
109 fITSReferences = new TClonesArray("AliTrackReference",10);
110 fTRDReferences = new TClonesArray("AliTrackReference",10);
111 fTRdecay.SetTrack(-1);
112}
113
114AliMCInfo::~AliMCInfo()
115{
116 if (fTPCReferences) {
117 delete fTPCReferences;
118 }
119 if (fITSReferences){
120 delete fITSReferences;
121 }
122 if (fTRDReferences){
123 delete fTRDReferences;
124 }
125}
126
127
128
129void AliMCInfo::Update()
130{
131 //
132 //
133 fMCtracks =1;
134 if (!fTPCReferences) {
135 fNTPCRef =0;
136 return;
137 }
138 Float_t direction=1;
139 //Float_t rlast=0;
140 fNTPCRef = fTPCReferences->GetEntriesFast();
141 fNITSRef = fITSReferences->GetEntriesFast();
142
143 for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
144 AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
145 //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
146 Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
147 if (iref==0) direction = newdirection;
148 if ( newdirection*direction<0){
149 //changed direction
150 direction = newdirection;
151 fMCtracks+=1;
152 }
153 //rlast=r;
154 }
155 //
156 // decay info
157 fTPCdecay=kFALSE;
158 if (fTRdecay.GetTrack()>0){
159 fDecayCoord[0] = fTRdecay.X();
160 fDecayCoord[1] = fTRdecay.Y();
161 fDecayCoord[2] = fTRdecay.Z();
162 if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
163 fTPCdecay=kTRUE;
164 }
165 else{
166 fDecayCoord[0] = 0;
167 fDecayCoord[1] = 0;
168 fDecayCoord[2] = 0;
169 }
170 }
171 //
172 //
173 //digits information update
174 fRowsWithDigits = fTPCRow.RowsOn();
175 fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
176 fRowsTrackLength = fTPCRow.Last() - fTPCRow.First();
177 //
178 //
179 // calculate primary ionization per cm
180 if (fParticle.GetPDG()){
181 fMass = fParticle.GetMass();
182 fCharge = fParticle.GetPDG()->Charge();
183 if (fTPCReferences->GetEntriesFast()>0){
184 fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
185 }
186 if (fMass>0){
187 Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
188 fTrackRef.Py()*fTrackRef.Py()+
189 fTrackRef.Pz()*fTrackRef.Pz());
190 if (p>0.001){
191 Float_t betagama = p /fMass;
192 fPrim = TPCBetheBloch(betagama);
193 }else fPrim=0;
194 }
195 }else{
196 fMass =0;
197 fPrim =0;
198 }
199}
200
201/////////////////////////////////////////////////////////////////////////////////
202void AliGenV0Info::Update()
203{
204 fMCPd[0] = fMCd.fParticle.Px();
205 fMCPd[1] = fMCd.fParticle.Py();
206 fMCPd[2] = fMCd.fParticle.Pz();
207 fMCPd[3] = fMCd.fParticle.P();
208 fMCX[0] = fMCd.fParticle.Vx();
209 fMCX[1] = fMCd.fParticle.Vy();
210 fMCX[2] = fMCd.fParticle.Vz();
211 fMCR = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
212 fPdg[0] = fMCd.fParticle.GetPdgCode();
213 fPdg[1] = fMCm.fParticle.GetPdgCode();
214 //
215 fLab[0] = fMCd.fParticle.GetUniqueID();
216 fLab[1] = fMCm.fParticle.GetUniqueID();
217
218}
219
220
221
222
223
224////////////////////////////////////////////////////////////////////////
225
226////////////////////////////////////////////////////////////////////////
227//
228// End of implementation of the class AliMCInfo
229//
230////////////////////////////////////////////////////////////////////////
231
232
233
234////////////////////////////////////////////////////////////////////////
235digitRow::digitRow()
236{
237 Reset();
238}
239////////////////////////////////////////////////////////////////////////
240digitRow & digitRow::operator=(const digitRow &digOld)
241{
242 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
243 return (*this);
244}
245////////////////////////////////////////////////////////////////////////
246void digitRow::SetRow(Int_t row)
247{
248 if (row >= 8*kgRowBytes) {
249 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
250 return;
251 }
252 Int_t iC = row/8;
253 Int_t iB = row%8;
254 SETBIT(fDig[iC],iB);
255}
256
257////////////////////////////////////////////////////////////////////////
258Bool_t digitRow::TestRow(Int_t row)
259{
260//
261// return kTRUE if row is on
262//
263 Int_t iC = row/8;
264 Int_t iB = row%8;
265 return TESTBIT(fDig[iC],iB);
266}
267////////////////////////////////////////////////////////////////////////
268Int_t digitRow::RowsOn(Int_t upto)
269{
270//
271// returns number of rows with a digit
272// count only rows less equal row number upto
273//
274 Int_t total = 0;
275 for (Int_t i = 0; i<kgRowBytes; i++) {
276 for (Int_t j = 0; j < 8; j++) {
277 if (i*8+j > upto) return total;
278 if (TESTBIT(fDig[i],j)) total++;
279 }
280 }
281 return total;
282}
283////////////////////////////////////////////////////////////////////////
284void digitRow::Reset()
285{
286//
287// resets all rows to zero
288//
289 for (Int_t i = 0; i<kgRowBytes; i++) {
290 fDig[i] <<= 8;
291 }
292}
293////////////////////////////////////////////////////////////////////////
294Int_t digitRow::Last()
295{
296//
297// returns the last row number with a digit
298// returns -1 if now digits
299//
300 for (Int_t i = kgRowBytes-1; i>=0; i--) {
301 for (Int_t j = 7; j >= 0; j--) {
302 if TESTBIT(fDig[i],j) return i*8+j;
303 }
304 }
305 return -1;
306}
307////////////////////////////////////////////////////////////////////////
308Int_t digitRow::First()
309{
310//
311// returns the first row number with a digit
312// returns -1 if now digits
313//
314 for (Int_t i = 0; i<kgRowBytes; i++) {
315 for (Int_t j = 0; j < 8; j++) {
316 if (TESTBIT(fDig[i],j)) return i*8+j;
317 }
318 }
319 return -1;
320}
321
322////////////////////////////////////////////////////////////////////////
323//
324// end of implementation of a class digitRow
325//
326////////////////////////////////////////////////////////////////////////
327
328////////////////////////////////////////////////////////////////////////
329AliGenInfoMaker::AliGenInfoMaker()
330{
331 Reset();
332}
333
334////////////////////////////////////////////////////////////////////////
335AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
336 Int_t nEvents, Int_t firstEvent)
337{
338 Reset();
339 fFirstEventNr = firstEvent;
340 fEventNr = firstEvent;
341 fNEvents = nEvents;
342 // fFnRes = fnRes;
343 sprintf(fFnRes,"%s",fnRes);
344 //
345 fLoader = AliRunLoader::Open(fnGalice);
346 if (gAlice){
347 delete gAlice->GetRunLoader();
348 delete gAlice;
349 gAlice = 0x0;
350 }
351 if (fLoader->LoadgAlice()){
352 cerr<<"Error occured while l"<<endl;
353 }
354 Int_t nall = fLoader->GetNumberOfEvents();
355 if (nEvents==0) {
356 nEvents =nall;
357 fNEvents=nall;
358 fFirstEventNr=0;
359 }
360
361 if (nall<=0){
362 cerr<<"no events available"<<endl;
363 fEventNr = 0;
364 return;
365 }
366 if (firstEvent+nEvents>nall) {
367 fEventNr = nall-firstEvent;
368 cerr<<"restricted number of events availaible"<<endl;
369 }
370 AliMagF * magf = gAlice->Field();
371 AliTracker::SetFieldMap(magf);
372}
373
374
375AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
376{
377 //
378 if (i<fNParticles) {
379 if (fGenInfo[i]) return fGenInfo[i];
380 fGenInfo[i] = new AliMCInfo;
381 fNInfos++;
382 return fGenInfo[i];
383 }
384 else
385 return 0;
386}
387
388////////////////////////////////////////////////////////////////////////
389void AliGenInfoMaker::Reset()
390{
391 fEventNr = 0;
392 fNEvents = 0;
393 fTreeGenTracks = 0;
394 fFileGenTracks = 0;
395 fGenInfo = 0;
396 fNInfos = 0;
397 //
398 //
399 fDebug = 0;
400 fVPrim[0] = -1000.;
401 fVPrim[1] = -1000.;
402 fVPrim[2] = -1000.;
403 fParamTPC = 0;
404}
405////////////////////////////////////////////////////////////////////////
406AliGenInfoMaker::~AliGenInfoMaker()
407{
408
409 if (fLoader){
410 fLoader->UnloadgAlice();
411 gAlice = 0;
412 delete fLoader;
413 }
414}
415
416Int_t AliGenInfoMaker::SetIO()
417{
418 //
419 //
420 CreateTreeGenTracks();
421 if (!fTreeGenTracks) return 1;
422 // AliTracker::SetFieldFactor();
423
424 fParamTPC = GetTPCParam();
425 //
426 return 0;
427}
428
429////////////////////////////////////////////////////////////////////////
430Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
431{
432 //
433 //
434 // SET INPUT
435 fLoader->SetEventNumber(eventNr);
436 //
437 fLoader->LoadHeader();
438 fLoader->LoadKinematics();
439 fStack = fLoader->Stack();
440 //
441 fLoader->LoadTrackRefs();
442 fTreeTR = fLoader->TreeTR();
443 //
444 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
445 tpcl->LoadDigits();
446 fTreeD = tpcl->TreeD();
447 return 0;
448}
449
450Int_t AliGenInfoMaker::CloseIOEvent()
451{
452 fLoader->UnloadHeader();
453 fLoader->UnloadKinematics();
454 fLoader->UnloadTrackRefs();
455 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
456 tpcl->UnloadDigits();
457 return 0;
458}
459
460Int_t AliGenInfoMaker::CloseIO()
461{
462 fLoader->UnloadgAlice();
463 return 0;
464}
465
466
467
468////////////////////////////////////////////////////////////////////////
469Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
470{
471 fNEvents = nEvents;
472 fFirstEventNr = firstEventNr;
473 return Exec();
474}
475
476////////////////////////////////////////////////////////////////////////
477Int_t AliGenInfoMaker::Exec()
478{
479 TStopwatch timer;
480 timer.Start();
481 Int_t status =SetIO();
482 if (status>0) return status;
483 //
484
485 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
486 fEventNr++) {
487 SetIO(fEventNr);
488 fNParticles = fStack->GetNtrack();
489 //
490 fGenInfo = new AliMCInfo*[fNParticles];
491 for (UInt_t i = 0; i<fNParticles; i++) {
492 fGenInfo[i]=0;
493 }
494 //
495 cout<<"Start to process event "<<fEventNr<<endl;
496 cout<<"\tfNParticles = "<<fNParticles<<endl;
497 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
498 if (TreeTRLoop()>0) return 1;
499 //
500 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
501 if (TreeDLoop()>0) return 1;
502 //
503 if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
504 if (TreeKLoop()>0) return 1;
505 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
506 for (UInt_t i = 0; i<fNParticles; i++) {
507 if (fGenInfo[i]) delete fGenInfo[i];
508 }
509 delete []fGenInfo;
510 CloseIOEvent();
511 }
512 //
513 CloseIO();
514 CloseOutputFile();
515
516 cerr<<"Exec finished"<<endl;
517
518 timer.Stop();
519 timer.Print();
520 return 0;
521}
522////////////////////////////////////////////////////////////////////////
523void AliGenInfoMaker::CreateTreeGenTracks()
524{
525 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
526 if (!fFileGenTracks) {
527 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
528 return;
529 }
530 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
531 AliMCInfo * info = new AliMCInfo;
532 //
533 fTreeGenTracks->Branch("MC","AliMCInfo",&info);
534 delete info;
535 fTreeGenTracks->AutoSave();
536}
537////////////////////////////////////////////////////////////////////////
538void AliGenInfoMaker::CloseOutputFile()
539{
540 if (!fFileGenTracks) {
541 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
542 return;
543 }
544 fFileGenTracks->cd();
545 fTreeGenTracks->Write();
546 delete fTreeGenTracks;
547 fFileGenTracks->Close();
548 delete fFileGenTracks;
549 return;
550}
551
552////////////////////////////////////////////////////////////////////////
553Int_t AliGenInfoMaker::TreeKLoop()
554{
555//
556// open the file with treeK
557// loop over all entries there and save information about some tracks
558//
559
560 AliStack * stack = fStack;
561 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
562
563 if (fDebug > 0) {
564 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
565 <<fEventNr<<endl;
566 }
567 Int_t ipdg = 0;
568 TParticlePDG *ppdg = 0;
569 // not all generators give primary vertex position. Take the vertex
570 // of the particle 0 as primary vertex.
571 TDatabasePDG pdg; //get pdg table
572 //thank you very much root for this
573 TBranch * br = fTreeGenTracks->GetBranch("MC");
574 TParticle *particle = stack->ParticleFromTreeK(0);
575 fVPrim[0] = particle->Vx();
576 fVPrim[1] = particle->Vy();
577 fVPrim[2] = particle->Vz();
578 for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
579 // load only particles with TR
580 AliMCInfo * info = GetInfo(iParticle);
581 if (!info) continue;
582 //////////////////////////////////////////////////////////////////////
583 info->fLabel = iParticle;
584 //
585 info->fParticle = *(stack->Particle(iParticle));
586 info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
587 info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
588 info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
589 info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
590 info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
591 //
592 //
593 ipdg = info->fParticle.GetPdgCode();
594 info->fPdg = ipdg;
595 ppdg = pdg.GetParticle(ipdg);
596 info->fEventNr = fEventNr;
597 info->Update();
598 //////////////////////////////////////////////////////////////////////
599 br->SetAddress(&info);
600 fTreeGenTracks->Fill();
601 }
602 fTreeGenTracks->AutoSave();
603 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
604 return 0;
605}
606
607
608
609
610////////////////////////////////////////////////////////////////////////
611Int_t AliGenInfoMaker::TreeDLoop()
612{
613 //
614 // open the file with treeD
615 // loop over all entries there and save information about some tracks
616 //
617
618 Int_t nInnerSector = fParamTPC->GetNInnerSector();
619 Int_t rowShift = 0;
620 Int_t zero=fParamTPC->GetZeroSup()+6;
621 //
622 //
623 AliSimDigits digitsAddress, *digits=&digitsAddress;
624 fTreeD->GetBranch("Segment")->SetAddress(&digits);
625
626 Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
627 if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
628 for (Int_t i=0; i<sectorsByRows; i++) {
629 if (!fTreeD->GetEvent(i)) continue;
630 Int_t sec,row;
631 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
632 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
633 // here I expect that upper sectors follow lower sectors
634 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
635 //
636 digits->ExpandTrackBuffer();
637 digits->First();
638 do {
639 Int_t iRow=digits->CurrentRow();
640 Int_t iColumn=digits->CurrentColumn();
641 Short_t digitValue = digits->CurrentDigit();
642 if (digitValue >= zero) {
643 Int_t label;
644 for (Int_t j = 0; j<3; j++) {
645 // label = digits->GetTrackID(iRow,iColumn,j);
646 label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
647 if (label >= (Int_t)fNParticles) { //don't label from bakground event
648 continue;
649 }
650 if (label >= 0 && label <= (Int_t)fNParticles) {
651 if (fDebug > 6 ) {
652 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
653 <<sec<<" "
654 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
655 <<" "<<row<<endl;
656 }
657 AliMCInfo * info = GetInfo(label);
658 if (info){
659 info->fTPCRow.SetRow(row+rowShift);
660 }
661 }
662 }
663 }
664 } while (digits->Next());
665 }
666
667 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
668 return 0;
669}
670
671
672////////////////////////////////////////////////////////////////////////
673Int_t AliGenInfoMaker::TreeTRLoop()
674{
675 //
676 // loop over TrackReferences and store the first one for each track
677 //
678 TTree * treeTR = fTreeTR;
679 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
680 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
681 //
682 //
683 //track references for TPC
684 TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
685 TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
686 TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
687 TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
688 //
689 if (treeTR->GetBranch("TPC")) treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
690 if (treeTR->GetBranch("ITS")) treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
691 if (treeTR->GetBranch("TRD")) treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
692 if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
693 //
694 //
695 //
696 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
697 treeTR->GetEntry(iPrimPart);
698 //
699 // Loop over TPC references
700 //
701 for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
702 AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
703 //
704 if (trackRef->Pt()<fgTPCPtCut) continue;
705 Int_t label = trackRef->GetTrack();
706 AliMCInfo * info = GetInfo(label);
707 if (!info) info = MakeInfo(label);
708 if (!info) continue;
709 TClonesArray & arr = *(info->fTPCReferences);
710 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
711 }
712 //
713 // Loop over ITS references
714 //
715 for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
716 AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);
717 //
718 if (trackRef->Pt()<fgTPCPtCut) continue;
719 Int_t label = trackRef->GetTrack();
720 AliMCInfo * info = GetInfo(label);
721 if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
722 if (!info) info = MakeInfo(label);
723 if (!info) continue;
724 TClonesArray & arr = *(info->fITSReferences);
725 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
726 }
727 //
728 // Loop over TRD references
729 //
730 for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
731 AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);
732 //
733 if (trackRef->Pt()<fgTPCPtCut) continue;
734 Int_t label = trackRef->GetTrack();
735 AliMCInfo * info = GetInfo(label);
736 if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
737 if (!info) info = MakeInfo(label);
738 if (!info) continue;
739 TClonesArray & arr = *(info->fTRDReferences);
740 new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
741 }
742 //
743 // get dacay position
744 for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
745 AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
746 //
747 Int_t label = trackRef->GetTrack();
748 AliMCInfo * info = GetInfo(label);
749 if (!info) continue;
750 if (!trackRef->TestBit(BIT(2))) continue; //if not decay
751 // if (TMath::Abs(trackRef.X());
752 info->fTRdecay = *trackRef;
753 }
754 }
755 //
756 TPCArrayTR->Delete();
757 delete TPCArrayTR;
758 TRDArrayTR->Delete();
759 delete TRDArrayTR;
760 ITSArrayTR->Delete();
761 delete ITSArrayTR;
762 RunArrayTR->Delete();
763 delete RunArrayTR;
764 //
765 return 0;
766}
767
768////////////////////////////////////////////////////////////////////////
769Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
770 AliTPCParam *paramTPC) {
771
772 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
773 Int_t index[4];
774 paramTPC->Transform0to1(x,index);
775 paramTPC->Transform1to2(x,index);
776 return x[0];
777}
778////////////////////////////////////////////////////////////////////////
779
780
781
782TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection,
783 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
784{
785 //
786 Double_t* bins = CreateLogBins(nbins, minx, maxx);
787 Int_t nBinsRes = 30;
788 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
789 char cut[1000];
790 sprintf(cut,"%s&&%s",selection,quality);
791 char expression[1000];
792 sprintf(expression,"%s:%s>>hRes2",chy,chx);
793 fTree->Draw(expression, cut, "groff");
794 TH1F* hMean=0;
795 TH1F* hRes = CreateResHisto(hRes2, &hMean);
796 AliLabelAxes(hRes, chx, chy);
797 //
798 delete hRes2;
799 delete[] bins;
800 fRes = hRes;
801 fMean = hMean;
802 return hRes;
803}
804
805TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
806 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
807{
808 //
809 Double_t* bins = CreateLogBins(nbins, minx, maxx);
810 Int_t nBinsRes = 30;
811 TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
812 char cut[1000];
813 sprintf(cut,"%s&&%s",selection,quality);
814 char expression[1000];
815 sprintf(expression,"%s:%s>>hRes2",chy,chx);
816 fTree->Draw(expression, cut, "groff");
817 TH1F* hMean=0;
818 TH1F* hRes = CreateResHisto(hRes2, &hMean);
819 AliLabelAxes(hRes, chx, chy);
820 //
821 delete hRes2;
822 delete[] bins;
823 fRes = hRes;
824 fMean = hMean;
825 return hRes;
826}
827
828///////////////////////////////////////////////////////////////////////////////////
829///////////////////////////////////////////////////////////////////////////////////
830TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality,
831 Int_t nbins, Float_t min, Float_t max)
832{
833 //
834 //
835 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
836 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
837 char inputGen[1000];
838 sprintf(inputGen,"%s>>hGen", variable);
839 fTree->Draw(inputGen, selection, "groff");
840 char selectionRec[256];
841 sprintf(selectionRec, "%s && %s", selection, quality);
842 char inputRec[1000];
843 sprintf(inputRec,"%s>>hRec", variable);
844 fTree->Draw(inputRec, selectionRec, "groff");
845 //
846 TH1F* hEff = CreateEffHisto(hGen, hRec);
847 AliLabelAxes(hEff, variable, "#epsilon [%]");
848 fRes = hEff;
849 delete hRec;
850 delete hGen;
851 return hEff;
852}
853
854
855
856///////////////////////////////////////////////////////////////////////////////////
857///////////////////////////////////////////////////////////////////////////////////
858TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality,
859 Int_t nbins, Float_t min, Float_t max)
860{
861 //
862 //
863 Double_t* bins = CreateLogBins(nbins, min, max);
864 TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
865 TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
866 char inputGen[1000];
867 sprintf(inputGen,"%s>>hGen", variable);
868 fTree->Draw(inputGen, selection, "groff");
869 char selectionRec[256];
870 sprintf(selectionRec, "%s && %s", selection, quality);
871 char inputRec[1000];
872 sprintf(inputRec,"%s>>hRec", variable);
873 fTree->Draw(inputRec, selectionRec, "groff");
874 //
875 TH1F* hEff = CreateEffHisto(hGen, hRec);
876 AliLabelAxes(hEff, variable, "#epsilon [%]");
877 fRes = hEff;
878 delete hRec;
879 delete hGen;
880 delete[] bins;
881 return hEff;
882}
883
884
885///////////////////////////////////////////////////////////////////////////////////
886///////////////////////////////////////////////////////////////////////////////////
887
888Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
889{
890 Double_t* bins = new Double_t[nBins+1];
891 bins[0] = xMin;
892 Double_t factor = pow(xMax/xMin, 1./nBins);
893 for (Int_t i = 1; i <= nBins; i++)
894 bins[i] = factor * bins[i-1];
895 return bins;
896}
897
898
899
900
901void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
902{
903 //
904 histo->GetXaxis()->SetTitle(xAxisTitle);
905 histo->GetYaxis()->SetTitle(yAxisTitle);
906}
907
908
909TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
910{
911 //
912 Int_t nBins = hGen->GetNbinsX();
913 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
914 hEff->SetTitle("");
915 hEff->SetStats(kFALSE);
916 hEff->SetMinimum(0.);
917 hEff->SetMaximum(110.);
918 //
919 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
920 Double_t nGen = hGen->GetBinContent(iBin);
921 Double_t nRec = hRec->GetBinContent(iBin);
922 if (nGen > 0) {
923 Double_t eff = nRec/nGen;
924 hEff->SetBinContent(iBin, 100. * eff);
925 Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
926 // if (error == 0) error = sqrt(0.1/nGen);
927 //
928 if (error == 0) error = 0.0001;
929 hEff->SetBinError(iBin, 100. * error);
930 } else {
931 hEff->SetBinContent(iBin, 100. * 0.5);
932 hEff->SetBinError(iBin, 100. * 0.5);
933 }
934 }
935 return hEff;
936}
937
938
939
940TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
941 Bool_t overflowBinFits)
942{
943 TVirtualPad* currentPad = gPad;
944 TAxis* axis = hRes2->GetXaxis();
945 Int_t nBins = axis->GetNbins();
946 TH1F* hRes, *hMean;
947 if (axis->GetXbins()->GetSize()){
948 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
949 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
950 }
951 else{
952 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
953 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
954
955 }
956 hRes->SetStats(false);
957 hRes->SetOption("E");
958 hRes->SetMinimum(0.);
959 //
960 hMean->SetStats(false);
961 hMean->SetOption("E");
962
963 // create the fit function
964 //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
965 // TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
966 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
967
968 fitFunc->SetLineWidth(2);
969 fitFunc->SetFillStyle(0);
970 // create canvas for fits
971 TCanvas* canBinFits = NULL;
972 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
973 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
974 Int_t ny = (nPads-1) / nx + 1;
975 if (drawBinFits) {
976 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
977 if (canBinFits) delete canBinFits;
978 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
979 canBinFits->Divide(nx, ny);
980 }
981
982 // loop over x bins and fit projection
983 Int_t dBin = ((overflowBinFits) ? 1 : 0);
984 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
985 if (drawBinFits) canBinFits->cd(bin + dBin);
986 TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
987 //
988 if (hBin->GetEntries() > 5) {
989 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
990 hBin->Fit(fitFunc,"s");
991 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
992
993 if (sigma > 0.){
994 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
995 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
996 }
997 else{
998 hRes->SetBinContent(bin, 0.);
999 hMean->SetBinContent(bin,0);
1000 }
1001 hRes->SetBinError(bin, fitFunc->GetParError(2));
1002 hMean->SetBinError(bin, fitFunc->GetParError(1));
1003
1004 //
1005 //
1006
1007 } else {
1008 hRes->SetBinContent(bin, 0.);
1009 hRes->SetBinError(bin, 0.);
1010 hMean->SetBinContent(bin, 0.);
1011 hMean->SetBinError(bin, 0.);
1012 }
1013
1014
1015 if (drawBinFits) {
1016 char name[256];
1017 if (bin == 0) {
1018 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1019 } else if (bin == nBins+1) {
1020 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1021 } else {
1022 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1023 axis->GetTitle(), axis->GetBinUpEdge(bin));
1024 }
1025 canBinFits->cd(bin + dBin);
1026 hBin->SetTitle(name);
1027 hBin->SetStats(kTRUE);
1028 hBin->DrawCopy("E");
1029 canBinFits->Update();
1030 canBinFits->Modified();
1031 canBinFits->Update();
1032 }
1033
1034 delete hBin;
1035 }
1036
1037 delete fitFunc;
1038 currentPad->cd();
1039 *phMean = hMean;
1040 return hRes;
1041}
1042
1043
1044