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