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 "AliDetector.h"
99 #include "AliTrackReference.h"
101 #include "AliTPCParamSR.h"
102 #include "AliTracker.h"
103 #include "AliComplexCluster.h"
104 #include "AliTPCComparisonMI.h"
117 AliTPCParam * GetTPCParam(){
118 AliTPCParamSR * par = new AliTPCParamSR;
124 //_____________________________________________________________________________
125 Float_t TPCBetheBloch(Float_t bg)
128 // Bethe-Bloch energy loss formula
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;
136 Double_t dbg = (Double_t) bg;
138 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
140 Double_t aa = TMath::Power(beta,kp4);
141 Double_t bb = TMath::Power(1./dbg,kp5);
143 bb=TMath::Log(kp3+bb);
145 return ((Float_t)((kp2-aa-bb)*kp1/aa));
152 ////////////////////////////////////////////////////////////////////////
153 AliTPCGenInfo::AliTPCGenInfo()
156 fTRdecay.SetTrack(-1);
159 AliTPCGenInfo::~AliTPCGenInfo()
161 if (fReferences) delete fReferences;
166 void AliTPCGenInfo::Update()
171 if (!fReferences) return;
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){
181 direction = newdirection;
188 ////////////////////////////////////////////////////////////////////////
190 ////////////////////////////////////////////////////////////////////////
192 // End of implementation of the class AliTPCGenInfo
194 ////////////////////////////////////////////////////////////////////////
198 ////////////////////////////////////////////////////////////////////////
203 ////////////////////////////////////////////////////////////////////////
204 digitRow & digitRow::operator=(const digitRow &digOld)
206 for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
209 ////////////////////////////////////////////////////////////////////////
210 void digitRow::SetRow(Int_t row)
212 if (row >= 8*kgRowBytes) {
213 cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
221 ////////////////////////////////////////////////////////////////////////
222 Bool_t digitRow::TestRow(Int_t row)
225 // return kTRUE if row is on
229 return TESTBIT(fDig[iC],iB);
231 ////////////////////////////////////////////////////////////////////////
232 Int_t digitRow::RowsOn(Int_t upto)
235 // returns number of rows with a digit
236 // count only rows less equal row number upto
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++;
247 ////////////////////////////////////////////////////////////////////////
248 void digitRow::Reset()
251 // resets all rows to zero
253 for (Int_t i = 0; i<kgRowBytes; i++) {
257 ////////////////////////////////////////////////////////////////////////
258 Int_t digitRow::Last()
261 // returns the last row number with a digit
262 // returns -1 if now digits
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;
271 ////////////////////////////////////////////////////////////////////////
272 Int_t digitRow::First()
275 // returns the first row number with a digit
276 // returns -1 if now digits
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;
286 ////////////////////////////////////////////////////////////////////////
288 // end of implementation of a class digitRow
290 ////////////////////////////////////////////////////////////////////////
292 ////////////////////////////////////////////////////////////////////////
293 TPCFindGenTracks::TPCFindGenTracks()
295 fMCInfo = new AliTPCGenInfo;
296 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
300 ////////////////////////////////////////////////////////////////////////
301 TPCFindGenTracks::TPCFindGenTracks(const char * fnGalice, const char* fnRes,
302 Int_t nEvents, Int_t firstEvent)
305 fFirstEventNr = firstEvent;
306 fEventNr = firstEvent;
309 sprintf(fFnRes,"%s",fnRes);
311 fLoader = AliRunLoader::Open(fnGalice);
313 delete gAlice->GetRunLoader();
317 if (fLoader->LoadgAlice()){
318 cerr<<"Error occured while l"<<endl;
320 Int_t nall = fLoader->GetNumberOfEvents();
328 cerr<<"no events available"<<endl;
332 if (firstEvent+nEvents>nall) {
333 fEventNr = nall-firstEvent;
334 cerr<<"restricted number of events availaible"<<endl;
338 ////////////////////////////////////////////////////////////////////////
339 void TPCFindGenTracks::Reset()
345 fContainerDigitRow = 0;
348 fReferenceIndex0 = 0;
349 fReferenceIndex1 = 0;
358 ////////////////////////////////////////////////////////////////////////
359 TPCFindGenTracks::~TPCFindGenTracks()
363 fLoader->UnloadgAlice();
370 Int_t TPCFindGenTracks::SetIO()
374 CreateTreeGenTracks();
375 if (!fTreeGenTracks) return 1;
376 AliTracker::SetFieldFactor();
377 fParamTPC = GetTPCParam();
382 ////////////////////////////////////////////////////////////////////////
383 Int_t TPCFindGenTracks::SetIO(Int_t eventNr)
388 fLoader->SetEventNumber(eventNr);
390 fLoader->LoadHeader();
391 fLoader->LoadKinematics();
392 fStack = fLoader->Stack();
394 fLoader->LoadTrackRefs();
395 fTreeTR = fLoader->TreeTR();
397 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
399 fTreeD = tpcl->TreeD();
403 Int_t TPCFindGenTracks::CloseIOEvent()
405 fLoader->UnloadHeader();
406 fLoader->UnloadKinematics();
407 fLoader->UnloadTrackRefs();
408 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
409 tpcl->UnloadDigits();
413 Int_t TPCFindGenTracks::CloseIO()
415 fLoader->UnloadgAlice();
421 ////////////////////////////////////////////////////////////////////////
422 Int_t TPCFindGenTracks::Exec(Int_t nEvents, Int_t firstEventNr)
425 fFirstEventNr = firstEventNr;
429 ////////////////////////////////////////////////////////////////////////
430 Int_t TPCFindGenTracks::Exec()
434 Int_t status =SetIO();
435 if (status>0) return status;
438 for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
441 fNParticles = fStack->GetNtrack();
442 fContainerDigitRow = new digitRow[fNParticles];
444 fReferences = new AliTrackReference[fgMaxTR];
445 fReferenceIndex0 = new Int_t[fNParticles];
446 fReferenceIndex1 = new Int_t[fNParticles];
447 fDecayRef = new AliTrackReference[fNParticles];
449 for (Int_t i = 0; i<fNParticles; i++) {
450 fReferenceIndex0[i] = -1;
451 fReferenceIndex1[i] = -1;
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;
464 delete [] fContainerDigitRow;
465 delete [] fReferences;
466 delete [] fReferenceIndex0;
467 delete [] fReferenceIndex1;
475 cerr<<"Exec finished"<<endl;
481 ////////////////////////////////////////////////////////////////////////
482 void TPCFindGenTracks::CreateTreeGenTracks()
484 fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
485 if (!fFileGenTracks) {
486 cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
489 fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
493 fMCInfo = new AliTPCGenInfo;
494 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
496 fTreeGenTracks->Branch("MC","AliTPCGenInfo",&fMCInfo);
498 fTreeGenTracks->Branch("fEventNr",&fEventNr,"fEventNr/I");
500 fTreeGenTracks->AutoSave();
502 ////////////////////////////////////////////////////////////////////////
503 void TPCFindGenTracks::CloseOutputFile()
505 if (!fFileGenTracks) {
506 cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
509 fFileGenTracks->cd();
510 fTreeGenTracks->Write();
511 delete fTreeGenTracks;
512 fFileGenTracks->Close();
513 delete fFileGenTracks;
517 ////////////////////////////////////////////////////////////////////////
518 Int_t TPCFindGenTracks::TreeKLoop()
521 // open the file with treeK
522 // loop over all entries there and save information about some tracks
525 AliStack * stack = fStack;
526 if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
529 cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
534 TParticlePDG *ppdg = 0;
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
541 //thank you very much root for this
542 TBranch * br = fTreeGenTracks->GetBranch("MC");
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++) {
550 // load only particles with TR
551 if (fReferenceIndex0[iParticle]<0) continue;
552 //////////////////////////////////////////////////////////////////////
553 fMCInfo->fLabel = iParticle;
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;
568 fMCInfo->fParticle = *(stack->Particle(iParticle));
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]);
579 Int_t index = fReferenceIndex0[iParticle];
580 AliTrackReference ref = fReferences[index];
581 // if (ref.GetTrack()!=iParticle)
582 // printf("problem2\n");
583 fMCInfo->fTrackRef = ref;
586 if (fMCInfo->fReferences !=0) delete fMCInfo->fReferences;
587 fMCInfo->fReferences =
588 new TClonesArray("AliTrackReference");
589 fMCInfo->fReferences->ExpandCreateFast(fReferenceIndex1[iParticle]-fReferenceIndex0[iParticle]+1);
591 for (Int_t i = fReferenceIndex0[iParticle];i<=fReferenceIndex1[iParticle];i++){
592 AliTrackReference ref = fReferences[i];
593 AliTrackReference *ref2 = (AliTrackReference*) fMCInfo->fReferences->At(rfindex);
594 if (ref.GetTrack()!=iParticle){
595 //printf("problem5\n");
603 ipdg = fMCInfo->fParticle.GetPdgCode();
604 fMCInfo->fPdg = ipdg;
605 ppdg = pdg.GetParticle(ipdg);
606 // calculate primary ionization per cm
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());
613 // Float_t betagama = fMCInfo->fParticle.P()/mass;
614 Float_t betagama = p /mass;
615 fMCInfo->fPrim = TPCBetheBloch(betagama);
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;
628 fMCInfo->fTRdecay.SetTrack(-1);
629 fMCInfo->fDecayCoord[0] = 0;
630 fMCInfo->fDecayCoord[1] = 0;
631 fMCInfo->fDecayCoord[2] = 0;
634 //////////////////////////////////////////////////////////////////////
636 br->SetAddress(&fMCInfo);
637 fTreeGenTracks->Fill();
640 fTreeGenTracks->AutoSave();
642 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
650 ////////////////////////////////////////////////////////////////////////
651 Int_t TPCFindGenTracks::TreeDLoop()
654 // open the file with treeD
655 // loop over all entries there and save information about some tracks
658 Int_t nInnerSector = fParamTPC->GetNInnerSector();
660 // Int_t zero=fParamTPC->GetZeroSup();
661 Int_t zero=fParamTPC->GetZeroSup()+6;
664 //char treeDName[100];
665 //sprintf(treeDName,"TreeD_75x40_100x60_150x60_%d",fEventNr);
666 //TTree *treeD=(TTree*)fFileTreeD->Get(treeDName);
669 AliSimDigits digitsAddress, *digits=&digitsAddress;
670 fTreeD->GetBranch("Segment")->SetAddress(&digits);
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;
678 fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
679 if (fDebug > 1) cout<<sec<<' '<<row<<" \r";
680 // cerr<<sec<<' '<<row<<endl;
682 // here I expect that upper sectors follow lower sectors
683 if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
685 //digits->ExpandTrackBuffer();
686 //Int_t *tracks = digits->GetTracks();
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) {
696 for (Int_t j = 0; j<3; j++) {
697 label = digits->GetTrackID(iRow,iColumn,j);
698 //label = digits->GetTrackIDFast(iRow,iColumn,j);
700 if (label >= fNParticles) { //don't label from bakground event
703 if (label >= 0 && label <= fNParticles) {
704 // if (label >= 0 && label <= fDebug) {
706 cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
708 <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
711 fContainerDigitRow[label].SetRow(row+rowShift);
715 } while (digits->Next());
718 if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
724 ////////////////////////////////////////////////////////////////////////
725 Int_t TPCFindGenTracks::TreeTRLoop()
728 // loop over TrackReferences and store the first one for each track
731 TTree * treeTR = fTreeTR;
733 cerr<<"TreeTR not found"<<endl;
736 Int_t nPrimaries = (Int_t) treeTR->GetEntries();
737 if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
740 //track references for TPC
741 TBranch *TPCBranchTR = treeTR->GetBranch("TPC");
743 cerr<<"TPC branch in TR not found"<<endl;
746 TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
747 TPCBranchTR->SetAddress(&TPCArrayTR);
748 //get decay point if exist
749 TBranch *runbranch = treeTR->GetBranch("AliRun");
751 cerr<<"Run branch in TR not found"<<endl;
754 TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
755 runbranch->SetAddress(&RunArrayTR);
760 for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
761 TPCBranchTR->GetEntry(iPrimPart);
763 for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
764 AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
766 Int_t label = trackRef->GetTrack();
767 if (label<0 || label > fNParticles) {
770 if (fReferenceIndex0[label]<0) ptstart = trackRef->Pt(); //store pt at the TPC entrance
771 if (ptstart<fgPtCut) continue;
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
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);
784 if (trackRef->Pt() < fgPtCut) continue;
785 Int_t label = trackRef->GetTrack();
786 if (label<0 || label > fNParticles) {
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
793 if (label>=fgMaxTR) continue; //restricted size of buffer for TR
794 fDecayRef[label] = *trackRef;
798 TPCArrayTR->Delete();
800 RunArrayTR->Delete();
806 ////////////////////////////////////////////////////////////////////////
807 Float_t TPCFindGenTracks::TR2LocalX(AliTrackReference *trackRef,
808 AliTPCParam *paramTPC) {
810 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
812 paramTPC->Transform0to1(x,index);
813 paramTPC->Transform1to2(x,index);
816 ////////////////////////////////////////////////////////////////////////
821 ////////////////////////////////////////////////////////////////////////
827 ////////////////////////////////////////////////////////////////////////
828 TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
830 const char* fnGalice,
831 Int_t nEvents, Int_t firstEvent)
834 // fFnGenTracks = fnGenTracks;
836 sprintf(fFnGenTracks,"%s",fnGenTracks);
837 sprintf(fFnCmp,"%s",fnCmp);
839 fFirstEventNr = firstEvent;
840 fEventNr = firstEvent;
844 fLoader = AliRunLoader::Open(fnGalice);
846 //delete gAlice->GetRunLoader();
850 if (fLoader->LoadgAlice()){
851 cerr<<"Error occured while l"<<endl;
853 Int_t nall = fLoader->GetNumberOfEvents();
861 cerr<<"no events available"<<endl;
865 if (firstEvent+nEvents>nall) {
866 fEventNr = nall-firstEvent;
867 cerr<<"restricted number of events availaible"<<endl;
872 ////////////////////////////////////////////////////////////////////////
873 TPCCmpTr::~TPCCmpTr()
880 //////////////////////////////////////////////////////////////
881 Int_t TPCCmpTr::SetIO()
886 if (!fTreeCmp) return 1;
887 AliTracker::SetFieldFactor();
888 fParamTPC = GetTPCParam();
890 if (!ConnectGenTree()) {
891 cerr<<"Cannot connect tree with generated tracks"<<endl;
897 //////////////////////////////////////////////////////////////
899 Int_t TPCCmpTr::SetIO(Int_t eventNr)
904 //gAlice->GetEvent(eventNr);
905 fLoader->SetEventNumber(eventNr);
907 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
909 fTreeRecTracks = tpcl->TreeT();
916 ////////////////////////////////////////////////////////////////////////
917 void TPCCmpTr::Reset()
922 // fFnCmp = "cmpTracks.root";
935 ////////////////////////////////////////////////////////////////////////
936 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
939 fFirstEventNr = firstEventNr;
943 ////////////////////////////////////////////////////////////////////////
944 Int_t TPCCmpTr::Exec()
952 fNextTreeGenEntryToRead = 0;
953 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
954 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
957 fNParticles = gAlice->GetEvent(fEventNr);
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;
968 fFakeRecTracks[i] = 0;
969 fMultiRecTracks[i] = 0;
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;
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;
981 delete fTreeRecTracks;
982 delete [] fIndexRecTracks;
983 delete [] fFakeRecTracks;
984 delete [] fMultiRecTracks;
989 cerr<<"Exec finished"<<endl;
995 ////////////////////////////////////////////////////////////////////////
996 Bool_t TPCCmpTr::ConnectGenTree()
999 // connect all branches from the genTracksTree
1000 // use the same variables as for the new cmp tree, it may work
1002 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
1003 if (!fFileGenTracks) {
1004 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
1007 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
1008 if (!fTreeGenTracks) {
1009 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
1010 <<fFnGenTracks<<endl;
1014 fMCInfo = new AliTPCGenInfo;
1015 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1016 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
1019 //fTreeGenTracks->SetBranchAddress("fEventNr",&fEventNr);
1023 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
1029 ////////////////////////////////////////////////////////////////////////
1030 void TPCCmpTr::CreateTreeCmp()
1032 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
1034 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
1039 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
1041 fMCInfo = new AliTPCGenInfo;
1042 fMCInfo->fReferences = new TClonesArray("AliTrackReference");
1043 fRecInfo = new AliTPCRecInfo;
1045 fTPCTrack = new AliTPCtrack;
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);
1052 fTreeCmp->AutoSave();
1055 ////////////////////////////////////////////////////////////////////////
1056 void TPCCmpTr::CloseOutputFile()
1059 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
1070 ////////////////////////////////////////////////////////////////////////
1072 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
1073 AliTPCParam *paramTPC) {
1075 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1077 paramTPC->Transform0to1(x,index);
1078 paramTPC->Transform1to2(x,index);
1081 ////////////////////////////////////////////////////////////////////////
1083 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
1086 // loop over all TPC reconstructed tracks and store info in memory
1091 if (!fTreeRecTracks) {
1092 cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
1095 //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
1097 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
1098 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
1100 TBranch * br= fTreeRecTracks->GetBranch("tracks");
1101 br->SetAddress(&fTPCTrack);
1103 if (fTreePoints) brp = fTreePoints->GetBranch("debug");
1110 fTrackPoints->Delete();
1111 delete fTrackPoints;
1114 fTracks = new TObjArray(nEntries);
1116 fTrackPoints = new TObjArray(nEntries);
1118 else fTrackPoints = 0;
1121 //load tracks to the memory
1122 for (Int_t i=0; i<nEntries;i++){
1123 AliTPCtrack * track = new AliTPCtrack;
1124 br->SetAddress(&track);
1126 fTracks->AddAt(track,i);
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);
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);
1146 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
1147 //br->GetEntry(iEntry);
1148 fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
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]++;
1156 if (fMultiRecTracks[absLabel]>0){
1157 if (fMultiRecTracks[absLabel]<4)
1158 fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
1161 fIndexRecTracks[absLabel*4] = iEntry;
1162 fMultiRecTracks[absLabel]++;
1165 printf("Time spended in TreeTLoop\n");
1168 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
1172 ////////////////////////////////////////////////////////////////////////
1173 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
1176 // loop over all entries for a given event, find corresponding
1177 // rec. track and store in the fTreeCmp
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");
1186 branch->SetAddress(&fRecInfo); // set all pointers
1188 while (entry < nParticlesTR) {
1189 fTreeGenTracks->GetEntry(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;
1201 if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
1202 fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
1204 // if (nBytes > 0) {
1207 TVector3 local = TR2Local(&(fMCInfo->fTrackRef),fParamTPC);
1208 local.GetXYZ(fRecInfo->fTRLocalCoord);
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)
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
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();
1235 fRecInfo->fTPCTrack.PropagateTo(local.X());
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]);
1252 fTreeCmp->AutoSave();
1256 fTrackPoints->Delete();
1257 delete fTrackPoints;
1260 printf("Time spended in TreeKLoop\n");
1262 if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1266 ////////////////////////////////////////////////////////////////////////
1267 ////////////////////////////////////////////////////////////////////////
1269 void AliTPCComparisonDraw::SetIO(const char *fname)
1272 TFile* file = TFile::Open(fname);
1274 printf("Could not open file - generated new one\n");
1275 TFile* filegen = TFile::Open("genTracks.root");
1277 printf("FILE with MC information is generated\n");
1278 TPCFindGenTracks *t = new TPCFindGenTracks("galice.root","genTracks.root",0);
1282 filegen = TFile::Open("genTracks.root");
1284 printf("COMPARISON FILE COULDNT BE GENERATED \n");
1287 printf("COMPARISON FILE IS GENERATED \n");
1288 TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",0);
1291 file = TFile::Open(fname);
1293 printf("Comparison file couldn't be generated\n");
1297 fTree = (TTree*) file->Get("TPCcmpTracks");
1299 printf("no track comparison tree found\n");
1305 void AliTPCComparisonDraw::ResPt()
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();
1317 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1318 p2->Draw();p2->cd();
1319 p2->SetGridx(); p2->SetGridy();
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");
1328 TH1F* hisp = ResPtvsPt("MC.fRowsWithDigits>100"+cteta1+cpos1,"1",0.15,2.,6);
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");
1347 void AliTPCComparisonDraw::Eff()
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();
1358 TPad *p2=new TPad("p2","p2",0.01,0.01,.99,.49);
1359 p2->Draw();p2->cd();
1360 p2->SetGridx(); p2->SetGridy();
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");
1369 TH1F* hisp = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cprim,"RC.fTPCTrack.fN>50");
1371 //hisp->DrawClone();
1372 TText * text = new TText(0.25,102.,"Primary particles");
1373 text->SetTextSize(0.05);
1377 TH1F* hiss = EffVsPt("MC.fRowsWithDigits>100"+cteta1+cpos1+cnprim,"RC.fTPCTrack.fN>50");
1379 //hiss->DrawClone();
1380 text = new TText(0.25,102.,"Secondary particles");
1381 text->SetTextSize(0.05);
1388 TH1F * AliTPCComparisonDraw::EffVsPt(const char* selection, const char * quality, Float_t min, Float_t max)
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);
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");
1402 TH1F* hEff = CreateEffHisto(hGen, hRec);
1403 AliLabelAxes(hEff, "p_{t} [GeV/c]", "#epsilon [%]");
1412 TH1F * AliTPCComparisonDraw::EffVsRM(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1415 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1416 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
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");
1423 TH1F* hEff = CreateEffHisto(hGen, hRec);
1424 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1431 TH1F * AliTPCComparisonDraw::EffVsRS(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1434 TH1F* hGen = new TH1F("hGen", "gen. tracks", nBins, min, max);
1435 TH1F* hRec = new TH1F("hRec", "rec. tracks", nBins, min, max);
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");
1442 TH1F* hEff = CreateEffHisto(hGen, hRec);
1443 AliLabelAxes(hEff, "r_{vertex} [cm]", "#epsilon [%]");
1451 TH1F * AliTPCComparisonDraw::ResPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
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);
1459 fTree->Draw("100.*(abs(1./fTPCTrack.Get1Pt())-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1462 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1463 AliLabelAxes(hRes, "p_{t} [GeV/c]", "#Delta p_{t} / p_{t} [%]");
1471 TH1F * AliTPCComparisonDraw::MeanPtvsPt(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
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);
1479 fTree->Draw("100.*(1./fTPCTrack.Get1Pt()-MC.fTrackRef.Pt())/MC.fTrackRef.Pt():MC.fTrackRef.Pt()>>hRes2", selection, "groff");
1482 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1483 AliLabelAxes(hRes, "p_{t} [GeV/c]", "mean p_{t} / p_{t} [%]");
1487 if (!hMean) return 0;
1493 TH1F * AliTPCComparisonDraw::ResdEdxvsN(const char* selection, const char * quality, Float_t min, Float_t max, Int_t nBins)
1496 Int_t nBinsRes = 15;
1498 TH2F* hRes2 = new TH2F("hRes2", "residuals", nBins, min, max, nBinsRes, 34, 60);
1500 fTree->Draw("RC.fTPCTrack.fdEdx/MC.fPrim:RC.fTPCTrack.fN>>hRes2", selection, "groff");
1503 TH1F* hRes = CreateResHisto(hRes2, &hMean);
1504 AliLabelAxes(hRes, "N points", "sigma dEdx/Nprim [%]");
1513 Double_t* AliTPCComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1515 Double_t* bins = new Double_t[nBins+1];
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];
1526 void AliTPCComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1529 histo->GetXaxis()->SetTitle(xAxisTitle);
1530 histo->GetYaxis()->SetTitle(yAxisTitle);
1534 TH1F* AliTPCComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1537 Int_t nBins = hGen->GetNbinsX();
1538 TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1540 hEff->SetStats(kFALSE);
1541 hEff->SetMinimum(0.);
1542 hEff->SetMaximum(110.);
1544 for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1545 Double_t nGen = hGen->GetBinContent(iBin);
1546 Double_t nRec = hRec->GetBinContent(iBin);
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);
1553 if (error == 0) error = 0.0001;
1554 hEff->SetBinError(iBin, 100. * error);
1556 hEff->SetBinContent(iBin, 100. * 0.5);
1557 hEff->SetBinError(iBin, 100. * 0.5);
1565 TH1F* AliTPCComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean, Bool_t drawBinFits,
1566 Bool_t overflowBinFits)
1568 TVirtualPad* currentPad = gPad;
1569 TAxis* axis = hRes2->GetXaxis();
1570 Int_t nBins = axis->GetNbins();
1572 if (axis->GetXbins()->GetSize()){
1573 hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1574 hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1577 hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1578 hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1581 hRes->SetStats(false);
1582 hRes->SetOption("E");
1583 hRes->SetMinimum(0.);
1585 hMean->SetStats(false);
1586 hMean->SetOption("E");
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);
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;
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);
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);
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));
1619 hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1620 hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1623 hRes->SetBinContent(bin, 0.);
1624 hMean->SetBinContent(bin,0);
1626 hRes->SetBinError(bin, fitFunc->GetParError(2));
1627 hMean->SetBinError(bin, fitFunc->GetParError(1));
1633 hRes->SetBinContent(bin, 0.);
1634 hRes->SetBinError(bin, 0.);
1635 hMean->SetBinContent(bin, 0.);
1636 hMean->SetBinError(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());
1647 sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1648 axis->GetTitle(), axis->GetBinUpEdge(bin));
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();