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 //
29 .L $ALICE_ROOT/STEER/AliGenInfo.C+
30 .L $ALICE_ROOT/STEER/AliTPCComparisonMI.C+
33 TPCCmpTr *t2 = new TPCCmpTr("genTracks.root","cmpTracks.root","galice.root",-1,1,0);
37 TCut cprim("cprim","MC.fVDist[3]<1");
38 TCut csec("cprim","MC.fVDist[3]>1");
39 TCut crec("crec","fReconstructed==1");
40 TCut cteta1("cteta1","abs(MC.fTrackRef.Theta()/3.1415-0.5)<0.25");
41 TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
42 TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
43 TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
44 AliTPCComparisonDraw comp;
49 comp.DrawXY("fTPCinP0[3]","fTPCDelta[4]/fTPCinP1[3]","fReconstructed==1&&fPdg==-211"+cprim,"1",4,0.2,1.5,-0.06,0.06)
52 comp.Eff("fTPCinP0[3]","fRowsWithDigits>120"+cteta1+cpos1,"fReconstructed>0",10,0,1.5)
60 #if !defined(__CINT__) || defined(__MAKECINT__)
70 #include "TBenchmark.h"
71 #include "TStopwatch.h"
72 #include "TParticle.h"
88 #include "AliTPCtrack.h"
89 #include "AliSimDigits.h"
90 #include "AliTPCParam.h"
92 #include "AliTPCLoader.h"
93 #include "AliDetector.h"
94 #include "AliTrackReference.h"
96 #include "AliTPCParamSR.h"
97 #include "AliTracker.h"
98 #include "AliComplexCluster.h"
101 #include "AliGenInfo.h"
102 #include "AliTPCComparisonMI.h"
109 void AliTPCRecInfo::Update(AliMCInfo* info,AliTPCParam * param, Bool_t reconstructed, Int_t direction)
113 //calculates derived variables
116 if (reconstructed==kFALSE){
117 fReconstructed = kFALSE;
120 fReconstructed = kTRUE;
121 // Find the nearest track reference
123 Double_t radius = TMath::Sqrt(fTrack.GetX()*fTrack.GetX()+fTrack.GetY()*fTrack.GetY());
124 TClonesArray * tpcreferences = info->fTPCReferences;
125 Int_t nentries = tpcreferences->GetEntriesFast();
126 AliTrackReference * ref = 0;
127 Double_t dr = 1000000;
128 for (Int_t i=0;i<nentries;i++){
129 AliTrackReference * refi = (AliTrackReference*)tpcreferences->UncheckedAt(i);
131 //if ( (direction<0) && (refi->R()<radius-10)) continue;
132 //if ( (direction>0) && (refi->R()>radius+10)) continue;
133 if (TMath::Abs(refi->R()-radius)<dr){
134 dr = TMath::Abs(refi->R()-radius);
139 fReconstructed = kFALSE;
144 fdEdx = fTrack.GetdEdx();
146 // propagate track to nearest reference in direction
147 TVector3 local = TPCCmpTr::TR2Local(ref,param);
148 local.GetXYZ(fTRLocalCoord);
150 // rotate if neccessary
151 Float_t ymax = local.X()*TMath::Tan(0.5*param->GetInnerAngle());
152 Float_t y = fTrack.GetYat(local.X());
154 fTrack.Rotate(param->GetInnerAngle());
155 } else if (y <-ymax) {
156 fTrack.Rotate(-param->GetInnerAngle());
158 Double_t xtrack = fTrack.GetX();
159 Double_t delta = (local.X()-xtrack)*0.1;
160 for (Int_t i=0;i<9;i++){
161 fTrack.PropagateTo(xtrack+delta*float(i));
163 fTrack.PropagateTo(local.X());
165 Double_t par[5], cov[15], localX;
166 fTrack.GetExternalParameters(localX,par);
167 fTrack.GetExternalCovariance(cov);
170 fRecPhi=TMath::ASin(par[2]) + fTrack.GetAlpha();
171 Double_t phi2 =TMath::ATan2(ref->Py(),ref->Px());
172 if (phi2<0) phi2+=2*TMath::Pi();
175 if (fRecPhi<0) fRecPhi+=2*TMath::Pi();
176 if (fRecPhi>=2*TMath::Pi()) fRecPhi-=2*TMath::Pi();
177 fLambda = TMath::ATan(par[3]);
178 fRecPt_1 = TMath::Abs(par[4]);
182 Double_t phi=TMath::ATan2(par[0],localX) + fTrack.GetAlpha();
183 Double_t r=TMath::Sqrt(localX*localX + par[0]*par[0]);
184 fTPCinR1[0]= r*TMath::Cos(phi);
185 fTPCinR1[1]= r*TMath::Sin(phi);
187 fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
188 fTPCinR1[4] = TMath::ATan2(fTPCinR1[1],fTPCinR1[0]);
190 fTPCinP1[0] = fTrack.Px();
191 fTPCinP1[1] = fTrack.Py();
192 fTPCinP1[2] = fTrack.Pz();
193 fTPCinP1[3] = TMath::Sqrt(fTPCinP1[0]*fTPCinP1[0]+fTPCinP1[1]*fTPCinP1[1]);
195 fTPCinR0[0] = ref->X();
196 fTPCinR0[1] = ref->Y();
197 fTPCinR0[2] = ref->Z();
198 fTPCinR0[3] = TMath::Sqrt(fTPCinR0[0]*fTPCinR0[0]+fTPCinR0[1]*fTPCinR0[1]);
199 fTPCinR0[4] = TMath::ATan2(fTPCinR0[1],fTPCinR0[0]);
201 fTPCinP0[0]=ref->Px();
202 fTPCinP0[1]=ref->Py();
203 fTPCinP0[2]=ref->Pz();
204 fTPCinP0[3]=ref->Pt();
207 if (fTPCinP1[3]>0.0000001){
208 fTPCAngle1[0] = TMath::ATan2(fTPCinP1[1],fTPCinP1[0]);
209 fTPCAngle1[1] = TMath::ATan(fTPCinP1[2]/fTPCinP1[3]);
211 fTPCAngle0[0] = TMath::ATan2(fTPCinP0[1],fTPCinP0[0]);
212 fTPCAngle0[1] = TMath::ATan(fTPCinP0[2]/fTPCinP0[3]);
215 fTPCDelta[0] = (fTPCinR0[4]-fTPCinR1[4])*fTPCinR1[3]; //delta rfi
216 fTPCPools[0] = fTPCDelta[0]/TMath::Sqrt(cov[0]);
217 fTPCDelta[1] = (fTPCinR0[2]-fTPCinR1[2]); //delta z
218 fTPCPools[1] = fTPCDelta[1]/TMath::Sqrt(cov[2]);
220 fTPCDelta[2] = (fTPCAngle0[0]-fTPCAngle1[0]);
221 fTPCPools[2] = fTPCDelta[2]/TMath::Sqrt(cov[5]);
223 fTPCDelta[3] = TMath::Tan(fTPCAngle0[1])-TMath::Tan(fTPCAngle1[1]);
224 fTPCPools[3] = fTPCDelta[3]/TMath::Sqrt(cov[9]);
225 fTPCDelta[4] = (fTPCinP0[3]-fTPCinP1[3]);
226 Double_t sign = (fTrack.GetC()>0) ? 1:-1;
227 fTPCPools[4] = sign*(1./fTPCinP0[3]-1./fTPCinP1[3])/TMath::Sqrt(cov[14]);
234 ////////////////////////////////////////////////////////////////////////
240 ////////////////////////////////////////////////////////////////////////
241 TPCCmpTr::TPCCmpTr(const char* fnGenTracks,
243 const char* fnGalice, Int_t direction,
244 Int_t nEvents, Int_t firstEvent)
247 // fFnGenTracks = fnGenTracks;
249 sprintf(fFnGenTracks,"%s",fnGenTracks);
250 sprintf(fFnCmp,"%s",fnCmp);
252 fFirstEventNr = firstEvent;
253 fEventNr = firstEvent;
255 fDirection = direction;
257 fLoader = AliRunLoader::Open(fnGalice);
259 //delete gAlice->GetRunLoader();
263 if (fLoader->LoadgAlice()){
264 cerr<<"Error occured while l"<<endl;
266 Int_t nall = fLoader->GetNumberOfEvents();
274 cerr<<"no events available"<<endl;
278 if (firstEvent+nEvents>nall) {
279 fEventNr = nall-firstEvent;
280 cerr<<"restricted number of events availaible"<<endl;
282 AliMagF * magf = gAlice->Field();
283 AliTracker::SetFieldMap(magf);
288 ////////////////////////////////////////////////////////////////////////
289 TPCCmpTr::~TPCCmpTr()
296 //////////////////////////////////////////////////////////////
297 Int_t TPCCmpTr::SetIO()
302 if (!fTreeCmp) return 1;
303 fParamTPC = GetTPCParam();
305 if (!ConnectGenTree()) {
306 cerr<<"Cannot connect tree with generated tracks"<<endl;
312 //////////////////////////////////////////////////////////////
314 Int_t TPCCmpTr::SetIO(Int_t eventNr)
319 //gAlice->GetEvent(eventNr);
320 // fLoader->SetEventNumber(eventNr);
321 //fLoader = AliRunLoader::Open("galice.root");
322 fLoader->GetEvent(eventNr);
324 AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
326 fTreeRecTracks = tpcl->TreeT();
334 ////////////////////////////////////////////////////////////////////////
335 void TPCCmpTr::Reset()
353 ////////////////////////////////////////////////////////////////////////
354 Int_t TPCCmpTr::Exec(Int_t nEvents, Int_t firstEventNr)
357 fFirstEventNr = firstEventNr;
361 ////////////////////////////////////////////////////////////////////////
362 Int_t TPCCmpTr::Exec()
370 cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
371 for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
375 fNParticles = gAlice->GetEvent(fEventNr);
377 fIndexRecTracks = new Int_t[fNParticles*4]; //write at maximum 4 tracks corresponding to particle
378 fFakeRecTracks = new Int_t[fNParticles];
379 fMultiRecTracks = new Int_t[fNParticles];
380 for (Int_t i = 0; i<fNParticles; i++) {
381 fIndexRecTracks[4*i] = -1;
382 fIndexRecTracks[4*i+1] = -1;
383 fIndexRecTracks[4*i+2] = -1;
384 fIndexRecTracks[4*i+3] = -1;
386 fFakeRecTracks[i] = 0;
387 fMultiRecTracks[i] = 0;
390 cout<<"Start to process event "<<fEventNr<<endl;
391 cout<<"\tfNParticles = "<<fNParticles<<endl;
392 if (fDebug>2) cout<<"\tStart loop over TreeT"<<endl;
393 if (TreeTLoop(eventNr)>0) return 1;
395 if (fDebug>2) cout<<"\tStart loop over tree genTracks"<<endl;
396 if (TreeGenLoop(eventNr)>0) return 1;
397 if (fDebug>2) cout<<"\tEnd loop over tree genTracks"<<endl;
399 delete [] fIndexRecTracks;
400 delete [] fFakeRecTracks;
401 delete [] fMultiRecTracks;
406 cerr<<"Exec finished"<<endl;
412 ////////////////////////////////////////////////////////////////////////
413 Bool_t TPCCmpTr::ConnectGenTree()
416 // connect all branches from the genTracksTree
417 // use the same variables as for the new cmp tree, it may work
419 fFileGenTracks = TFile::Open(fFnGenTracks,"READ");
420 if (!fFileGenTracks) {
421 cerr<<"Error in ConnectGenTree: cannot open file "<<fFnGenTracks<<endl;
424 fTreeGenTracks = (TTree*)fFileGenTracks->Get("genTracksTree");
425 if (!fTreeGenTracks) {
426 cerr<<"Error in ConnectGenTree: cannot find genTracksTree in the file "
427 <<fFnGenTracks<<endl;
431 fMCInfo = new AliMCInfo;
432 fTreeGenTracks->SetBranchAddress("MC",&fMCInfo);
436 cout<<"Number of gen. tracks with TR: "<<fTreeGenTracks->GetEntries()<<endl;
442 ////////////////////////////////////////////////////////////////////////
443 void TPCCmpTr::CreateTreeCmp()
445 fFileCmp = TFile::Open(fFnCmp,"RECREATE");
447 cerr<<"Error in CreateTreeCmp: cannot open file "<<fFnCmp<<endl;
452 fTreeCmp = new TTree("TPCcmpTracks","TPCcmpTracks");
454 fMCInfo = new AliMCInfo;
455 fRecInfo = new AliTPCRecInfo;
457 fTPCTrack = new AliTPCtrack;
459 fTreeCmp->Branch("MC","AliMCInfo",&fMCInfo);
460 fTreeCmp->Branch("RC","AliTPCRecInfo",&fRecInfo);
461 fTreeCmp->Branch("fTPCTrack","AliTPCtrack",&fTPCTrack);
463 fTreeCmp->AutoSave();
466 ////////////////////////////////////////////////////////////////////////
467 void TPCCmpTr::CloseOutputFile()
470 cerr<<"File "<<fFnCmp<<" not found as an open file."<<endl;
481 ////////////////////////////////////////////////////////////////////////
483 TVector3 TPCCmpTr::TR2Local(AliTrackReference *trackRef,
484 AliTPCParam *paramTPC) {
486 Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
488 paramTPC->Transform0to1(x,index);
489 paramTPC->Transform1to2(x,index);
492 ////////////////////////////////////////////////////////////////////////
494 Int_t TPCCmpTr::TreeTLoop(Int_t eventNr)
497 // loop over all TPC reconstructed tracks and store info in memory
502 if (!fTreeRecTracks) {
503 cerr<<"Can't get a tree with TPC rec. tracks "<<endl;
506 //fTreePoints=(TTree*)fFileRecTracks->Get("trackDebug");
508 Int_t nEntries = (Int_t) fTreeRecTracks->GetEntries();
509 if (fDebug > 2) cout<<"Event, rec. tracks: "<<eventNr<<" "
511 TBranch * br= fTreeRecTracks->GetBranch("tracks");
512 br->SetAddress(&fTPCTrack);
514 if (fTreePoints) brp = fTreePoints->GetBranch("debug");
521 fTrackPoints->Delete();
525 fTracks = new TObjArray(nEntries);
527 fTrackPoints = new TObjArray(nEntries);
529 else fTrackPoints = 0;
532 //load tracks to the memory
533 for (Int_t i=0; i<nEntries;i++){
534 AliTPCtrack * track = new AliTPCtrack;
535 br->SetAddress(&track);
537 fTracks->AddAt(track,i);
540 //load track points to the memory
541 if (brp) for (Int_t i=0; i<nEntries;i++){
542 TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
543 brp->SetAddress(&arr);
546 for (Int_t j=0;j<arr->GetEntriesFast();j++){
547 AliTPCTrackPoint2 * point = (AliTPCTrackPoint2*)arr->UncheckedAt(j);
548 if (point && point->fID>=0){
549 fTrackPoints->AddAt(arr,point->fID);
557 for (Int_t iEntry=0; iEntry<nEntries;iEntry++){
558 //br->GetEntry(iEntry);
559 fTPCTrack = (AliTPCtrack*)fTracks->UncheckedAt(iEntry);
561 Int_t label = fTPCTrack->GetLabel();
562 Int_t absLabel = abs(label);
563 if (absLabel < fNParticles) {
564 // fIndexRecTracks[absLabel] = iEntry;
565 if (label < 0) fFakeRecTracks[absLabel]++;
567 if (fMultiRecTracks[absLabel]>0){
568 if (fMultiRecTracks[absLabel]<4)
569 fIndexRecTracks[absLabel*4+fMultiRecTracks[absLabel]] = iEntry;
572 fIndexRecTracks[absLabel*4] = iEntry;
573 fMultiRecTracks[absLabel]++;
576 printf("Time spended in TreeTLoop\n");
579 if (fDebug > 2) cerr<<"end of TreeTLoop"<<endl;
584 ////////////////////////////////////////////////////////////////////////
585 Int_t TPCCmpTr::TreeGenLoop(Int_t eventNr)
588 // loop over all entries for a given event, find corresponding
589 // rec. track and store in the fTreeCmp
594 Double_t nParticlesTR = fTreeGenTracks->GetEntriesFast();
595 TBranch * branch = fTreeCmp->GetBranch("RC");
596 branch->SetAddress(&fRecInfo); // set all pointers
598 while (entry < nParticlesTR) {
599 fTreeGenTracks->GetEntry(entry);
601 if (eventNr < fMCInfo->fEventNr) continue;
602 if (eventNr > fMCInfo->fEventNr) continue;
604 if (fDebug > 2 && fMCInfo->fLabel < 10) {
605 cerr<<"Fill track with a label "<<fMCInfo->fLabel<<endl;
609 if (fIndexRecTracks[fMCInfo->fLabel*4] >= 0) {
610 fTPCTrack= (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4]);
616 //if multifound find track with biggest pt - returning tracks - loos energy
618 if (fIndexRecTracks[fMCInfo->fLabel*4]+1){
619 Double_t pt = fTPCTrack-Pt();
620 for (Int_t i=1;i<4;i++){
621 if (fIndexRecTracks[fMCInfo->fLabel*4+i]>=0){
622 AliTPCtrack * track = (AliTPCtrack*)fTracks->UncheckedAt(fIndexRecTracks[fMCInfo->fLabel*4+i]);
623 Double_t pt2 = TMath::Abs(1/track->Get1Pt());
631 fRecInfo->fTrack =*fTPCTrack;
632 fRecInfo->fReconstructed = 1;
633 fRecInfo->fFake = fFakeRecTracks[fMCInfo->fLabel];
634 fRecInfo->fMultiple = fMultiRecTracks[fMCInfo->fLabel];
635 fRecInfo->Update(fMCInfo, fParamTPC, kTRUE, 1);
639 fRecInfo->Update(fMCInfo, fParamTPC, kFALSE, 1);
643 fTreeCmp->AutoSave();
647 fTrackPoints->Delete();
651 printf("Time spended in TreeGenLoop\n");
653 if (fDebug > 2) cerr<<"end of TreeGenLoop"<<endl;
658 ////////////////////////////////////////////////////////////////////////
659 ////////////////////////////////////////////////////////////////////////
661 void AliTPCComparisonDraw::SetIO(const char *fname)
663 TFile *file = TFile::Open(fname);
665 printf("Comparison file couldn't be generated\n");
669 fTree = (TTree*) file->Get("TPCcmpTracks");
671 printf("no track comparison tree found\n");
678 void AliTPCComparisonDraw::Eff()
683 void AliTPCComparisonDraw::ResPt()