]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCComparisonMI.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Formatting changes.
[u/mrichter/AliRoot.git] / TPC / AliTPCComparisonMI.C
CommitLineData
e0d98ba6 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// Time Projection Chamber //
20// Comparison macro for TPC //
21// responsible:
22// marian.ivanov@cern.ch //
23//
24//
25// usage:
26
27//
28// TO CREATE OF TREE WITH MC INFO
29
30// .L AliTPCComparisonMI.C+
31// TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",1,0)
32// t->Exec()
33//
34// TO CREATE COMPARISON TREE
35// TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",1,0);
36// t2->Exec()
37//
38// EXAMPLE OF COMPARISON VISUALIZATION SESSION
39//
40// .L AliTPCComparisonMI.C+
41// TCut cprim("cprim","MC.fVDist[3]<1");
42// TCut cprim("cprim","MC.fVDist[3]<1");
43// TCut crec("crec","fReconstructed==1");
44// TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
45// TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
46// TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
47// TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
48//
49//
50// AliTPCComparisonDraw comp;
51// comp.SetIO();
52// (comp.EffVsPt("MC.fRowsWithDigits>120"+cteta1+cprim,"1"))->Draw()
53// (comp.EffVsRM("MC.fRowsWithDigits>20"+cteta1+cprim,"1"))->Draw()
54// (comp.EffVsRS("MC.fRowsWithDigits>20"+cteta1+cpos1,"1"))->Draw()
55//
56// (comp.ResPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
57// (comp.MeanPtvsPt("MC.fRowsWithDigits>20"+cteta1+cpos1,"1",0.15,2.))->Draw()
58// (comp.ResdEdxvsN("RC.fReconstructed==1&&MC.fPrim<1.5&&abs(MC.fParticle.fPdgCode)==211&&MC.fParticle.P()>0.35"+cteta1+cpos1+cprim,"1",100,160,4))->Draw()
59
60
61///////////////////////////////////////////////////////////////////////////////
62
63
64
65
7699f61f 66#if !defined(__CINT__) || defined(__MAKECINT__)
7699f61f 67#include <stdio.h>
68#include <string.h>
e0d98ba6 69//ROOT includes
7699f61f 70#include "Rtypes.h"
71#include "TFile.h"
72#include "TTree.h"
e0d98ba6 73#include "TChain.h"
74#include "TCut.h"
7699f61f 75#include "TString.h"
76#include "TBenchmark.h"
77#include "TStopwatch.h"
78#include "TParticle.h"
e0d98ba6 79#include "TSystem.h"
80#include "TTimer.h"
81#include "TVector3.h"
82#include "TPad.h"
83#include "TCanvas.h"
84#include "TH1F.h"
85#include "TH2F.h"
86#include "TF1.h"
87#include "TText.h"
88#include "Getline.h"
89#include "TStyle.h"
90
91//ALIROOT includes
7699f61f 92#include "AliRun.h"
93#include "AliStack.h"
94#include "AliTPCtrack.h"
95#include "AliSimDigits.h"
96#include "AliTPCParam.h"
7699f61f 97#include "AliTPC.h"
927d593f 98#include "AliTPCLoader.h"
7699f61f 99#include "AliDetector.h"
100#include "AliTrackReference.h"
7699f61f 101#include "AliRun.h"
7699f61f 102#include "AliTPCParamSR.h"
103#include "AliTracker.h"
104#include "AliComplexCluster.h"
927d593f 105#include "AliMagF.h"
e0d98ba6 106#endif
d8d3b5b8 107#include "AliTPCComparisonMI.h"
e0d98ba6 108
109
110
111
112
7699f61f 113
114
7699f61f 115
116//
117//
7699f61f 118
119AliTPCParam * GetTPCParam(){
120 AliTPCParamSR * par = new AliTPCParamSR;
121 par->Update();
122 return par;
123}
124
125
7699f61f 126//_____________________________________________________________________________
127Float_t TPCBetheBloch(Float_t bg)
128{
129 //
130 // Bethe-Bloch energy loss formula
131 //
132 const Double_t kp1=0.76176e-1;
133 const Double_t kp2=10.632;
134 const Double_t kp3=0.13279e-4;
135 const Double_t kp4=1.8631;
136 const Double_t kp5=1.9479;
137
138 Double_t dbg = (Double_t) bg;
139
140 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
141
142 Double_t aa = TMath::Power(beta,kp4);
143 Double_t bb = TMath::Power(1./dbg,kp5);
144
145 bb=TMath::Log(kp3+bb);
146
147 return ((Float_t)((kp2-aa-bb)*kp1/aa));
148}
149
150
151
152
153
154////////////////////////////////////////////////////////////////////////
155AliTPCGenInfo::AliTPCGenInfo()
156{
157 fReferences = 0;
e0d98ba6 158 fTRdecay.SetTrack(-1);
7699f61f 159}
160
161AliTPCGenInfo::~AliTPCGenInfo()
162{
163 if (fReferences) delete fReferences;
164 fReferences =0;
e0d98ba6 165
166}
167
168void AliTPCGenInfo::Update()
169{
170 //
171 //
172 fMCtracks =1;
173 if (!fReferences) return;
174 Float_t direction=1;
175 //Float_t rlast=0;
176 for (Int_t iref =0;iref<fReferences->GetEntriesFast();iref++){
177 AliTrackReference * ref = (AliTrackReference *) fReferences->At(iref);
178 //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
179 Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
180 if (iref==0) direction = newdirection;
181 if ( newdirection*direction<0){
182 //changed direction
183 direction = newdirection;
184 fMCtracks+=1;
185 }
186 //rlast=r;
187 }
7699f61f 188}
e0d98ba6 189
7699f61f 190////////////////////////////////////////////////////////////////////////
191
192////////////////////////////////////////////////////////////////////////
193//
194// End of implementation of the class AliTPCGenInfo
195//
196////////////////////////////////////////////////////////////////////////
197
198
199
200////////////////////////////////////////////////////////////////////////
201digitRow::digitRow()
202{
203 Reset();
204}
205////////////////////////////////////////////////////////////////////////
206digitRow & digitRow::operator=(const digitRow &digOld)
207{
208 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
209 return (*this);
210}
211////////////////////////////////////////////////////////////////////////
212void digitRow::SetRow(Int_t row)
213{
214 if (row >= 8*kgRowBytes) {
215 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
216 return;
217 }
218 Int_t iC = row/8;
219 Int_t iB = row%8;
220 SETBIT(fDig[iC],iB);
221}
222
223////////////////////////////////////////////////////////////////////////
224Bool_t digitRow::TestRow(Int_t row)
225{
226//
227// return kTRUE if row is on
228//
229 Int_t iC = row/8;
230 Int_t iB = row%8;
231 return TESTBIT(fDig[iC],iB);
232}
233////////////////////////////////////////////////////////////////////////
234Int_t digitRow::RowsOn(Int_t upto)
235{
236//
237// returns number of rows with a digit
238// count only rows less equal row number upto
239//
240 Int_t total = 0;
241 for (Int_t i = 0; i<kgRowBytes; i++) {
242 for (Int_t j = 0; j < 8; j++) {
243 if (i*8+j > upto) return total;
244 if (TESTBIT(fDig[i],j)) total++;
245 }
246 }
247 return total;
248}
249////////////////////////////////////////////////////////////////////////
250void digitRow::Reset()
251{
252//
253// resets all rows to zero
254//
255 for (Int_t i = 0; i<kgRowBytes; i++) {
256 fDig[i] <<= 8;
257 }
258}
259////////////////////////////////////////////////////////////////////////
260Int_t digitRow::Last()
261{
262//
263// returns the last row number with a digit
264// returns -1 if now digits
265//
266 for (Int_t i = kgRowBytes-1; i>=0; i--) {
267 for (Int_t j = 7; j >= 0; j--) {
268 if TESTBIT(fDig[i],j) return i*8+j;
269 }
270 }
271 return -1;
272}
273////////////////////////////////////////////////////////////////////////
274Int_t digitRow::First()
275{
276//
277// returns the first row number with a digit
278// returns -1 if now digits
279//
280 for (Int_t i = 0; i<kgRowBytes; i++) {
281 for (Int_t j = 0; j < 8; j++) {
282 if (TESTBIT(fDig[i],j)) return i*8+j;
283 }
284 }
285 return -1;
286}
287
288////////////////////////////////////////////////////////////////////////
289//
290// end of implementation of a class digitRow
291//
292////////////////////////////////////////////////////////////////////////
293
294////////////////////////////////////////////////////////////////////////
295TPCFindGenTracks::TPCFindGenTracks()
296{
297 fMCInfo = new AliTPCGenInfo;
298 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
927d593f 299
7699f61f 300 Reset();
301}
302
303////////////////////////////////////////////////////////////////////////
e0d98ba6 304TPCFindGenTracks::TPCFindGenTracks(const char * fnGalice, const char* fnRes,
7699f61f 305 Int_t nEvents, Int_t firstEvent)
306{
307 Reset();
308 fFirstEventNr = firstEvent;
309 fEventNr = firstEvent;
310 fNEvents = nEvents;
e0d98ba6 311 // fFnRes = fnRes;
312 sprintf(fFnRes,"%s",fnRes);
7699f61f 313 //
314 fLoader = AliRunLoader::Open(fnGalice);
315 if (gAlice){
316 delete gAlice->GetRunLoader();
317 delete gAlice;
318 gAlice = 0x0;
319 }
320 if (fLoader->LoadgAlice()){
321 cerr<<"Error occured while l"<<endl;
322 }
323 Int_t nall = fLoader->GetNumberOfEvents();
e0d98ba6 324 if (nEvents==0) {
325 nEvents =nall;
326 fNEvents=nall;
327 fFirstEventNr=0;
328 }
329
7699f61f 330 if (nall<=0){
331 cerr<<"no events available"<<endl;
332 fEventNr = 0;
333 return;
334 }
335 if (firstEvent+nEvents>nall) {
336 fEventNr = nall-firstEvent;
337 cerr<<"restricted number of events availaible"<<endl;
338 }
927d593f 339 AliMagF * magf = gAlice->Field();
340 AliTracker::SetFieldMap(magf);
7699f61f 341}
342
343////////////////////////////////////////////////////////////////////////
344void TPCFindGenTracks::Reset()
345{
346 fEventNr = 0;
347 fNEvents = 0;
348 fTreeGenTracks = 0;
7699f61f 349 fFileGenTracks = 0;
350 fContainerDigitRow = 0;
351 //
352 fReferences = 0;
353 fReferenceIndex0 = 0;
354 fReferenceIndex1 = 0;
e0d98ba6 355 fDecayRef = 0;
7699f61f 356 //
357 fDebug = 0;
358 fVPrim[0] = -1000.;
359 fVPrim[1] = -1000.;
360 fVPrim[2] = -1000.;
361 fParamTPC = 0;
362}
363////////////////////////////////////////////////////////////////////////
364TPCFindGenTracks::~TPCFindGenTracks()
365{
e0d98ba6 366
367 if (fLoader){
368 fLoader->UnloadgAlice();
369 gAlice = 0;
370 delete fLoader;
371 }
7699f61f 372}
373
374
375Int_t TPCFindGenTracks::SetIO()
376{
377 //
378 //
379 CreateTreeGenTracks();
380 if (!fTreeGenTracks) return 1;
927d593f 381 // AliTracker::SetFieldFactor();
382
7699f61f 383 fParamTPC = GetTPCParam();
384 //
385 return 0;
386}
387
388////////////////////////////////////////////////////////////////////////
389Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
390{
391 //
392 //
393 // SET INPUT
7699f61f 394 fLoader->SetEventNumber(eventNr);
395 //
396 fLoader->LoadHeader();
397 fLoader->LoadKinematics();
398 fStack = fLoader->Stack();
399 //
400 fLoader->LoadTrackRefs();
401 fTreeTR = fLoader->TreeTR();
402 //
403 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
404 tpcl->LoadDigits();
e0d98ba6 405 fTreeD = tpcl->TreeD();
406 return 0;
407}
408
409Int_t TPCFindGenTracks::CloseIOEvent()
410{
411 fLoader->UnloadHeader();
412 fLoader->UnloadKinematics();
413 fLoader->UnloadTrackRefs();
414 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
415 tpcl->UnloadDigits();
416 return 0;
417}
418
419Int_t TPCFindGenTracks::CloseIO()
420{
421 fLoader->UnloadgAlice();
7699f61f 422 return 0;
423}
424
e0d98ba6 425
426
7699f61f 427////////////////////////////////////////////////////////////////////////
428Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
429{
430 fNEvents = nEvents;
431 fFirstEventNr = firstEventNr;
432 return Exec();
433}
434
435////////////////////////////////////////////////////////////////////////
436Int_t TPCFindGenTracks::Exec()
437{
438 TStopwatch timer;
439 timer.Start();
440 Int_t status =SetIO();
441 if (status>0) return status;
442 //
443
444 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
445 fEventNr++) {
446 SetIO(fEventNr);
447 fNParticles = fStack->GetNtrack();
448 fContainerDigitRow = new digitRow[fNParticles];
449 //
e0d98ba6 450 fReferences = new AliTrackReference[fgMaxTR];
7699f61f 451 fReferenceIndex0 = new Int_t[fNParticles];
452 fReferenceIndex1 = new Int_t[fNParticles];
e0d98ba6 453 fDecayRef = new AliTrackReference[fNParticles];
7699f61f 454
455 for (Int_t i = 0; i<fNParticles; i++) {
456 fReferenceIndex0[i] = -1;
457 fReferenceIndex1[i] = -1;
458 }
459
460 cout<<"Start to process event "<<fEventNr<<endl;
461 cout<<"\tfNParticles = "<<fNParticles<<endl;
462 if (fDebug>2) cout<<"\tStart loop over TreeD"<<endl;
463 if (TreeDLoop()>0) return 1;
464 if (fDebug>2) cout<<"\tStart loop over TreeTR"<<endl;
465 if (TreeTRLoop()>0) return 1;
466 if (fDebug>2) cout<<"\tStart loop over TreeK"<<endl;
467 if (TreeKLoop()>0) return 1;
468 if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
469
470 delete [] fContainerDigitRow;
471 delete [] fReferences;
472 delete [] fReferenceIndex0;
473 delete [] fReferenceIndex1;
e0d98ba6 474 delete [] fDecayRef;
475 CloseIOEvent();
7699f61f 476 }
e0d98ba6 477 //
478 CloseIO();
7699f61f 479 CloseOutputFile();
480
481 cerr<<"Exec finished"<<endl;
7699f61f 482
483 timer.Stop();
484 timer.Print();
485 return 0;
486}
487////////////////////////////////////////////////////////////////////////
488void TPCFindGenTracks::CreateTreeGenTracks()
489{
490 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
491 if (!fFileGenTracks) {
492 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
493 return;
494 }
495 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
496
497
498
499 fMCInfo = new AliTPCGenInfo;
500 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
501
502 fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
503
504 fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
505
506 fTreeGenTracks->AutoSave();
507}
508////////////////////////////////////////////////////////////////////////
509void TPCFindGenTracks::CloseOutputFile()
510{
511 if (!fFileGenTracks) {
512 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
513 return;
514 }
515 fFileGenTracks->cd();
516 fTreeGenTracks->Write();
517 delete fTreeGenTracks;
518 fFileGenTracks->Close();
519 delete fFileGenTracks;
520 return;
521}
522
523////////////////////////////////////////////////////////////////////////
524Int_t TPCFindGenTracks::TreeKLoop()
525{
526//
527// open the file with treeK
528// loop over all entries there and save information about some tracks
529//
530
531 AliStack * stack = fStack;
532 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
533
534 if (fDebug > 0) {
535 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
536 <<fEventNr<<endl;
537 }
538
539 Int_t ipdg = 0;
540 TParticlePDG *ppdg = 0;
541
542
543 // not all generators give primary vertex position. Take the vertex
544 // of the particle 0 as primary vertex.
545 TDatabasePDG pdg; //get pdg table
546
547 //thank you very much root for this
548 TBranch * br = fTreeGenTracks->GetBranch("MC");
549
550 TParticle *particle = stack->ParticleFromTreeK(0);
551 fVPrim[0] = particle->Vx();
552 fVPrim[1] = particle->Vy();
553 fVPrim[2] = particle->Vz();
554 for (Int_t iParticle = 0; iParticle < fNParticles; iParticle++) {
555
556// load only particles with TR
557 if (fReferenceIndex0[iParticle]<0) continue;
558 //////////////////////////////////////////////////////////////////////
559 fMCInfo->fLabel = iParticle;
560
561 fMCInfo->fRow = (fContainerDigitRow[iParticle]);
562 fMCInfo->fRowsWithDigits = fMCInfo->fRow.RowsOn();
563 if (fMCInfo->fRowsWithDigits<10) continue;
564 fMCInfo->fRowsWithDigitsInn = fMCInfo->fRow.RowsOn(63); // 63 = number of inner rows
565 fMCInfo->fRowsTrackLength = fMCInfo->fRow.Last() - fMCInfo->fRow.First();
566 fMCInfo->fDigitsInSeed = 0;
567 if (fMCInfo->fRow.TestRow(seedRow11) && fMCInfo->fRow.TestRow(seedRow12))
568 fMCInfo->fDigitsInSeed = 1;
569 if (fMCInfo->fRow.TestRow(seedRow21) && fMCInfo->fRow.TestRow(seedRow22))
570 fMCInfo->fDigitsInSeed += 10;
571 //
572 //
573 //
574 fMCInfo->fParticle = *(stack->Particle(iParticle));
7699f61f 575 //
576 //
577 fMCInfo->fLabel = iParticle;
578 fMCInfo->fVDist[0] = fMCInfo->fParticle.Vx()-fVPrim[0];
579 fMCInfo->fVDist[1] = fMCInfo->fParticle.Vy()-fVPrim[1];
580 fMCInfo->fVDist[2] = fMCInfo->fParticle.Vz()-fVPrim[2];
581 fMCInfo->fVDist[3] = TMath::Sqrt(fMCInfo->fVDist[0]*fMCInfo->fVDist[0]+
582 fMCInfo->fVDist[1]*fMCInfo->fVDist[1]+fMCInfo->fVDist[2]*fMCInfo->fVDist[2]);
583
584 //
585 Int_t index = fReferenceIndex0[iParticle];
586 AliTrackReference ref = fReferences[index];
587 // if (ref.GetTrack()!=iParticle)
588 // printf("problem2\n");
589 fMCInfo->fTrackRef = ref;
590
591 Int_t rfindex =0;
592 if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
593 fMCInfo->fReferences =
594 new TClonesArray("AliTrackReference");
595 fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
596 //
597 for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
598 AliTrackReference ref = fReferences[i];
599 AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
e0d98ba6 600 if (ref.GetTrack()!=iParticle){
601 //printf("problem5\n");
602 continue;
603 }
7699f61f 604 *ref2 = ref;
605 rfindex++;
606 }
607 //
608 //
609 ipdg = fMCInfo->fParticle.GetPdgCode();
e0d98ba6 610 fMCInfo->fPdg = ipdg;
7699f61f 611 ppdg = pdg.GetParticle(ipdg);
612 // calculate primary ionization per cm
613 if (ppdg){
614 Float_t mass = ppdg->Mass();
615 Float_t p = TMath::Sqrt(fMCInfo->fTrackRef.Px()*fMCInfo->fTrackRef.Px()+
616 fMCInfo->fTrackRef.Py()*fMCInfo->fTrackRef.Py()+
617 fMCInfo->fTrackRef.Pz()*fMCInfo->fTrackRef.Pz());
618
619 // Float_t betagama = fMCInfo->fParticle.P()/mass;
620 Float_t betagama = p /mass;
621 fMCInfo->fPrim = TPCBetheBloch(betagama);
e0d98ba6 622 }
623 fMCInfo->fTPCdecay=kFALSE;
624 if (fDecayRef[iParticle].GetTrack()>0){
625 fMCInfo->fTRdecay = fDecayRef[iParticle];
626 fMCInfo->fDecayCoord[0] = fMCInfo->fTRdecay.X();
627 fMCInfo->fDecayCoord[1] = fMCInfo->fTRdecay.Y();
628 fMCInfo->fDecayCoord[2] = fMCInfo->fTRdecay.Z();
629 if ( (fMCInfo->fTRdecay.R()<250)&&(fMCInfo->fTRdecay.R()>85) && (TMath::Abs(fMCInfo->fTRdecay.Z())<250) ){
630 fMCInfo->fTPCdecay=kTRUE;
631 }
632 }
633 else{
634 fMCInfo->fTRdecay.SetTrack(-1);
635 fMCInfo->fDecayCoord[0] = 0;
636 fMCInfo->fDecayCoord[1] = 0;
637 fMCInfo->fDecayCoord[2] = 0;
7699f61f 638 }
e0d98ba6 639 fMCInfo->Update();
7699f61f 640 //////////////////////////////////////////////////////////////////////
e0d98ba6 641
7699f61f 642 br->SetAddress(&fMCInfo);
643 fTreeGenTracks->Fill();
644
645 }
646 fTreeGenTracks->AutoSave();
647
648 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
649
650 return 0;
651}
652
653
654
655
656////////////////////////////////////////////////////////////////////////
657Int_t TPCFindGenTracks::TreeDLoop()
658{
659//
660// open the file with treeD
661// loop over all entries there and save information about some tracks
662//
663
664 Int_t nInnerSector = fParamTPC->GetNInnerSector();
665 Int_t rowShift = 0;
666 // Int_t zero=fParamTPC->GetZeroSup();
667 Int_t zero=fParamTPC->GetZeroSup()+6;
668
669
670 //char treeDName[100];
671 //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
672 //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
673 //
674 //
675 AliSimDigits digitsAddress, *digits=&digitsAddress;
676 fTreeD->GetBranch("Segment")->SetAddress(&digits);
677
678 Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
679 if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
680 for (Int_t i=0; i<sectorsByRows; i++) {
681// for (Int_t i=5720; i<sectorsByRows; i++) {
682 if (!fTreeD->GetEvent(i)) continue;
683 Int_t sec,row;
684 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
685 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
686// cerr<<sec<<' '<<row<<endl;
687
688// here I expect that upper sectors follow lower sectors
689 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
690 //
691 //digits->ExpandTrackBuffer();
692 //Int_t *tracks = digits->GetTracks();
693 digits->First();
694 do {
695 Int_t iRow=digits->CurrentRow();
696 Int_t iColumn=digits->CurrentColumn();
697 Short_t digitValue = digits->CurrentDigit();
698// cout<<"Inner loop: sector, iRow, iColumn "
699// <<sec<<" "<<iRow<<" "<<iColumn<<endl;
700 if (digitValue >= zero) {
701 Int_t label;
702 for (Int_t j = 0; j<3; j++) {
703 label = digits->GetTrackID(iRow,iColumn,j);
704 //label = digits->GetTrackIDFast(iRow,iColumn,j);
705
706 if (label >= fNParticles) { //don't label from bakground event
707 continue;
708 }
709 if (label >= 0 && label <= fNParticles) {
710// if (label >= 0 && label <= fDebug) {
711 if (fDebug > 6 ) {
712 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
713 <<sec<<" "
714 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
715 <<" "<<row<<endl;
716 }
717 fContainerDigitRow[label].SetRow(row+rowShift);
718 }
719 }
720 }
721 } while (digits->Next());
722 }
723
724 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
725
726 return 0;
727}
728
729
730////////////////////////////////////////////////////////////////////////
731Int_t TPCFindGenTracks::TreeTRLoop()
732{
733 //
734 // loop over TrackReferences and store the first one for each track
735 //
736
737 TTree * treeTR = fTreeTR;
738 if (!treeTR) {
739 cerr<<"TreeTR not found"<<endl;
740 return 1;
741 }
742 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
743 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
744 //
745 //
e0d98ba6 746 //track references for TPC
7699f61f 747 TBranch *TPCBranchTR = treeTR->GetBranch("TPC");
748 if (!TPCBranchTR) {
749 cerr<<"TPC branch in TR not found"<<endl;
750 return 1;
751 }
752 TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
753 TPCBranchTR->SetAddress(&TPCArrayTR);
e0d98ba6 754 //get decay point if exist
755 TBranch *runbranch = treeTR->GetBranch("AliRun");
756 if (!runbranch) {
757 cerr<<"Run branch in TR not found"<<endl;
758 return 1;
759 }
760 TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
761 runbranch->SetAddress(&RunArrayTR);
7699f61f 762 //
763 //
764 //
765 Int_t index = 0;
766 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
767 TPCBranchTR->GetEntry(iPrimPart);
e0d98ba6 768 Float_t ptstart = 0;
7699f61f 769 for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
e0d98ba6 770 AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
7699f61f 771 //
7699f61f 772 Int_t label = trackRef->GetTrack();
773 if (label<0 || label > fNParticles) {
774 continue;
775 }
e0d98ba6 776 if (fReferenceIndex0[label]<0) ptstart = trackRef->Pt(); //store pt at the TPC entrance
777 if (ptstart<fgPtCut) continue;
778
7699f61f 779 if (index>=fgMaxTR) continue; //restricted size of buffer for TR
780 fReferences[index] = *trackRef;
781 fReferenceIndex1[label] = index; // the last ref with given label
782 if (fReferenceIndex0[label]==-1) fReferenceIndex0[label] = index; //the first index with label
e0d98ba6 783 index++;
7699f61f 784 }
e0d98ba6 785 // get dacay position
786 runbranch->GetEntry(iPrimPart);
787 for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
788 AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
789 //
790 if (trackRef->Pt() < fgPtCut) continue;
791 Int_t label = trackRef->GetTrack();
792 if (label<0 || label > fNParticles) {
793 continue;
794 }
795 if (trackRef->R()>450) continue; //not decay in TPC
796 if (trackRef->Z()>450) continue; //not decay in TPC
797 if (!trackRef->TestBit(BIT(2))) continue; //if not decay
798
799 if (label>=fgMaxTR) continue; //restricted size of buffer for TR
800 fDecayRef[label] = *trackRef;
801 }
802
7699f61f 803 }
e0d98ba6 804 TPCArrayTR->Delete();
805 delete TPCArrayTR;
806 RunArrayTR->Delete();
807 delete RunArrayTR;
808
7699f61f 809 return 0;
810}
811
812////////////////////////////////////////////////////////////////////////
813Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
814 AliTPCParam *paramTPC) {
815
816 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
817 Int_t index[4];
818 paramTPC->Transform0to1(x,index);
819 paramTPC->Transform1to2(x,index);
820 return x[0];
821}
822////////////////////////////////////////////////////////////////////////
823
824
825
826
827////////////////////////////////////////////////////////////////////////
828TPCCmpTr::TPCCmpTr()
829{
830 Reset();
831}
832
833////////////////////////////////////////////////////////////////////////
e0d98ba6 834TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
835 const char* fnCmp,
836 const char* fnGalice,
7699f61f 837 Int_t nEvents, Int_t firstEvent)
838{
839 Reset();
e0d98ba6 840 // fFnGenTracks = fnGenTracks;
841 // fFnCmp = fnCmp;
842 sprintf(fFnGenTracks,"%s",fnGenTracks);
843 sprintf(fFnCmp,"%s",fnCmp);
844
7699f61f 845 fFirstEventNr = firstEvent;
846 fEventNr = firstEvent;
847 fNEvents = nEvents;
848
849 //
850 fLoader = AliRunLoader::Open(fnGalice);
851 if (gAlice){
e0d98ba6 852 //delete gAlice->GetRunLoader();
7699f61f 853 delete gAlice;
854 gAlice = 0x0;
855 }
856 if (fLoader->LoadgAlice()){
857 cerr<<"Error occured while l"<<endl;
858 }
859 Int_t nall = fLoader->GetNumberOfEvents();
e0d98ba6 860 if (nEvents==0) {
861 nEvents =nall;
862 fNEvents=nall;
863 fFirstEventNr=0;
864 }
865
7699f61f 866 if (nall<=0){
867 cerr<<"no events available"<<endl;
868 fEventNr = 0;
869 return;
870 }
871 if (firstEvent+nEvents>nall) {
872 fEventNr = nall-firstEvent;
873 cerr<<"restricted number of events availaible"<<endl;
874 }
927d593f 875 AliMagF * magf = gAlice->Field();
876 AliTracker::SetFieldMap(magf);
877
7699f61f 878}
879
e0d98ba6 880
881////////////////////////////////////////////////////////////////////////
882TPCCmpTr::~TPCCmpTr()
883{
884 if (fLoader) {
885 delete fLoader;
886 }
887}
888
7699f61f 889//////////////////////////////////////////////////////////////
890Int_t TPCCmpTr::SetIO()
891{
892 //
893 //
894 CreateTreeCmp();
895 if (!fTreeCmp) return 1;
7699f61f 896 fParamTPC = GetTPCParam();
897 //
898 if (!ConnectGenTree()) {
899 cerr<<"Cannot connect tree with generated tracks"<<endl;
900 return 1;
901 }
902 return 0;
903}
904
905//////////////////////////////////////////////////////////////
906
907Int_t TPCCmpTr::SetIO(Int_t eventNr)
908{
909 //
910 //
911 // SET INPUT
e0d98ba6 912 //gAlice->GetEvent(eventNr);
7699f61f 913 fLoader->SetEventNumber(eventNr);
914 //
915 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
916 tpcl->LoadTracks();
917 fTreeRecTracks = tpcl->TreeT();
918 //
919 return 0;
920}
921
922
923
924////////////////////////////////////////////////////////////////////////
925void TPCCmpTr::Reset()
926{
927 fEventNr = 0;
928 fNEvents = 0;
929 fTreeCmp = 0;
e0d98ba6 930 // fFnCmp = "cmpTracks.root";
7699f61f 931 fFileGenTracks = 0;
932 fDebug = 0;
933 //
934 fParamTPC = 0;
935 fTreeRecTracks = 0;
936 fTreePoints =0;
937
938 fTPCTrack = 0;
939 fTracks = 0;
940 fTrackPoints =0;
941}
e0d98ba6 942
7699f61f 943////////////////////////////////////////////////////////////////////////
944Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
945{
946 fNEvents = nEvents;
947 fFirstEventNr = firstEventNr;
948 return Exec();
949}
950
951////////////////////////////////////////////////////////////////////////
952Int_t TPCCmpTr::Exec()
953{
954 TStopwatch timer;
955 timer.Start();
956
957 if (SetIO()==1)
958 return 1;
959
960 fNextTreeGenEntryToRead = 0;
961 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
962 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
963 eventNr++) {
964 SetIO(fEventNr);
965 fNParticles = gAlice->GetEvent(fEventNr);
966
967 fIndexRecTracks = new Int_t[fNParticles*4]; //write at maximum 4 tracks corresponding to particle
968 fFakeRecTracks = new Int_t[fNParticles];
969 fMultiRecTracks = new Int_t[fNParticles];
970 for (Int_t i = 0; i<fNParticles; i++) {
971 fIndexRecTracks[4*i] = -1;
972 fIndexRecTracks[4*i+1] = -1;
973 fIndexRecTracks[4*i+2] = -1;
974 fIndexRecTracks[4*i+3] = -1;
975
976 fFakeRecTracks[i] = 0;
977 fMultiRecTracks[i] = 0;
978 }
979
980 cout<<"Start to process event "<<fEventNr<<endl;
981 cout<<"\tfNParticles = "<<fNParticles<<endl;
982 if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
983 if (TreeTLoop(eventNr)>0) return 1;
984
985 if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
986 if (TreeGenLoop(eventNr)>0) return 1;
987 if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
988
989 delete fTreeRecTracks;
990 delete [] fIndexRecTracks;
991 delete [] fFakeRecTracks;
992 delete [] fMultiRecTracks;
993 }
994
995 CloseOutputFile();
996
997 cerr<<"Exec finished"<<endl;
7699f61f 998 timer.Stop();
999 timer.Print();
1000 return 0;
e0d98ba6 1001
7699f61f 1002}
1003////////////////////////////////////////////////////////////////////////
1004Bool_t TPCCmpTr::ConnectGenTree()
1005{
1006//
1007// connect all branches from the genTracksTree
1008// use the same variables as for the new cmp tree, it may work
1009//
1010 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1011 if (!fFileGenTracks) {
1012 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1013 return kFALSE;
1014 }
1015 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1016 if (!fTreeGenTracks) {
1017 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1018 <<fFnGenTracks<<endl;
1019 return kFALSE;
1020 }
1021 //
1022 fMCInfo = new AliTPCGenInfo;
1023 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1024 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1025
1026 //
1027 //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1028
1029
1030 if (fDebug > 1) {
1031 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1032 }
1033 return kTRUE;
1034}
1035
1036
1037////////////////////////////////////////////////////////////////////////
1038void TPCCmpTr::CreateTreeCmp()
1039{
1040 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1041 if (!fFileCmp) {
1042 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1043 return;
1044 }
1045
1046
1047 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1048 //
1049 fMCInfo = new AliTPCGenInfo;
1050 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1051 fRecInfo = new AliTPCRecInfo;
1052 //
1053 fTPCTrack = new AliTPCtrack;
1054 //
1055 fTreeCmp->Branch("MC","AliTPCGenInfo",&fMCInfo);
1056 fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
1057 fTreeCmp->Branch("fEventNr",&fEventNr,"fEventNr/I");
1058 fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
1059 //
1060 fTreeCmp->AutoSave();
1061
1062}
1063////////////////////////////////////////////////////////////////////////
1064void TPCCmpTr::CloseOutputFile()
1065{
1066 if (!fFileCmp) {
1067 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1068 return;
1069 }
1070 fFileCmp->cd();
1071 fTreeCmp->Write();
1072 delete fTreeCmp;
1073
1074 fFileCmp->Close();
1075 delete fFileCmp;
1076 return;
1077}
1078////////////////////////////////////////////////////////////////////////
1079
1080TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
1081 AliTPCParam *paramTPC) {
1082
1083 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1084 Int_t index[4];
1085 paramTPC->Transform0to1(x,index);
1086 paramTPC->Transform1to2(x,index);
1087 return TVector3(x);
1088}
1089////////////////////////////////////////////////////////////////////////
1090
1091Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1092{
1093 //
1094 // loop over all TPC reconstructed tracks and store info in memory
1095 //
e0d98ba6 1096 TStopwatch timer;
1097 timer.Start();
7699f61f 1098
1099 if (!fTreeRecTracks) {
1100 cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
1101 return 1;
1102 }
1103 //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
1104
1105 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1106 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1107 <<nEntries<<endl;
1108 TBranch * br= fTreeRecTracks->GetBranch("tracks");
1109 br->SetAddress(&fTPCTrack);
1110 TBranch *brp = 0;
1111 if (fTreePoints) brp = fTreePoints->GetBranch("debug");
1112
1113 if (fTracks){
1114 fTracks->Delete();
1115 delete fTracks;
1116 }
1117 if (fTrackPoints){
1118 fTrackPoints->Delete();
1119 delete fTrackPoints;
1120 fTrackPoints =0;
1121 }
1122 fTracks = new TObjArray(nEntries);
1123 if (brp){
1124 fTrackPoints = new TObjArray(nEntries);
1125 }
1126 else fTrackPoints = 0;
1127
1128 //
1129 //load tracks to the memory
1130 for (Int_t i=0; i<nEntries;i++){
1131 AliTPCtrack * track = new AliTPCtrack;
1132 br->SetAddress(&track);
1133 br->GetEntry(i);
1134 fTracks->AddAt(track,i);
1135 }
1136 //
1137 //load track points to the memory
1138 if (brp) for (Int_t i=0; i<nEntries;i++){
1139 TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
1140 brp->SetAddress(&arr);
1141 brp->GetEntry(i);
1142 if (arr!=0)
1143 for (Int_t j=0;j<arr->GetEntriesFast();j++){
1144 AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
1145 if (point && point->fID>=0){
1146 fTrackPoints->AddAt(arr,point->fID);
1147 break;
1148 }
1149 }
1150 }
1151 //
1152
1153 //
1154 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1155 //br->GetEntry(iEntry);
1156 fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
1157 //
1158 Int_t label = fTPCTrack->GetLabel();
1159 Int_t absLabel = abs(label);
1160 if (absLabel < fNParticles) {
1161 // fIndexRecTracks[absLabel] = iEntry;
1162 if (label < 0) fFakeRecTracks[absLabel]++;
1163
1164 if (fMultiRecTracks[absLabel]>0){
1165 if (fMultiRecTracks[absLabel]<4)
1166 fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
1167 }
1168 else
1169 fIndexRecTracks[absLabel*4] = iEntry;
1170 fMultiRecTracks[absLabel]++;
1171 }
1172 }
e0d98ba6 1173 printf("Time spended in TreeTLoop\n");
1174 timer.Print();
1175
7699f61f 1176 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1177
1178 return 0;
1179}
1180////////////////////////////////////////////////////////////////////////
1181Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1182{
1183//
1184// loop over all entries for a given event, find corresponding
1185// rec. track and store in the fTreeCmp
1186//
e0d98ba6 1187 TStopwatch timer;
1188 timer.Start();
7699f61f 1189 Int_t entry = fNextTreeGenEntryToRead;
1190 Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
1191 cerr<<"fNParticles, nParticlesTR, fNextTreeGenEntryToRead: "<<fNParticles<<" "
1192 <<nParticlesTR<<" "<<fNextTreeGenEntryToRead<<endl;
1193 TBranch * branch = fTreeCmp->GetBranch("RC");
e0d98ba6 1194 branch->SetAddress(&fRecInfo); // set all pointers
1195
7699f61f 1196 while (entry < nParticlesTR) {
1197 fTreeGenTracks->GetEntry(entry);
1198 entry++;
1199 if (fEventNr < eventNr) continue;
1200 if (fEventNr > eventNr) break;
1201 fNextTreeGenEntryToRead = entry-1;
1202 if (fDebug > 2 && fMCInfo->fLabel < 10) {
1203 cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
1204 }
1205
1206 fRecInfo->Reset();
1207
1208
1209 if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1210 fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1211
1212 // if (nBytes > 0) {
1213 if (fTPCTrack) {
1214 //
1215 TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1216 local.GetXYZ(fRecInfo->fTRLocalCoord);
1217 //
1218 // find nearest track if multifound
1219 if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
1220 Float_t dz = TMath::Abs(local.Z()-fTPCTrack->GetZ());
1221 for (Int_t i=1;i<4;i++){
1222 if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
1223 AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
1224 if (TMath::Abs(local.Z()-track->GetZ())<dz)
1225 fTPCTrack = track;
1226 }
1227 }
1228 }
1229 fRecInfo->fTP=0;
1230 if (fTrackPoints){
1231 Int_t id = fTPCTrack->GetUniqueID();
1232 if (fTrackPoints->UncheckedAt(id)){
1233 fRecInfo->fTP = (TClonesArray*)fTrackPoints->UncheckedAt(id);
1234 // fTrackPoints->AddAt(0,id); //not owner anymore
1235 }
1236 }
1237 fRecInfo->fTPCTrack =*fTPCTrack;
1238 fRecInfo->fReconstructed = 1;
1239 fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
1240 fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
1241 fRecInfo->fdEdx = fTPCTrack->GetdEdx();
1242 //
1243 fRecInfo->fTPCTrack.PropagateTo(local.X());
1244 Double_t par[5];
1245 //
1246 Double_t localX = local.X();
1247 fTPCTrack->GetExternalParameters(localX,par);
1248 fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
927d593f 1249 if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*TMath::Pi();
1250 if (fRecInfo->fRecPhi>=2*TMath::Pi()) fRecInfo->fRecPhi-=2*TMath::Pi();
7699f61f 1251// fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1252 fRecInfo->fLambda = TMath::ATan(par[3]);
1253 fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1254 }
1255
1256 }
1257
7699f61f 1258 fTreeCmp->Fill();
7699f61f 1259 }
1260 fTreeCmp->AutoSave();
1261 fTracks->Delete();
1262 fTPCTrack =0;
1263 if (fTrackPoints){
1264 fTrackPoints->Delete();
1265 delete fTrackPoints;
1266 fTrackPoints =0;
e0d98ba6 1267 }
1268 printf("Time spended in TreeKLoop\n");
1269 timer.Print();
7699f61f 1270 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1271
1272 return 0;
1273}
1274////////////////////////////////////////////////////////////////////////
e0d98ba6 1275////////////////////////////////////////////////////////////////////////
1276
1277void AliTPCComparisonDraw::SetIO(const char *fname)
1278{
1279 //
1280 TFile* file = TFile::Open(fname);
1281 if (!file) {
1282 printf("Could not open file - generated new one\n");
1283 TFile* filegen = TFile::Open("genTracks.root");
1284 if (!filegen){
1285 printf("FILE with MC information is generated\n");
1286 TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",0);
1287 t->Exec();
1288 delete t;
1289 }
1290 filegen = TFile::Open("genTracks.root");
1291 if (!filegen){
1292 printf("COMPARISON FILE COULDNT BE GENERATED \n");
1293 return;
1294 }
1295 printf("COMPARISON FILE IS GENERATED \n");
1296 TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",0);
1297 t2->Exec();
1298 }
1299 file = TFile::Open(fname);
1300 if (!file){
1301 printf("Comparison file couldn't be generated\n");
1302 return;
1303 }
1304 //
1305 fTree = (TTree*) file->Get("TPCcmpTracks");
1306 if (!fTree) {
1307 printf("no track comparison tree found\n");
1308 file->Close();
1309 delete file;
1310 }
1311}
1312
1313void AliTPCComparisonDraw::ResPt()
1314{
1315 //
1316 //
1317 gStyle->SetOptFit();
1318 TCanvas *c1=new TCanvas("TPC pt resolution","TPC pt resolution",0,0,700,850);
1319 c1->Draw(); c1->cd();
1320 TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99);
1321 p1->Draw();p1->cd();
1322 p1->SetGridx(); p1->SetGridy();
1323 //
1324 c1->cd();
1325 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1326 p2->Draw();p2->cd();
1327 p2->SetGridx(); p2->SetGridy();
1328 //
1329 //Default cuts
1330 TCut cprim("cprim","MC.fVDist[3]<1");
1331 TCut cnprim("cnprim","MC.fVDist[3]>1");
1332 TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1333 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1334 //
1335 c1->cd(); p1->cd();
1336 TH1F* hisp = ResPtvsPt("MC.fRowsWithDigits>100"+cteta1+cpos1,"1",0.15,2.,6);
1337 c1->cd();
1338 c1->Draw();
1339 p1->cd();
1340 p1->Draw();
1341 hisp->Draw();
1342 //
1343 c1->cd();
1344 c1->Draw();
1345 p2->cd();
1346 p2->Draw();
1347 TH1F* his2 = new TH1F("Ptresolution","Ptresolution",40,-5.,5.);
1348 fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt()>>Ptresolution","MC.fRowsWithDigits>100&&RC.fTPCTrack.fN>50&&RC.fMultiple==1"+cteta1+cpos1+cprim);
1349 AliLabelAxes(his2, "#Delta p_{t} / p_{t} [%]", "entries");
1350 his2->Fit("gaus");
1351 his2->Draw();
1352
1353}
1354
1355void AliTPCComparisonDraw::Eff()
1356{
1357 //
1358 //
1359 TCanvas *c1=new TCanvas("TPC efficiency","TPC efficiency",0,0,700,850);
1360 c1->Draw(); c1->cd();
1361 TPad *p1=new TPad("p1","p1",0.01,0.51,.99,.99);
1362 p1->Draw();p1->cd();
1363 p1->SetGridx(); p1->SetGridy();
1364 //
1365 c1->cd();
1366 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1367 p2->Draw();p2->cd();
1368 p2->SetGridx(); p2->SetGridy();
1369 //
1370 //Default cuts
1371 TCut cprim("cprim","MC.fVDist[3]<1");
1372 TCut cnprim("cnprim","MC.fVDist[3]>1");
1373 TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
1374 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
1375 //
1376 c1->cd(); p1->cd();
1377 TH1F* hisp = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cprim,"RC.fTPCTrack.fN>50");
1378 hisp->Draw();
1379 //hisp->DrawClone();
1380 TText * text = new TText(0.25,102.,"Primary particles");
1381 text->SetTextSize(0.05);
1382 text->Draw();
1383 //
1384 c1->cd(); p2->cd();
1385 TH1F* hiss = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cnprim,"RC.fTPCTrack.fN>50");
1386 hiss->Draw();
1387 //hiss->DrawClone();
1388 text = new TText(0.25,102.,"Secondary particles");
1389 text->SetTextSize(0.05);
1390 text->Draw();
1391
1392
1393}
1394
1395
1396TH1F * AliTPCComparisonDraw::EffVsPt(const char* selection, const char * quality, Float_t min, Float_t max)
1397{
1398 //
1399 //
1400 Int_t nBins = 10;
1401 Double_t* bins = CreateLogBins(nBins, min, max);
1402 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, bins);
1403 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, bins);
1404
1405 fTree->Draw("MC.fParticle.Pt()>>hGen", selection, "groff");
1406 char selectionRec[256];
1407 sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1408 fTree->Draw("MC.fParticle.Pt()>>hRec", selectionRec, "groff");
1409
1410 TH1F* hEff = CreateEffHisto(hGen, hRec);
1411 AliLabelAxes(hEff, "p_{t} [GeV/c]", "#epsilon [%]");
1412
1413 delete hRec;
1414 delete hGen;
1415 delete[] bins;
1416 return hEff;
1417}
1418
1419
1420TH1F * AliTPCComparisonDraw::EffVsRM(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1421{
1422 //
1423 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1424 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1425 //
1426 fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hGen", selection, "groff");
1427 char selectionRec[256];
1428 sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1429 fTree->Draw("sqrt(MC.fDecayCoord[0]*MC.fDecayCoord[0] + MC.fDecayCoord[1]*MC.fDecayCoord[1])>>hRec", selectionRec, "groff");
1430 //
1431 TH1F* hEff = CreateEffHisto(hGen, hRec);
1432 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1433 //
1434 delete hRec;
1435 delete hGen;
1436 return hEff;
1437}
1438
1439TH1F * AliTPCComparisonDraw::EffVsRS(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1440{
1441 //
1442 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1443 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
1444 //
1445 fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hGen", selection, "groff");
1446 char selectionRec[256];
1447 sprintf(selectionRec, "%s && RC.fReconstructed && %s", selection, quality);
1448 fTree->Draw("sqrt(MC.fVDist[0]*MC.fVDist[0] + MC.fVDist[1]*MC.fVDist[1])>>hRec", selectionRec, "groff");
1449 //
1450 TH1F* hEff = CreateEffHisto(hGen, hRec);
1451 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1452 //
1453 delete hRec;
1454 delete hGen;
1455 return hEff;
1456
1457}
1458
1459TH1F * AliTPCComparisonDraw::ResPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1460{
1461 //
1462 Double_t* bins = CreateLogBins(nBins, min, max);
1463 Int_t nBinsRes = 30;
1464 Double_t maxRes = 10.;
1465 TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1466
1467 fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1468
1469 TH1F* hMean=0;
1470 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1471 AliLabelAxes(hRes, "p_{t} [GeV/c]", "#Delta p_{t} / p_{t} [%]");
1472 //
1473 delete hRes2;
1474 delete[] bins;
1475 return hRes;
1476
1477}
1478
1479TH1F * AliTPCComparisonDraw::MeanPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1480{
1481 //
1482 Double_t* bins = CreateLogBins(nBins, min, max);
1483 Int_t nBinsRes = 30;
1484 Double_t maxRes = 10.;
1485 TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, bins, nBinsRes, -maxRes, maxRes);
1486
1487 fTree->Draw("100.*(1./fTPCTrack.Get1Pt()-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1488
1489 TH1F* hMean=0;
1490 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1491 AliLabelAxes(hRes, "p_{t} [GeV/c]", "mean p_{t} / p_{t} [%]");
1492 //
1493 delete hRes2;
1494 delete[] bins;
1495 if (!hMean) return 0;
1496 return hMean;
1497
1498}
1499
1500
1501TH1F * AliTPCComparisonDraw::ResdEdxvsN(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1502{
1503 //
1504 Int_t nBinsRes = 15;
1505
1506 TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, min, max, nBinsRes, 34, 60);
1507
1508 fTree->Draw("RC.fTPCTrack.fdEdx/MC.fPrim:RC.fTPCTrack.fN>>hRes2", selection, "groff");
1509
1510 TH1F* hMean=0;
1511 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1512 AliLabelAxes(hRes, "N points", "sigma dEdx/Nprim [%]");
1513 //
1514 delete hRes2;
1515 return hRes;
1516
1517}
1518
1519
1520
1521Double_t* AliTPCComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1522{
1523 Double_t* bins = new Double_t[nBins+1];
1524 bins[0] = xMin;
1525 Double_t factor = pow(xMax/xMin, 1./nBins);
1526 for (Int_t i = 1; i <= nBins; i++)
1527 bins[i] = factor * bins[i-1];
1528 return bins;
1529}
1530
1531
1532
1533
1534void AliTPCComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1535{
1536 //
1537 histo->GetXaxis()->SetTitle(xAxisTitle);
1538 histo->GetYaxis()->SetTitle(yAxisTitle);
1539}
1540
1541
1542TH1F* AliTPCComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1543{
1544 //
1545 Int_t nBins = hGen->GetNbinsX();
1546 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1547 hEff->SetTitle("");
1548 hEff->SetStats(kFALSE);
1549 hEff->SetMinimum(0.);
1550 hEff->SetMaximum(110.);
1551 //
1552 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1553 Double_t nGen = hGen->GetBinContent(iBin);
1554 Double_t nRec = hRec->GetBinContent(iBin);
1555 if (nGen > 0) {
1556 Double_t eff = nRec/nGen;
1557 hEff->SetBinContent(iBin, 100. * eff);
1558 Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
1559 // if (error == 0) error = sqrt(0.1/nGen);
1560 //
1561 if (error == 0) error = 0.0001;
1562 hEff->SetBinError(iBin, 100. * error);
1563 } else {
1564 hEff->SetBinContent(iBin, 100. * 0.5);
1565 hEff->SetBinError(iBin, 100. * 0.5);
1566 }
1567 }
1568 return hEff;
1569}
1570
1571
1572
1573TH1F* AliTPCComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
1574 Bool_t overflowBinFits)
1575{
1576 TVirtualPad* currentPad = gPad;
1577 TAxis* axis = hRes2->GetXaxis();
1578 Int_t nBins = axis->GetNbins();
1579 TH1F* hRes, *hMean;
1580 if (axis->GetXbins()->GetSize()){
1581 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1582 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1583 }
1584 else{
1585 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1586 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1587
1588 }
1589 hRes->SetStats(false);
1590 hRes->SetOption("E");
1591 hRes->SetMinimum(0.);
1592 //
1593 hMean->SetStats(false);
1594 hMean->SetOption("E");
1595
1596 // create the fit function
1597 //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1598 // TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1599 TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1600
1601 fitFunc->SetLineWidth(2);
1602 fitFunc->SetFillStyle(0);
1603 // create canvas for fits
1604 TCanvas* canBinFits = NULL;
1605 Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1606 Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1607 Int_t ny = (nPads-1) / nx + 1;
1608 if (drawBinFits) {
1609 canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1610 if (canBinFits) delete canBinFits;
1611 canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1612 canBinFits->Divide(nx, ny);
1613 }
1614
1615 // loop over x bins and fit projection
1616 Int_t dBin = ((overflowBinFits) ? 1 : 0);
1617 for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1618 if (drawBinFits) canBinFits->cd(bin + dBin);
1619 TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1620 //
982aff31 1621 if (hBin->GetEntries() > 5) {
e0d98ba6 1622 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1623 hBin->Fit(fitFunc,"s");
1624 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1625
1626 if (sigma > 0.){
1627 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1628 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1629 }
1630 else{
1631 hRes->SetBinContent(bin, 0.);
1632 hMean->SetBinContent(bin,0);
1633 }
1634 hRes->SetBinError(bin, fitFunc->GetParError(2));
1635 hMean->SetBinError(bin, fitFunc->GetParError(1));
1636
1637 //
1638 //
1639
1640 } else {
1641 hRes->SetBinContent(bin, 0.);
1642 hRes->SetBinError(bin, 0.);
1643 hMean->SetBinContent(bin, 0.);
1644 hMean->SetBinError(bin, 0.);
1645 }
1646
1647
1648 if (drawBinFits) {
1649 char name[256];
1650 if (bin == 0) {
1651 sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1652 } else if (bin == nBins+1) {
1653 sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1654 } else {
1655 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1656 axis->GetTitle(), axis->GetBinUpEdge(bin));
1657 }
1658 canBinFits->cd(bin + dBin);
1659 hBin->SetTitle(name);
1660 hBin->SetStats(kTRUE);
1661 hBin->DrawCopy("E");
1662 canBinFits->Update();
1663 canBinFits->Modified();
1664 canBinFits->Update();
1665 }
1666
1667 delete hBin;
1668 }
1669
1670 delete fitFunc;
1671 currentPad->cd();
1672 *phMean = hMean;
1673 return hRes;
1674}
1675
1676
1677