1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Time Projection Chamber //
20 // Comparison macro for TPC //
22 // marian.ivanov@cern.ch //
28 // TO CREATE OF TREE WITH MC INFO
30 // .L AliTPCComparisonMI.C+
31 // TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",1,0)
34 // TO CREATE COMPARISON TREE
35 // TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",1,0);
38 // EXAMPLE OF COMPARISON VISUALIZATION SESSION
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)");
50 // AliTPCComparisonDraw comp;
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()
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()
61 ///////////////////////////////////////////////////////////////////////////////
66 #if !defined(__CINT__) || defined(__MAKECINT__)
76 #include "TBenchmark.h"
77 #include "TStopwatch.h"
78 #include "TParticle.h"
94 #include "AliTPCtrack.h"
95 #include "AliSimDigits.h"
96 #include "AliTPCParam.h"
98 #include "AliTPCLoader.h"
99 #include "AliDetector.h"
100 #include "AliTrackReference.h"
102 #include "AliTPCParamSR.h"
103 #include "AliTracker.h"
104 #include "AliComplexCluster.h"
107 #include "AliTPCComparisonMI.h"
119 AliTPCParam * GetTPCParam(){
120 AliTPCParamSR * par = new AliTPCParamSR;
126 //_____________________________________________________________________________
127 Float_t TPCBetheBloch(Float_t bg)
130 // Bethe-Bloch energy loss formula
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;
138 Double_t dbg = (Double_t) bg;
140 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
142 Double_t aa = TMath::Power(beta,kp4);
143 Double_t bb = TMath::Power(1./dbg,kp5);
145 bb=TMath::Log(kp3+bb);
147 return ((Float_t)((kp2-aa-bb)*kp1/aa));
154 ////////////////////////////////////////////////////////////////////////
155 AliTPCGenInfo::AliTPCGenInfo()
158 fTRdecay.SetTrack(-1);
161 AliTPCGenInfo::~AliTPCGenInfo()
163 if (fReferences) delete fReferences;
168 void AliTPCGenInfo::Update()
173 if (!fReferences) return;
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){
183 direction = newdirection;
190 ////////////////////////////////////////////////////////////////////////
192 ////////////////////////////////////////////////////////////////////////
194 // End of implementation of the class AliTPCGenInfo
196 ////////////////////////////////////////////////////////////////////////
200 ////////////////////////////////////////////////////////////////////////
205 ////////////////////////////////////////////////////////////////////////
206 digitRow & digitRow::operator=(const digitRow &digOld)
208 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
211 ////////////////////////////////////////////////////////////////////////
212 void digitRow::SetRow(Int_t row)
214 if (row >= 8*kgRowBytes) {
215 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
223 ////////////////////////////////////////////////////////////////////////
224 Bool_t digitRow::TestRow(Int_t row)
227 // return kTRUE if row is on
231 return TESTBIT(fDig[iC],iB);
233 ////////////////////////////////////////////////////////////////////////
234 Int_t digitRow::RowsOn(Int_t upto)
237 // returns number of rows with a digit
238 // count only rows less equal row number upto
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++;
249 ////////////////////////////////////////////////////////////////////////
250 void digitRow::Reset()
253 // resets all rows to zero
255 for (Int_t i = 0; i<kgRowBytes; i++) {
259 ////////////////////////////////////////////////////////////////////////
260 Int_t digitRow::Last()
263 // returns the last row number with a digit
264 // returns -1 if now digits
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;
273 ////////////////////////////////////////////////////////////////////////
274 Int_t digitRow::First()
277 // returns the first row number with a digit
278 // returns -1 if now digits
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;
288 ////////////////////////////////////////////////////////////////////////
290 // end of implementation of a class digitRow
292 ////////////////////////////////////////////////////////////////////////
294 ////////////////////////////////////////////////////////////////////////
295 TPCFindGenTracks::TPCFindGenTracks()
297 fMCInfo = new AliTPCGenInfo;
298 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
303 ////////////////////////////////////////////////////////////////////////
304 TPCFindGenTracks::TPCFindGenTracks(const char * fnGalice, const char* fnRes,
305 Int_t nEvents, Int_t firstEvent)
308 fFirstEventNr = firstEvent;
309 fEventNr = firstEvent;
312 sprintf(fFnRes,"%s",fnRes);
314 fLoader = AliRunLoader::Open(fnGalice);
316 delete gAlice->GetRunLoader();
320 if (fLoader->LoadgAlice()){
321 cerr<<"Error occured while l"<<endl;
323 Int_t nall = fLoader->GetNumberOfEvents();
331 cerr<<"no events available"<<endl;
335 if (firstEvent+nEvents>nall) {
336 fEventNr = nall-firstEvent;
337 cerr<<"restricted number of events availaible"<<endl;
339 AliMagF * magf = gAlice->Field();
340 AliTracker::SetFieldMap(magf);
343 ////////////////////////////////////////////////////////////////////////
344 void TPCFindGenTracks::Reset()
350 fContainerDigitRow = 0;
353 fReferenceIndex0 = 0;
354 fReferenceIndex1 = 0;
363 ////////////////////////////////////////////////////////////////////////
364 TPCFindGenTracks::~TPCFindGenTracks()
368 fLoader->UnloadgAlice();
375 Int_t TPCFindGenTracks::SetIO()
379 CreateTreeGenTracks();
380 if (!fTreeGenTracks) return 1;
381 // AliTracker::SetFieldFactor();
383 fParamTPC = GetTPCParam();
388 ////////////////////////////////////////////////////////////////////////
389 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
394 fLoader->SetEventNumber(eventNr);
396 fLoader->LoadHeader();
397 fLoader->LoadKinematics();
398 fStack = fLoader->Stack();
400 fLoader->LoadTrackRefs();
401 fTreeTR = fLoader->TreeTR();
403 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
405 fTreeD = tpcl->TreeD();
409 Int_t TPCFindGenTracks::CloseIOEvent()
411 fLoader->UnloadHeader();
412 fLoader->UnloadKinematics();
413 fLoader->UnloadTrackRefs();
414 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
415 tpcl->UnloadDigits();
419 Int_t TPCFindGenTracks::CloseIO()
421 fLoader->UnloadgAlice();
427 ////////////////////////////////////////////////////////////////////////
428 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
431 fFirstEventNr = firstEventNr;
435 ////////////////////////////////////////////////////////////////////////
436 Int_t TPCFindGenTracks::Exec()
440 Int_t status =SetIO();
441 if (status>0) return status;
444 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
447 fNParticles = fStack->GetNtrack();
448 fContainerDigitRow = new digitRow[fNParticles];
450 fReferences = new AliTrackReference[fgMaxTR];
451 fReferenceIndex0 = new Int_t[fNParticles];
452 fReferenceIndex1 = new Int_t[fNParticles];
453 fDecayRef = new AliTrackReference[fNParticles];
455 for (Int_t i = 0; i<fNParticles; i++) {
456 fReferenceIndex0[i] = -1;
457 fReferenceIndex1[i] = -1;
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;
470 delete [] fContainerDigitRow;
471 delete [] fReferences;
472 delete [] fReferenceIndex0;
473 delete [] fReferenceIndex1;
481 cerr<<"Exec finished"<<endl;
487 ////////////////////////////////////////////////////////////////////////
488 void TPCFindGenTracks::CreateTreeGenTracks()
490 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
491 if (!fFileGenTracks) {
492 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
495 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
499 fMCInfo = new AliTPCGenInfo;
500 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
502 fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
504 fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
506 fTreeGenTracks->AutoSave();
508 ////////////////////////////////////////////////////////////////////////
509 void TPCFindGenTracks::CloseOutputFile()
511 if (!fFileGenTracks) {
512 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
515 fFileGenTracks->cd();
516 fTreeGenTracks->Write();
517 delete fTreeGenTracks;
518 fFileGenTracks->Close();
519 delete fFileGenTracks;
523 ////////////////////////////////////////////////////////////////////////
524 Int_t TPCFindGenTracks::TreeKLoop()
527 // open the file with treeK
528 // loop over all entries there and save information about some tracks
531 AliStack * stack = fStack;
532 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
535 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
540 TParticlePDG *ppdg = 0;
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
547 //thank you very much root for this
548 TBranch * br = fTreeGenTracks->GetBranch("MC");
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++) {
556 // load only particles with TR
557 if (fReferenceIndex0[iParticle]<0) continue;
558 //////////////////////////////////////////////////////////////////////
559 fMCInfo->fLabel = iParticle;
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;
574 fMCInfo->fParticle = *(stack->Particle(iParticle));
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]);
585 Int_t index = fReferenceIndex0[iParticle];
586 AliTrackReference ref = fReferences[index];
587 // if (ref.GetTrack()!=iParticle)
588 // printf("problem2\n");
589 fMCInfo->fTrackRef = ref;
592 if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
593 fMCInfo->fReferences =
594 new TClonesArray("AliTrackReference");
595 fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
597 for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
598 AliTrackReference ref = fReferences[i];
599 AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
600 if (ref.GetTrack()!=iParticle){
601 //printf("problem5\n");
609 ipdg = fMCInfo->fParticle.GetPdgCode();
610 fMCInfo->fPdg = ipdg;
611 ppdg = pdg.GetParticle(ipdg);
612 // calculate primary ionization per cm
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());
619 // Float_t betagama = fMCInfo->fParticle.P()/mass;
620 Float_t betagama = p /mass;
621 fMCInfo->fPrim = TPCBetheBloch(betagama);
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;
634 fMCInfo->fTRdecay.SetTrack(-1);
635 fMCInfo->fDecayCoord[0] = 0;
636 fMCInfo->fDecayCoord[1] = 0;
637 fMCInfo->fDecayCoord[2] = 0;
640 //////////////////////////////////////////////////////////////////////
642 br->SetAddress(&fMCInfo);
643 fTreeGenTracks->Fill();
646 fTreeGenTracks->AutoSave();
648 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
656 ////////////////////////////////////////////////////////////////////////
657 Int_t TPCFindGenTracks::TreeDLoop()
660 // open the file with treeD
661 // loop over all entries there and save information about some tracks
664 Int_t nInnerSector = fParamTPC->GetNInnerSector();
666 // Int_t zero=fParamTPC->GetZeroSup();
667 Int_t zero=fParamTPC->GetZeroSup()+6;
670 //char treeDName[100];
671 //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
672 //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
675 AliSimDigits digitsAddress, *digits=&digitsAddress;
676 fTreeD->GetBranch("Segment")->SetAddress(&digits);
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;
684 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
685 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
686 // cerr<<sec<<' '<<row<<endl;
688 // here I expect that upper sectors follow lower sectors
689 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
691 //digits->ExpandTrackBuffer();
692 //Int_t *tracks = digits->GetTracks();
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) {
702 for (Int_t j = 0; j<3; j++) {
703 label = digits->GetTrackID(iRow,iColumn,j);
704 //label = digits->GetTrackIDFast(iRow,iColumn,j);
706 if (label >= fNParticles) { //don't label from bakground event
709 if (label >= 0 && label <= fNParticles) {
710 // if (label >= 0 && label <= fDebug) {
712 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
714 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
717 fContainerDigitRow[label].SetRow(row+rowShift);
721 } while (digits->Next());
724 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
730 ////////////////////////////////////////////////////////////////////////
731 Int_t TPCFindGenTracks::TreeTRLoop()
734 // loop over TrackReferences and store the first one for each track
737 TTree * treeTR = fTreeTR;
739 cerr<<"TreeTR not found"<<endl;
742 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
743 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
746 //track references for TPC
747 TBranch *TPCBranchTR = treeTR->GetBranch("TPC");
749 cerr<<"TPC branch in TR not found"<<endl;
752 TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
753 TPCBranchTR->SetAddress(&TPCArrayTR);
754 //get decay point if exist
755 TBranch *runbranch = treeTR->GetBranch("AliRun");
757 cerr<<"Run branch in TR not found"<<endl;
760 TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
761 runbranch->SetAddress(&RunArrayTR);
766 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
767 TPCBranchTR->GetEntry(iPrimPart);
769 for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
770 AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
772 Int_t label = trackRef->GetTrack();
773 if (label<0 || label > fNParticles) {
776 if (fReferenceIndex0[label]<0) ptstart = trackRef->Pt(); //store pt at the TPC entrance
777 if (ptstart<fgPtCut) continue;
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
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);
790 if (trackRef->Pt() < fgPtCut) continue;
791 Int_t label = trackRef->GetTrack();
792 if (label<0 || label > fNParticles) {
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
799 if (label>=fgMaxTR) continue; //restricted size of buffer for TR
800 fDecayRef[label] = *trackRef;
804 TPCArrayTR->Delete();
806 RunArrayTR->Delete();
812 ////////////////////////////////////////////////////////////////////////
813 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
814 AliTPCParam *paramTPC) {
816 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
818 paramTPC->Transform0to1(x,index);
819 paramTPC->Transform1to2(x,index);
822 ////////////////////////////////////////////////////////////////////////
827 ////////////////////////////////////////////////////////////////////////
833 ////////////////////////////////////////////////////////////////////////
834 TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
836 const char* fnGalice,
837 Int_t nEvents, Int_t firstEvent)
840 // fFnGenTracks = fnGenTracks;
842 sprintf(fFnGenTracks,"%s",fnGenTracks);
843 sprintf(fFnCmp,"%s",fnCmp);
845 fFirstEventNr = firstEvent;
846 fEventNr = firstEvent;
850 fLoader = AliRunLoader::Open(fnGalice);
852 //delete gAlice->GetRunLoader();
856 if (fLoader->LoadgAlice()){
857 cerr<<"Error occured while l"<<endl;
859 Int_t nall = fLoader->GetNumberOfEvents();
867 cerr<<"no events available"<<endl;
871 if (firstEvent+nEvents>nall) {
872 fEventNr = nall-firstEvent;
873 cerr<<"restricted number of events availaible"<<endl;
875 AliMagF * magf = gAlice->Field();
876 AliTracker::SetFieldMap(magf);
881 ////////////////////////////////////////////////////////////////////////
882 TPCCmpTr::~TPCCmpTr()
889 //////////////////////////////////////////////////////////////
890 Int_t TPCCmpTr::SetIO()
895 if (!fTreeCmp) return 1;
896 fParamTPC = GetTPCParam();
898 if (!ConnectGenTree()) {
899 cerr<<"Cannot connect tree with generated tracks"<<endl;
905 //////////////////////////////////////////////////////////////
907 Int_t TPCCmpTr::SetIO(Int_t eventNr)
912 //gAlice->GetEvent(eventNr);
913 fLoader->SetEventNumber(eventNr);
915 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
917 fTreeRecTracks = tpcl->TreeT();
924 ////////////////////////////////////////////////////////////////////////
925 void TPCCmpTr::Reset()
930 // fFnCmp = "cmpTracks.root";
943 ////////////////////////////////////////////////////////////////////////
944 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
947 fFirstEventNr = firstEventNr;
951 ////////////////////////////////////////////////////////////////////////
952 Int_t TPCCmpTr::Exec()
960 fNextTreeGenEntryToRead = 0;
961 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
962 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
965 fNParticles = gAlice->GetEvent(fEventNr);
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;
976 fFakeRecTracks[i] = 0;
977 fMultiRecTracks[i] = 0;
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;
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;
989 delete fTreeRecTracks;
990 delete [] fIndexRecTracks;
991 delete [] fFakeRecTracks;
992 delete [] fMultiRecTracks;
997 cerr<<"Exec finished"<<endl;
1003 ////////////////////////////////////////////////////////////////////////
1004 Bool_t TPCCmpTr::ConnectGenTree()
1007 // connect all branches from the genTracksTree
1008 // use the same variables as for the new cmp tree, it may work
1010 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1011 if (!fFileGenTracks) {
1012 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1015 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1016 if (!fTreeGenTracks) {
1017 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1018 <<fFnGenTracks<<endl;
1022 fMCInfo = new AliTPCGenInfo;
1023 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1024 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1027 //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1031 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1037 ////////////////////////////////////////////////////////////////////////
1038 void TPCCmpTr::CreateTreeCmp()
1040 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1042 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1047 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1049 fMCInfo = new AliTPCGenInfo;
1050 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1051 fRecInfo = new AliTPCRecInfo;
1053 fTPCTrack = new AliTPCtrack;
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);
1060 fTreeCmp->AutoSave();
1063 ////////////////////////////////////////////////////////////////////////
1064 void TPCCmpTr::CloseOutputFile()
1067 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1078 ////////////////////////////////////////////////////////////////////////
1080 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
1081 AliTPCParam *paramTPC) {
1083 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1085 paramTPC->Transform0to1(x,index);
1086 paramTPC->Transform1to2(x,index);
1089 ////////////////////////////////////////////////////////////////////////
1091 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1094 // loop over all TPC reconstructed tracks and store info in memory
1099 if (!fTreeRecTracks) {
1100 cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
1103 //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
1105 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1106 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1108 TBranch * br= fTreeRecTracks->GetBranch("tracks");
1109 br->SetAddress(&fTPCTrack);
1111 if (fTreePoints) brp = fTreePoints->GetBranch("debug");
1118 fTrackPoints->Delete();
1119 delete fTrackPoints;
1122 fTracks = new TObjArray(nEntries);
1124 fTrackPoints = new TObjArray(nEntries);
1126 else fTrackPoints = 0;
1129 //load tracks to the memory
1130 for (Int_t i=0; i<nEntries;i++){
1131 AliTPCtrack * track = new AliTPCtrack;
1132 br->SetAddress(&track);
1134 fTracks->AddAt(track,i);
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);
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);
1154 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1155 //br->GetEntry(iEntry);
1156 fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
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]++;
1164 if (fMultiRecTracks[absLabel]>0){
1165 if (fMultiRecTracks[absLabel]<4)
1166 fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
1169 fIndexRecTracks[absLabel*4] = iEntry;
1170 fMultiRecTracks[absLabel]++;
1173 printf("Time spended in TreeTLoop\n");
1176 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1180 ////////////////////////////////////////////////////////////////////////
1181 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1184 // loop over all entries for a given event, find corresponding
1185 // rec. track and store in the fTreeCmp
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");
1194 branch->SetAddress(&fRecInfo); // set all pointers
1196 while (entry < nParticlesTR) {
1197 fTreeGenTracks->GetEntry(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;
1209 if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1210 fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1212 // if (nBytes > 0) {
1215 TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1216 local.GetXYZ(fRecInfo->fTRLocalCoord);
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)
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
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();
1243 fRecInfo->fTPCTrack.PropagateTo(local.X());
1246 Double_t localX = local.X();
1247 fTPCTrack->GetExternalParameters(localX,par);
1248 fRecInfo->fRecPhi=TMath::ASin(par[2]) + fTPCTrack->GetAlpha();
1249 if (fRecInfo->fRecPhi<0) fRecInfo->fRecPhi+=2*TMath::Pi();
1250 if (fRecInfo->fRecPhi>=2*TMath::Pi()) fRecInfo->fRecPhi-=2*TMath::Pi();
1251 // fRecInfo->fRecPhi = (fRecInfo->fRecPhi)*kRaddeg;
1252 fRecInfo->fLambda = TMath::ATan(par[3]);
1253 fRecInfo->fRecPt_1 = TMath::Abs(par[4]);
1260 fTreeCmp->AutoSave();
1264 fTrackPoints->Delete();
1265 delete fTrackPoints;
1268 printf("Time spended in TreeKLoop\n");
1270 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1274 ////////////////////////////////////////////////////////////////////////
1275 ////////////////////////////////////////////////////////////////////////
1277 void AliTPCComparisonDraw::SetIO(const char *fname)
1280 TFile* file = TFile::Open(fname);
1282 printf("Could not open file - generated new one\n");
1283 TFile* filegen = TFile::Open("genTracks.root");
1285 printf("FILE with MC information is generated\n");
1286 TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",0);
1290 filegen = TFile::Open("genTracks.root");
1292 printf("COMPARISON FILE COULDNT BE GENERATED \n");
1295 printf("COMPARISON FILE IS GENERATED \n");
1296 TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",0);
1299 file = TFile::Open(fname);
1301 printf("Comparison file couldn't be generated\n");
1305 fTree = (TTree*) file->Get("TPCcmpTracks");
1307 printf("no track comparison tree found\n");
1313 void AliTPCComparisonDraw::ResPt()
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();
1325 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1326 p2->Draw();p2->cd();
1327 p2->SetGridx(); p2->SetGridy();
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");
1336 TH1F* hisp = ResPtvsPt("MC.fRowsWithDigits>100"+cteta1+cpos1,"1",0.15,2.,6);
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");
1355 void AliTPCComparisonDraw::Eff()
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();
1366 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1367 p2->Draw();p2->cd();
1368 p2->SetGridx(); p2->SetGridy();
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");
1377 TH1F* hisp = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cprim,"RC.fTPCTrack.fN>50");
1379 //hisp->DrawClone();
1380 TText * text = new TText(0.25,102.,"Primary particles");
1381 text->SetTextSize(0.05);
1385 TH1F* hiss = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cnprim,"RC.fTPCTrack.fN>50");
1387 //hiss->DrawClone();
1388 text = new TText(0.25,102.,"Secondary particles");
1389 text->SetTextSize(0.05);
1396 TH1F * AliTPCComparisonDraw::EffVsPt(const char* selection, const char * quality, Float_t min, Float_t max)
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);
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");
1410 TH1F* hEff = CreateEffHisto(hGen, hRec);
1411 AliLabelAxes(hEff, "p_{t} [GeV/c]", "#epsilon [%]");
1420 TH1F * AliTPCComparisonDraw::EffVsRM(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1423 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1424 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
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");
1431 TH1F* hEff = CreateEffHisto(hGen, hRec);
1432 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1439 TH1F * AliTPCComparisonDraw::EffVsRS(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1442 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1443 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
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");
1450 TH1F* hEff = CreateEffHisto(hGen, hRec);
1451 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1459 TH1F * AliTPCComparisonDraw::ResPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
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);
1467 fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1470 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1471 AliLabelAxes(hRes, "p_{t} [GeV/c]", "#Delta p_{t} / p_{t} [%]");
1479 TH1F * AliTPCComparisonDraw::MeanPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
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);
1487 fTree->Draw("100.*(1./fTPCTrack.Get1Pt()-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1490 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1491 AliLabelAxes(hRes, "p_{t} [GeV/c]", "mean p_{t} / p_{t} [%]");
1495 if (!hMean) return 0;
1501 TH1F * AliTPCComparisonDraw::ResdEdxvsN(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1504 Int_t nBinsRes = 15;
1506 TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, min, max, nBinsRes, 34, 60);
1508 fTree->Draw("RC.fTPCTrack.fdEdx/MC.fPrim:RC.fTPCTrack.fN>>hRes2", selection, "groff");
1511 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1512 AliLabelAxes(hRes, "N points", "sigma dEdx/Nprim [%]");
1521 Double_t* AliTPCComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1523 Double_t* bins = new Double_t[nBins+1];
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];
1534 void AliTPCComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1537 histo->GetXaxis()->SetTitle(xAxisTitle);
1538 histo->GetYaxis()->SetTitle(yAxisTitle);
1542 TH1F* AliTPCComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1545 Int_t nBins = hGen->GetNbinsX();
1546 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1548 hEff->SetStats(kFALSE);
1549 hEff->SetMinimum(0.);
1550 hEff->SetMaximum(110.);
1552 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1553 Double_t nGen = hGen->GetBinContent(iBin);
1554 Double_t nRec = hRec->GetBinContent(iBin);
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);
1561 if (error == 0) error = 0.0001;
1562 hEff->SetBinError(iBin, 100. * error);
1564 hEff->SetBinContent(iBin, 100. * 0.5);
1565 hEff->SetBinError(iBin, 100. * 0.5);
1573 TH1F* AliTPCComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
1574 Bool_t overflowBinFits)
1576 TVirtualPad* currentPad = gPad;
1577 TAxis* axis = hRes2->GetXaxis();
1578 Int_t nBins = axis->GetNbins();
1580 if (axis->GetXbins()->GetSize()){
1581 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1582 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1585 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1586 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1589 hRes->SetStats(false);
1590 hRes->SetOption("E");
1591 hRes->SetMinimum(0.);
1593 hMean->SetStats(false);
1594 hMean->SetOption("E");
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);
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;
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);
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);
1621 if (hBin->GetEntries() > 5) {
1622 fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1623 hBin->Fit(fitFunc,"s");
1624 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1627 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1628 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1631 hRes->SetBinContent(bin, 0.);
1632 hMean->SetBinContent(bin,0);
1634 hRes->SetBinError(bin, fitFunc->GetParError(2));
1635 hMean->SetBinError(bin, fitFunc->GetParError(1));
1641 hRes->SetBinContent(bin, 0.);
1642 hRes->SetBinError(bin, 0.);
1643 hMean->SetBinContent(bin, 0.);
1644 hMean->SetBinError(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());
1655 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1656 axis->GetTitle(), axis->GetBinUpEdge(bin));
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();