]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MFT/AliMuonForwardTrackFinder.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes for coverity: SELF_ASSIGN
[u/mrichter/AliRoot.git] / MFT / AliMuonForwardTrackFinder.cxx
CommitLineData
820b4d9e 1// ROOT includes
2#include "TObject.h"
3#include "TClonesArray.h"
4#include "TObjArray.h"
5#include "TH2D.h"
6#include "TH1D.h"
7#include "TFile.h"
8#include "TGeoManager.h"
9#include "TMatrixD.h"
10#include "TParticle.h"
11#include "TMath.h"
12#include "TGraph.h"
13#include "TEllipse.h"
14#include "TCanvas.h"
15#include "TString.h"
16#include "TLatex.h"
17#include "TMarker.h"
18#include "TNtuple.h"
19#include "TRandom.h"
20#include "TIterator.h"
21
22// STEER includes
23#include "AliLog.h"
24#include "AliRun.h"
25#include "AliRunLoader.h"
26#include "AliLoader.h"
27#include "AliHeader.h"
28#include "AliMC.h"
29#include "AliStack.h"
30#include "AliMagF.h"
31#include "AliTracker.h"
32#include "AliGRPObject.h"
33#include "AliCDBEntry.h"
34#include "AliCDBManager.h"
35
36// MUON includes
37#include "AliMUONConstants.h"
38#include "AliMUONTrack.h"
39#include "AliMUONRecoCheck.h"
40#include "AliMUONTrackParam.h"
41#include "AliMUONTrackExtrap.h"
42#include "AliMUONVTrackStore.h"
43#include "AliMUONVCluster.h"
44
45// MFT includes
46#include "AliMuonForwardTrack.h"
47#include "AliMFTCluster.h"
48#include "AliMFT.h"
49#include "AliMFTSegmentation.h"
50
51#include "AliMuonForwardTrackFinder.h"
52
53//====================================================================================================================================================
54//
55// Class for the creation of the "global muon tracks" built from the clusters in the
56// muon spectrometer and the clusters of the Muon Forward Tracker. QA histograms are also created
57//
58// Contact author: antonio.uras@cern.ch
59//
60//====================================================================================================================================================
61
62ClassImp(AliMuonForwardTrackFinder)
63
64//=====================================================================================================
65
66AliMuonForwardTrackFinder::AliMuonForwardTrackFinder():
67 TObject(),
68 fRun(0),
69 fNEventsToAnalyze(0),
70 fSigmaClusterCut(0),
71 fChi2GlobalCut(0),
72 fSigmaSpectrometerCut(0),
73 fExtrapOriginTransvError(0),
74 fGaussianBlurZVert(0),
75 fNFinalCandidatesCut(0),
76 fReadDir(0),
77 fOutDir(0),
78 fDrawOption(0),
79
80 fDistanceFromGoodClusterAndTrackAtLastPlane(-1),
81 fDistanceFromBestClusterAndTrackAtLastPlane(-1),
82
83 fRAbsorberCut(0),
84 fLowPtCut(0),
85 fNPlanesMFT(0),
86 fNPlanesMFTAnalyzed(0),
87 fNMaxMissingMFTClusters(0),
88
89 fEv(0),
90 fLabelMC(0),
91
92 fHistPtSpectrometer(0),
93 fHistPtMuonTrackWithGoodMatch(0),
94 fHistPtMuonTrackWithBadMatch(0),
95 fHistRadiusEndOfAbsorber(0),
96 fHistNGoodClustersForFinalTracks(0),
97 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane(0),
98 fHistDistanceGoodClusterFromTrackAtLastPlane(0),
99
100 fNtuFinalCandidates1(0),
101 fNtuFinalBestCandidates1(0),
102 fNtuFinalCandidates2(0),
103 fNtuFinalBestCandidates2(0),
104
105 fCanvas(0),
106
107 fTxtMuonHistory(0),
108 fTxtTrackGoodClusters(0),
109 fTxtTrackFinalChi2(0),
110 fTxtFinalCandidates(0),
111 fTxtDummy(0),
112 fTxtAllClust(0),
113 fTxtClustGoodChi2(0),
114 fTxtClustMC(0),
115 fTxtClustOfTrack(0),
116 fMrkAllClust(0),
117 fMrkClustGoodChi2(0),
118 fMrkClustMC(0),
119 fMrkClustOfTrack(0),
120
121 fCountRealTracksAnalyzed(0),
122 fCountRealTracksWithRefMC(0),
123 fCountRealTracksWithRefMC_andTrigger(0),
124 fCountRealTracksWithRefMC_andTrigger_andGoodPt(0),
125 fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta(0),
126 fCountRealTracksAnalyzedOfEvent(0),
127 fCountRealTracksAnalyzedWithFinalCandidates(0),
128
129 fFileCluster(0),
130 fFileESD(0),
131 fFile_gAlice(0),
132
133 fRunLoader(0),
134 fMFTLoader(0),
135 fMuonRecoCheck(0),
136 fMFTClusterTree(0),
137 fMuonTrackReco(0),
138 fCurrentTrack(0),
139 fIsCurrentMuonTrackable(0),
140 fCandidateTracks(0),
141 fTrackStore(0),
142 fTrackRefStore(0),
143 fNextTrack(0),
144 fStack(0),
145 fMFT(0),
146 fSegmentation(0),
147 fOutputTreeFile(0),
148 fOutputEventTree(0),
149 fMuonForwardTracks(0),
150 fMatchingMode(-1),
151 fMinResearchRadiusAtLastPlane(0),
152 fGRPData(0),
153 fRunInfo(0)
154
155{
156
157 // Default constructor
158
159 for (Int_t iPlane=0; iPlane<fMaxNPlanesMFT; iPlane++) {
160
161 fHistNTracksAfterExtrapolation[iPlane] = 0;
162 fHistChi2Cluster_GoodCluster[iPlane] = 0;
163 fHistChi2Cluster_BadCluster[iPlane] = 0;
164 fHistResearchRadius[iPlane] = 0;
165
166 fIsGoodClusterInPlane[iPlane] = kFALSE;
167
168 fHistChi2Cluster_GoodCluster[iPlane] = 0;
169 fHistChi2Cluster_BadCluster[iPlane] = 0;
170
171 fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = 0;
172 fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = 0;
173
174 fZPlane[iPlane] = 0.;
175 fRPlaneMax[iPlane] = 0.;
176 fRPlaneMin[iPlane] = 0.;
177
178 for (Int_t i=0; i<4; i++) fGrMFTPlane[iPlane][i] = 0;
179 fCircleExt[iPlane] = 0;
180 fCircleInt[iPlane] = 0;
181
182 fTxtTrackChi2[iPlane] = 0;
183
184 fIsClusterCompatible[iPlane] = 0;
185
186 fMFTClusterArray[iPlane] = 0;
187 fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
188 fMFTClusterArrayBack[iPlane] = new TClonesArray("AliMFTCluster");
189
190 fIsPlaneMandatory[iPlane] = kFALSE;
191
192 }
193
194 // fNextTrack = 0;
195
196 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane = 0;
197 fHistDistanceGoodClusterFromTrackAtLastPlane = 0;
198
199 fMFTClusterTree = 0;
200 fCandidateTracks = 0;
201
202 fOutputTreeFile = new TFile("MuonGlobalTracks.root", "recreate");
203 fOutputEventTree = new TTree("AliMuonForwardTracks", "Tree of AliMuonForwardTracks");
204 fMuonForwardTracks = new TClonesArray("AliMuonForwardTrack");
205 fOutputEventTree -> Branch("tracks", &fMuonForwardTracks);
206
207}
208
209//=====================================================================================================
210
211void AliMuonForwardTrackFinder::Init(Int_t nRun,
212 Char_t *readDir,
213 Char_t *outDir,
214 Int_t nEventsToAnalyze) {
215
216 if (fRunLoader) {
217 printf("WARNING: run already initialized!!\n");
218 }
219
220 SetRun(nRun);
221 SetReadDir(readDir);
222 SetOutDir(outDir);
223
224 printf("input dir = %s\n", fReadDir.Data());
225 printf("output dir = %s\n", fOutDir.Data());
226
227 // -------------------------- initializing files...
228
229 printf("initializing files for run %d...\n", fRun);
230
231 Char_t geoFileName[300];
232 Char_t esdFileName[300];
233 Char_t gAliceName[300];
234 Char_t clusterName[300];
235
236 sprintf(geoFileName , "%s/geometry.root", fReadDir.Data());
237 sprintf(esdFileName , "%s/AliESDs.root" , fReadDir.Data());
238 sprintf(gAliceName , "%s/galice.root" , fReadDir.Data());
239 sprintf(clusterName , "%s/MFT.RecPoints.root", fReadDir.Data());
240
241 // Import TGeo geometry (needed by AliMUONTrackExtrap::ExtrapToVertex)
242 if (!gGeoManager) {
243 TGeoManager::Import(geoFileName);
244 if (!gGeoManager) {
245 printf("getting geometry from file %s failed", geoFileName);
246 return;
247 }
248 }
249
250 fFileESD = new TFile(esdFileName);
251 if (!fFileESD || !fFileESD->IsOpen()) return;
252 else printf("file %s successfully opened\n", fFileESD->GetName());
253
254 fMuonRecoCheck = new AliMUONRecoCheck(esdFileName, Form("%s/generated/", fReadDir.Data())); // Utility class to check reconstruction
255 fFile_gAlice = new TFile(gAliceName);
256 if (!fFile_gAlice || !fFile_gAlice->IsOpen()) return;
257 else printf("file %s successfully opened\n", fFile_gAlice->GetName());
258
259 fRunLoader = AliRunLoader::Open(gAliceName);
260 gAlice = fRunLoader->GetAliRun();
261 if (!gAlice) fRunLoader->LoadgAlice();
262 fMFT = (AliMFT*) gAlice->GetDetector("MFT");
263 fSegmentation = fMFT->GetSegmentation();
264 SetNPlanesMFT(fSegmentation->GetNPlanes());
265
266 if (!SetRunNumber()) return;
267 if (!InitGRP()) return;
268 AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations
269
270 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
271 fZPlane[iPlane] = fSegmentation->GetPlane(iPlane)->GetZCenter();
272 fRPlaneMax[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMaxSupport();
273 fRPlaneMin[iPlane] = fSegmentation->GetPlane(iPlane)->GetRMinSupport();
274 }
275
276 // Loading MFT clusters
277 fMFTLoader = fRunLoader->GetDetectorLoader("MFT");
278 fMFTLoader->LoadRecPoints("READ");
279
280 fMFTClusterTree = fMFTLoader->TreeR();
281
282
283 Int_t nEventsInFile = fMuonRecoCheck->NumberOfEvents();
284 if (!nEventsInFile) {
285 printf("no events available!!!\n");
286 return;
287 }
288 if (nEventsInFile<nEventsToAnalyze || nEventsToAnalyze<0) fNEventsToAnalyze = nEventsInFile;
289 else fNEventsToAnalyze = nEventsToAnalyze;
290
291 fCandidateTracks = new TClonesArray("AliMuonForwardTrack",200000);
292
293 // -------------------------- initializing histograms...
294
295 printf("\ninitializing histograms...\n");
296 BookHistos();
297 SetTitleHistos();
298 printf("... done!\n\n");
299
300 // -------------------------- initializing graphics...
301
302 printf("initializing graphics...\n");
303 BookPlanes();
304 printf("... done!\n\n");
305
306 SetSigmaSpectrometerCut(4.0);
307 SetSigmaClusterCut(4.5);
308 SetChi2GlobalCut(2.0);
309 SetNFinalCandidatesCut(10);
310 SetRAbsorberCut(26.4);
311 SetLowPtCut(0.5);
312 SetExtrapOriginTransvError(1.0);
313
314}
315
316//======================================================================================================================================
317
318Bool_t AliMuonForwardTrackFinder::LoadNextEvent() {
319
320 // load next reconstructed event from the tree
321
322 if (fEv) FillOutputTree();
323
324 if (fEv>=fNEventsToAnalyze) return kFALSE;
325
326 fCountRealTracksAnalyzedOfEvent = 0;
327
328 printf(" **** analyzing event # %d \n", fEv);
329
330 fTrackStore = fMuonRecoCheck->ReconstructedTracks(fEv);
331 fTrackRefStore = fMuonRecoCheck->ReconstructibleTracks(fEv);
332
333 fRunLoader->GetEvent(fEv);
334 if (!fMFTLoader->TreeR()->GetEvent()) return kFALSE;
335 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
336 printf("plane %02d: nClusters = %d\n", iPlane, (fMFT->GetRecPointsList(iPlane))->GetEntries());
337 fMFTClusterArray[iPlane] = fMFT->GetRecPointsList(iPlane);
338 }
339 SeparateFrontBackClusters();
340
341 fRunLoader -> LoadKinematics();
342 fStack = fRunLoader->Stack();
343 fNextTrack = fTrackStore->CreateIterator();
344 fMuonForwardTracks->Clear();
345
346 fEv++;
347
348 return kTRUE;
349
350}
351
352//======================================================================================================================================
353
354Int_t AliMuonForwardTrackFinder::LoadNextTrack() {
355
356 fNPlanesMFTAnalyzed = 0;
357
358 // load next muon track from the reconstructed event
359
360 if (!fCountRealTracksAnalyzed) if (!LoadNextEvent()) return kFALSE;
361
362 while ( !(fMuonTrackReco = static_cast<AliMUONTrack*>(fNextTrack->Next())) ) if (!LoadNextEvent()) return kFALSE;
363
364 printf("**************************************************************************************\n");
365 printf("*************************** MUON TRACK %3d ***************************************\n", fCountRealTracksAnalyzedOfEvent);
366 printf("**************************************************************************************\n");
367
368 fCountRealTracksAnalyzed++;
369
370 fCandidateTracks -> Clear();
371
372 fLabelMC = -1;
373 fDistanceFromGoodClusterAndTrackAtLastPlane = -1.;
374 fDistanceFromBestClusterAndTrackAtLastPlane = -1.;
375 ResetPlanes();
376
377 TIter nextTrackRef(fTrackRefStore->CreateIterator());
378 AliMUONTrack *trackRef=0;
379
380 // --------------------------------------- loop on MC generated tracks to find the MC reference...
381
382 while ( (trackRef = static_cast<AliMUONTrack*>(nextTrackRef())) ) {
383 // number of compatible clusters between trackReco and trackRef
384 Int_t nMatchCluster = fMuonTrackReco->FindCompatibleClusters(*trackRef, fSigmaSpectrometerCut, fIsClusterCompatible);
385 if ( (fIsClusterCompatible[0] || fIsClusterCompatible[1] || fIsClusterCompatible[2] || fIsClusterCompatible[3]) && // before the dipole
386 (fIsClusterCompatible[6] || fIsClusterCompatible[7] || fIsClusterCompatible[8] || fIsClusterCompatible[9]) && // after the dipole
387 2*nMatchCluster>fMuonTrackReco->GetNClusters() ) {
388 fMuonTrackReco->SetMCLabel(trackRef->GetUniqueID()); // MC reference has been found for trackReco!
389 break;
390 }
391 }
392
393 // ------------------------------------- ...done!
394
395 if (fMuonTrackReco->GetMCLabel()>=0) fCountRealTracksWithRefMC++;
396
397 fLabelMC = fMuonTrackReco->GetMCLabel();
398
399 CheckCurrentMuonTrackable();
400
401 PrintParticleHistory();
402
403 if (fMuonTrackReco->GetMatchTrigger()) fCountRealTracksWithRefMC_andTrigger++;
404
405 // the track we are going to build, starting from fMuonTrackReco and adding the MFT clusters
406 AliMuonForwardTrack *track = new ((*fCandidateTracks)[0]) AliMuonForwardTrack();
407 track -> SetMUONTrack(fMuonTrackReco);
408 if (fLabelMC>=0) track -> SetMCTrackRef(fStack->Particle(fLabelMC));
409 track -> SetMCLabel(fMuonTrackReco->GetMCLabel());
410 track -> SetMatchTrigger(fMuonTrackReco->GetMatchTrigger());
411
412 // track parameters at the first tracking station in the Muon Spectrometer
413 AliMUONTrackParam *param = (AliMUONTrackParam*) (fMuonTrackReco->GetTrackParamAtCluster()->First());
414 Double_t ptSpectrometer = TMath::Sqrt(param->Px()*param->Px() + param->Py()*param->Py());
415 Double_t thetaSpectrometer = TMath::ATan(ptSpectrometer/param->Pz());
416 if (thetaSpectrometer<0.) thetaSpectrometer += TMath::Pi();
417 Double_t etaSpectrometer = -1.*TMath::Log(TMath::Tan(0.5*thetaSpectrometer));
418 fHistPtSpectrometer -> Fill(ptSpectrometer);
419
420 // if the transverse momentum in the Muon Spectrometer is smaller than the threshold, skip to the next track
421 if (ptSpectrometer < fLowPtCut) return 3;
422
423 // track parameters linearly extrapolated from the first tracking station to the end of the absorber
424 AliMUONTrackParam trackParamEndOfAbsorber(*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
425 AliMUONTrackExtrap::ExtrapToZCov(&trackParamEndOfAbsorber, -503.); // absorber extends from -90 to -503 cm
426 Double_t xEndOfAbsorber = trackParamEndOfAbsorber.GetNonBendingCoor();
427 Double_t yEndOfAbsorber = trackParamEndOfAbsorber.GetBendingCoor();
428 Double_t rAbsorber = TMath::Sqrt(xEndOfAbsorber*xEndOfAbsorber + yEndOfAbsorber*yEndOfAbsorber);
429 fHistRadiusEndOfAbsorber -> Fill(rAbsorber);
430
431 // if the radial distance of the track at the end of the absorber is smaller than a radius corresponding to
432 // 3 degrees as seen from the interaction point, skip to the next track
433 if (rAbsorber < fRAbsorberCut) return 4;
434
435 //------------------------- NOW THE CYCLE OVER THE MFT PLANES STARTS ---------------------------------------
436
437 for (Int_t iPlane=fNPlanesMFT-1; iPlane>=0; iPlane--) { // *** do not reverse the order of this cycle!!!
438 // *** this reflects the fact that the extrapolation is performed
439 // *** starting from the last MFT plane back to the origin
440
441 // --------- updating the array of tracks according to the clusters available in the i-th plane ---------
442
443 fNPlanesMFTAnalyzed++;
444
445 if (fMatchingMode==kRealMatching) {
446 Int_t nTracksToBeAnalyzed = fCandidateTracks->GetEntriesFast();
447 for (Int_t iTrack=0; iTrack<nTracksToBeAnalyzed; iTrack++) {
448 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
449 // if the old track is compatible with the new cluster, the track is updated and inserted as new track in the array
450 // (several new tracks can be created for one old track)
451 FindClusterInPlane(iPlane);
452 if ((fNPlanesMFTAnalyzed-fCurrentTrack->GetNMFTClusters())>fNMaxMissingMFTClusters || fIsPlaneMandatory[iPlane]) {
453 fCandidateTracks->Remove(fCurrentTrack); // the old track is removed after the check;
454 }
455 }
456 fCandidateTracks->Compress();
457 if (fIsCurrentMuonTrackable) fHistNTracksAfterExtrapolation[iPlane] -> Fill(fCandidateTracks->GetEntriesFast());
458 }
459
460 else if (fMatchingMode==kIdealMatching) {
461 fCurrentTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
462 printf("plane %02d: fCandidateTracks->GetEntriesFast() = %d fCandidateTracks->UncheckedAt(0) = %p fCurrentTrack = %p\n",
463 iPlane, fCandidateTracks->GetEntriesFast(), fCandidateTracks->UncheckedAt(0), fCurrentTrack);
464 AttachGoodClusterInPlane(iPlane);
465 }
466
467 }
468
469 // -------------------------- END OF THE CYCLE OVER THE MFT PLANES --------------------------------------------
470
471 if (fMatchingMode==kIdealMatching) {
472 printf("Adding track to output tree...\n");
473 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(0);
474 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack); // AU
475 printf("...track added!\n");
476 fCandidateTracks->Clear();
477 fCountRealTracksAnalyzedOfEvent++;
478 return 5;
479 }
480
481 // If we have several final tracks, we must find the best candidate:
482
483 Int_t nFinalTracks = fCandidateTracks->GetEntriesFast();
484
485 if (nFinalTracks) fCountRealTracksAnalyzedWithFinalCandidates++;
486
487 Double_t theVariable_Best = -1.; // variable defining the best candidate
488 Bool_t bestCandidateExists = kFALSE;
489 Int_t nGoodClustersBestCandidate = 0;
490 Int_t idBestCandidate = 0;
491 Double_t chi2HistoryForBestCandidate[fMaxNPlanesMFT]={0}; // chi2 on each plane, for the best candidate
492 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = -1.;
493
494 fTxtFinalCandidates = new TLatex(0.10, 0.78, Form("N_{FinalCandidates} = %d", nFinalTracks));
495
496 Int_t nClustersMC = 0;
497 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) nClustersMC += fIsGoodClusterInPlane[iPlane];
498
499 for (Int_t iTrack=0; iTrack<nFinalTracks; iTrack++) {
500
501 AliMuonForwardTrack *finalTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(iTrack);
502
503 Double_t chi2AtPlane[fMaxNPlanesMFT]={0};
504 Int_t nGoodClusters = 0;
505 Int_t nMFTClusters = finalTrack->GetNMFTClusters();
506 Int_t nMUONClusters = finalTrack->GetNMUONClusters();
507
508 Int_t plane = 0;
509 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) {
510 while (!finalTrack->PlaneExists(plane)) plane++;
511 AliMFTCluster *localCluster = finalTrack->GetMFTCluster(iCluster);
512 chi2AtPlane[plane++] = localCluster->GetTrackChi2();
513 if (IsCorrectMatch(localCluster)) nGoodClusters++;
514 Int_t nClustersGlobalTrack = nMUONClusters + (nMFTClusters-iCluster); // Muon Spectrometer clusters + clusters in the Vertex Telescope
515 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
516 chi2AtPlane[plane] /= Double_t(ndfGlobalTrack);
517 }
518
519 if (fIsCurrentMuonTrackable) fHistNGoodClustersForFinalTracks -> Fill(nGoodClusters);
520
521 fNtuFinalCandidates1 -> Fill(Double_t(fRun),
522 Double_t(fEv),
523 Double_t(fCountRealTracksAnalyzedOfEvent),
524 Double_t(nFinalTracks),
525 Double_t(nClustersMC),
526 Double_t(nGoodClusters),
527 ptSpectrometer,
528 thetaSpectrometer,
529 etaSpectrometer);
530
531 fNtuFinalCandidates2 -> Fill(chi2AtPlane[0],
532 chi2AtPlane[1],
533 chi2AtPlane[2],
534 chi2AtPlane[3],
535 chi2AtPlane[4],
536 chi2AtPlane[5],
537 chi2AtPlane[6],
538 chi2AtPlane[7],
539 chi2AtPlane[8],
540 chi2AtPlane[9],
541 chi2AtPlane[10],
542 chi2AtPlane[11],
543 chi2AtPlane[12],
544 chi2AtPlane[13],
545 chi2AtPlane[14]);
546
547 // now comparing the tracks with various criteria, in order to find the best one
548
549 Double_t theVariable = 0.;
550 for (Int_t iCluster=0; iCluster<nMFTClusters; iCluster++) theVariable += chi2AtPlane[iCluster];
551 theVariable /= Double_t(nMFTClusters);
552
553 if (theVariable<theVariable_Best || theVariable_Best<0.) {
554 nGoodClustersBestCandidate = nGoodClusters;
555 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) chi2HistoryForBestCandidate[iPlane] = chi2AtPlane[iPlane];
556 theVariable_Best = theVariable;
557 fTxtTrackFinalChi2 = new TLatex(0.20, 0.52, Form("#chi^{2}_{final} = %3.1f", chi2HistoryForBestCandidate[0]));
558 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
559 fTxtTrackChi2[iPlane] = new TLatex(0.55*fRPlaneMax[fNPlanesMFT-1],
560 0.90*fRPlaneMax[fNPlanesMFT-1],
561 Form("#chi^{2} = %3.1f", chi2AtPlane[iPlane]));
562 }
563 idBestCandidate = iTrack;
564 bestCandidateExists=kTRUE;
565 }
566
567 // ----------------------------------------------------------
568
569 }
570
571 if (nFinalTracks) {
572 FillPlanesWithTrackHistory((AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate));
573 AliMuonForwardTrack *newTrack = (AliMuonForwardTrack*) fCandidateTracks->UncheckedAt(idBestCandidate);
574 new ((*fMuonForwardTracks)[fMuonForwardTracks->GetEntries()]) AliMuonForwardTrack(*newTrack); // AU
575 fNtuFinalBestCandidates1 -> Fill(Double_t(fRun),
576 Double_t(fEv),
577 Double_t(fCountRealTracksAnalyzedOfEvent),
578 Double_t(nFinalTracks),
579 Double_t(nClustersMC),
580 Double_t(nGoodClustersBestCandidate),
581 ptSpectrometer,
582 thetaSpectrometer,
583 etaSpectrometer);
584
585 fNtuFinalBestCandidates2 -> Fill(chi2HistoryForBestCandidate[0],
586 chi2HistoryForBestCandidate[1],
587 chi2HistoryForBestCandidate[2],
588 chi2HistoryForBestCandidate[3],
589 chi2HistoryForBestCandidate[4],
590 chi2HistoryForBestCandidate[5],
591 chi2HistoryForBestCandidate[6],
592 chi2HistoryForBestCandidate[7],
593 chi2HistoryForBestCandidate[8],
594 chi2HistoryForBestCandidate[9],
595 chi2HistoryForBestCandidate[10],
596 chi2HistoryForBestCandidate[11],
597 chi2HistoryForBestCandidate[12],
598 chi2HistoryForBestCandidate[13],
599 chi2HistoryForBestCandidate[14]);
600
601 }
602
603 if (fDrawOption && bestCandidateExists) {
604 fTxtTrackGoodClusters = new TLatex(0.20, 0.59, Form("N_{GoodClusters} = %d", nGoodClustersBestCandidate));
605 DrawPlanes();
606 }
607
608 if (fIsCurrentMuonTrackable) {
609 if (nGoodClustersBestCandidate==5) fHistPtMuonTrackWithGoodMatch -> Fill(ptSpectrometer);
610 else fHistPtMuonTrackWithBadMatch -> Fill(ptSpectrometer);
611 }
612
613 // -------------------------------------------------------------------------------------------
614
615 fCandidateTracks->Clear();
616
617 fCountRealTracksAnalyzedOfEvent++;
618
619 return 5;
620
621}
622
623//===========================================================================================================================================
624
625void AliMuonForwardTrackFinder::FindClusterInPlane(Int_t planeId) {
626
627 printf(">>>> executing AliMuonForwardTrackFinder::FindClusterInPlane(%d)\n", planeId);
628
629 // !!!!!!!!! coordinates and errors on the interaction vertex should be taken from the event itself (ITS) if available
630
631 // propagate track to plane #planeId (both to front and back active sensors)
632 // look for compatible clusters
633 // update TrackParam at found cluster (if any) using Kalman Filter
634
635 AliMUONTrackParam currentParamFront, currentParamBack, currentParamForResearchFront, currentParamForResearchBack;
636
637 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
638 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
639 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
640 currentParamForResearchFront = currentParamFront;
641 currentParamForResearchBack = currentParamBack;
642 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, 0.);
643 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack, 0.);
644 AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchFront, 0., 0., gRandom->Gaus(0,fGaussianBlurZVert), fExtrapOriginTransvError, fExtrapOriginTransvError);
645 AliMUONTrackExtrap::ExtrapToVertex(&currentParamForResearchBack, 0., 0., gRandom->Gaus(0,fGaussianBlurZVert), fExtrapOriginTransvError, fExtrapOriginTransvError);
646 }
647 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
648 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
649 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
650 currentParamForResearchFront = currentParamFront;
651 currentParamForResearchBack = currentParamBack;
652 AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
653 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
654 AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchFront,(fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
655 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
656 AliMUONTrackExtrap::AddMCSEffect(&currentParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
657 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
658 AliMUONTrackExtrap::AddMCSEffect(&currentParamForResearchBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
659 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
660 }
661 // for all planes: extrapolation to the Z of the plane
662 AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
663 AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
664 AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
665 AliMUONTrackExtrap::ExtrapToZCov(&currentParamForResearchBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
666
667 //---------------------------------------------------------------------------------------
668
669 TMatrixD covFront(5,5); covFront = currentParamForResearchFront.GetCovariances();
670 TMatrixD covBack(5,5); covBack = currentParamForResearchBack.GetCovariances();
671
672 Double_t squaredError_X_Front = covFront(0,0);
673 Double_t squaredError_Y_Front = covFront(2,2);
674 Double_t squaredError_X_Back = covBack(0,0);
675 Double_t squaredError_Y_Back = covBack(2,2);
676
677 Double_t corrFact = 1.0;
678
679 Double_t researchRadiusFront = TMath::Sqrt(squaredError_X_Front + squaredError_Y_Front);
680 Double_t researchRadiusBack = TMath::Sqrt(squaredError_X_Back + squaredError_Y_Back);
681 if (planeId==fNPlanesMFT-1 && 0.5*(researchRadiusFront+researchRadiusBack)<fMinResearchRadiusAtLastPlane) {
682 corrFact = fMinResearchRadiusAtLastPlane/(0.5*(researchRadiusFront+researchRadiusBack));
683 }
684 if (fIsCurrentMuonTrackable) fHistResearchRadius[planeId] -> Fill(corrFact*0.5*(researchRadiusFront+researchRadiusBack));
685
686 Double_t position_X_Front = currentParamForResearchFront.GetNonBendingCoor();
687 Double_t position_Y_Front = currentParamForResearchFront.GetBendingCoor();
688 Double_t position_X_Back = currentParamForResearchBack.GetNonBendingCoor();
689 Double_t position_Y_Back = currentParamForResearchBack.GetBendingCoor();
690 Double_t radialPositionOfTrackFront = TMath::Sqrt(position_X_Front*position_X_Front + position_Y_Front*position_Y_Front);
691 Double_t radialPositionOfTrackBack = TMath::Sqrt(position_X_Back*position_X_Back + position_Y_Back*position_Y_Back);
692
693 //---------------------------------------------------------------------------------------
694
695 Double_t chi2cut = 2.*fSigmaClusterCut*fSigmaClusterCut; // depends on the number of variables (here, 2)
696
697 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
698
699 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
700 printf("There are %3d clusters in plane %02d FRONT\n", nClustersFront, planeId);
701
702 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
703
704 Bool_t isGoodChi2 = kFALSE;
705
706 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->At(iCluster);
707 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchFront, cluster); // describes the compatibility between the track and the cluster
708 if (chi2<chi2cut) isGoodChi2 = kTRUE;
709
710 Double_t radialPositionOfClusterFront = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
711 if (planeId == fNPlanesMFT-1) {
712 if (TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront)<fDistanceFromBestClusterAndTrackAtLastPlane ||
713 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
714 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
715 }
716 if (IsCorrectMatch(cluster)) {
717 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackFront-radialPositionOfClusterFront);
718 }
719 }
720
721 if (fIsCurrentMuonTrackable) {
722 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
723 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
724 }
725
726 if (isGoodChi2) {
727 printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
728 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
729 newTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
730 newTrack->SetPlaneExists(planeId);
731 printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
732 if (fIsCurrentMuonTrackable) {
733 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
734 printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
735 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
736 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
737 if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
738 else fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
739 }
740 fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
741 }
742 else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
743
744 }
745
746 // Analyizing the clusters: BACK ACTIVE ELEMENTS
747
748 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
749 printf("There are %3d clusters in plane %02d BACK\n", nClustersBack, planeId);
750
751 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
752
753 Bool_t isGoodChi2 = kFALSE;
754
755 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->At(iCluster);
756 Double_t chi2 = (1./(corrFact*corrFact)) * TryOneCluster(currentParamForResearchBack, cluster); // describes the compatibility between the track and the cluster
757 if (chi2<chi2cut) isGoodChi2 = kTRUE;
758
759 Double_t radialPositionOfClusterBack = TMath::Sqrt(cluster->GetX()*cluster->GetX() + cluster->GetY()*cluster->GetY());
760 if (planeId == fNPlanesMFT-1) {
761 if (TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack)<fDistanceFromBestClusterAndTrackAtLastPlane ||
762 fDistanceFromBestClusterAndTrackAtLastPlane<0.) {
763 fDistanceFromBestClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
764 }
765 if (IsCorrectMatch(cluster)) {
766 fDistanceFromGoodClusterAndTrackAtLastPlane = TMath::Abs(radialPositionOfTrackBack-radialPositionOfClusterBack);
767 }
768 }
769
770 if (fIsCurrentMuonTrackable) {
771 if (IsCorrectMatch(cluster)) fHistChi2Cluster_GoodCluster[planeId]->Fill(chi2/2.); // chi2/ndf
772 else fHistChi2Cluster_BadCluster[planeId] ->Fill(chi2/2.); // chi2/ndf
773 }
774
775 if (isGoodChi2) {
776 printf("accepting cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
777 AliMuonForwardTrack *newTrack = new ((*fCandidateTracks)[fCandidateTracks->GetEntriesFast()]) AliMuonForwardTrack(*fCurrentTrack);
778 newTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
779 newTrack->SetPlaneExists(planeId);
780 printf("current muon is trackable: %d\n", fIsCurrentMuonTrackable);
781 if (fIsCurrentMuonTrackable) {
782 Double_t newGlobalChi2 = ((AliMUONTrackParam*) newTrack->GetTrackParamAtCluster()->First())->GetTrackChi2();
783 printf("new chi2 = %f (= %f)\n", newGlobalChi2, newTrack->GetMFTCluster(0)->GetTrackChi2());
784 Int_t nClustersGlobalTrack = newTrack->GetNMUONClusters() + newTrack->GetNMFTClusters(); // Muon Spectrometer clusters + clusters in the Vertex Telescope
785 Int_t ndfGlobalTrack = GetNDF(nClustersGlobalTrack);
786 if (IsCorrectMatch(cluster)) fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[planeId]->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
787 else fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[planeId] ->Fill(newGlobalChi2/Double_t(ndfGlobalTrack));
788 }
789 fGrMFTPlane[planeId][kClustersGoodChi2] -> SetPoint(fGrMFTPlane[planeId][kClustersGoodChi2]->GetN(), cluster->GetX(), cluster->GetY());
790 }
791 else printf("discarding cluster: chi2=%f (cut = %f)\n", chi2, chi2cut);
792
793 }
794
795 //---------------------------------------------------------------------------------------------
796
797 if (planeId == fNPlanesMFT-1) {
798 if (fIsCurrentMuonTrackable && fDistanceFromGoodClusterAndTrackAtLastPlane>0.) {
799 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Fill(TMath::Abs(fDistanceFromBestClusterAndTrackAtLastPlane-
800 fDistanceFromGoodClusterAndTrackAtLastPlane));
801 fHistDistanceGoodClusterFromTrackAtLastPlane -> Fill(fDistanceFromGoodClusterAndTrackAtLastPlane);
802 }
803 }
804
805}
806
807//==========================================================================================================================================
808
809void AliMuonForwardTrackFinder::AttachGoodClusterInPlane(Int_t planeId) {
810
811 printf(">>>> executing AliMuonForwardTrackFinder::AttachGoodClusterInPlane(%d)\n", planeId);
812
813 AliMUONTrackParam currentParamFront, currentParamBack;
814
815 if (planeId == fNPlanesMFT-1) { // last plane of the telecope
816 currentParamFront = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
817 currentParamBack = (*((AliMUONTrackParam*)(fMuonTrackReco->GetTrackParamAtCluster()->First())));
818 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamFront, 0.);
819 AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&currentParamBack, 0.);
820 }
821 else { // MFT planes others than the last one: mult. scattering correction because of the upstream MFT planes is performed
822 printf("fCurrentTrack = %p\n", fCurrentTrack);
823 currentParamFront = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
824 currentParamBack = (*((AliMUONTrackParam*)(fCurrentTrack->GetTrackParamAtCluster()->First())));
825 AliMUONTrackExtrap::AddMCSEffect(&currentParamFront, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
826 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeFront())/radLengthSi,-1.);
827 AliMUONTrackExtrap::AddMCSEffect(&currentParamBack, (fSegmentation->GetPlane(planeId+1)->GetEquivalentSilicon()+
828 fSegmentation->GetPlane(planeId)->GetEquivalentSiliconBeforeBack())/radLengthSi,-1.);
829 }
830 // for all planes: linear extrapolation to the Z of the plane
831 AliMUONTrackExtrap::ExtrapToZCov(&currentParamFront, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveFront());
832 AliMUONTrackExtrap::ExtrapToZCov(&currentParamBack, -1.*fSegmentation->GetPlane(planeId)->GetZCenterActiveBack());
833
834 Bool_t goodClusterFound = kFALSE;
835
836 // Analyizing the clusters: FRONT ACTIVE ELEMENTS
837
838 Int_t nClustersFront = fMFTClusterArrayFront[planeId]->GetEntries();
839
840 printf("nClustersFront = %d\n", nClustersFront);
841 for (Int_t iCluster=0; iCluster<nClustersFront; iCluster++) {
842 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayFront[planeId]->UncheckedAt(iCluster);
843 printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersFront, cluster, fCurrentTrack);
844 if (IsCorrectMatch(cluster)) {
845 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamFront, *cluster); // creating new track param and attaching the cluster
846 fCurrentTrack->SetPlaneExists(planeId);
847 goodClusterFound = kTRUE;
848 break;
849 }
850 }
851
852 if (goodClusterFound) return;
853
854 // Analyizing the clusters: BACK ACTIVE ELEMENTS
855
856 Int_t nClustersBack = fMFTClusterArrayBack[planeId]->GetEntries();
857
858 printf("nClustersBack = %d\n", nClustersBack);
859 for (Int_t iCluster=0; iCluster<nClustersBack; iCluster++) {
860 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArrayBack[planeId]->UncheckedAt(iCluster);
861 printf("checking cluster %02d of %02d: cluter=%p, fCurrentTrack=%p\n", iCluster, nClustersBack, cluster, fCurrentTrack);
862 if (IsCorrectMatch(cluster)) {
863 fCurrentTrack->AddTrackParamAtMFTCluster(currentParamBack, *cluster); // creating new track param and attaching the cluster
864 fCurrentTrack->SetPlaneExists(planeId);
865 goodClusterFound = kTRUE;
866 break;
867 }
868 }
869
870}
871
872//==========================================================================================================================================
873
874void AliMuonForwardTrackFinder::CheckCurrentMuonTrackable() {
875
876 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
877 fIsGoodClusterInPlane[iPlane] = kFALSE;
878 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
879 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
880 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
881 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
882 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
883 fIsGoodClusterInPlane[iPlane] = kTRUE;
884 break;
885 }
886 }
887 }
888 }
889
890 fIsCurrentMuonTrackable = kTRUE;
891 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) fIsCurrentMuonTrackable = (fIsCurrentMuonTrackable&&fIsGoodClusterInPlane[iPlane]);
892
893}
894
895//==========================================================================================================================================
896
897void AliMuonForwardTrackFinder::FillPlanesWithTrackHistory(AliMuonForwardTrack *track) {
898
899 // recover track parameters on each planes and look for the corresponding clusters
900
901 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
902
903 AliMFTCluster *trackCluster = (AliMFTCluster *) track->GetMFTCluster(iPlane);
904
905 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetPoint(fGrMFTPlane[iPlane][kClusterOfTrack]->GetN(), trackCluster->GetX(), trackCluster->GetY());
906
907 Int_t nClusters = fMFTClusterArray[iPlane]->GetEntriesFast();
908
909 for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
910
911 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->UncheckedAt(iCluster);
912
913 fGrMFTPlane[iPlane][kAllClusters] -> SetPoint(fGrMFTPlane[iPlane][kAllClusters]->GetN(), cluster->GetX(), cluster->GetY());
914
915 if (IsCorrectMatch(cluster)) {
916 fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetPoint(fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN(), cluster->GetX(), cluster->GetY());
917 }
918
919 }
920
921 }
922
923}
924
925//======================================================================================================================================
926
927Bool_t AliMuonForwardTrackFinder::IsCorrectMatch(AliMFTCluster *cluster) {
928
929 Bool_t result = kFALSE;
930
931 // check if the cluster belongs to the correct MC track
932
933 for (Int_t iTrack=0; iTrack<cluster->GetNMCTracks(); iTrack++) {
934 if (cluster->GetMCLabel(iTrack)==fLabelMC) {
935 result = kTRUE;
936 break;
937 }
938 }
939
940 printf("returning %d\n", result);
941
942 return result;
943
944}
945
946//======================================================================================================================================
947
948Double_t AliMuonForwardTrackFinder::TryOneCluster(const AliMUONTrackParam &trackParam, AliMFTCluster *cluster) {
949
950 // Test the compatibility between the track and the cluster (using trackParam's covariance matrix):
951 // return the corresponding Chi2
952 // assume the track parameters are given at the Z of the cluster
953
954 // Set differences between trackParam and cluster in the bending and non bending directions
955 Double_t dX = cluster->GetX() - trackParam.GetNonBendingCoor();
956 Double_t dY = cluster->GetY() - trackParam.GetBendingCoor();
957 printf("dX = %f, dY = %f\n", dX, dY);
958
959 // Calculate errors and covariances
960 const TMatrixD& kParamCov = trackParam.GetCovariances();
961 Double_t sigmaX2 = kParamCov(0,0) + cluster->GetErrX2();
962 Double_t sigmaY2 = kParamCov(2,2) + cluster->GetErrY2();
963 printf("dX2 = %f, dY2 = %f\n", sigmaX2, sigmaY2);
964 Double_t covXY = kParamCov(0,2);
965 Double_t det = sigmaX2 * sigmaY2 - covXY * covXY;
966
967 // Compute chi2
968 if (det==0.) return 1.e10;
969 return (dX*dX*sigmaY2 + dY*dY*sigmaX2 - 2.*dX*dY*covXY) / det;
970
971}
972
973//=========================================================================================================================================
974
975void AliMuonForwardTrackFinder::SeparateFrontBackClusters() {
976
977 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
978 fMFTClusterArrayFront[iPlane]->Clear();
979 fMFTClusterArrayBack[iPlane] ->Clear();
980 for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
981 AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
982 if (TMath::Abs(cluster->GetZ())<TMath::Abs(fSegmentation->GetPlane(iPlane)->GetZCenter())) {
983 new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
984 }
985 else {
986 new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
987 }
988 }
989 }
990
991}
992
993//=========================================================================================================================================
994
995Int_t AliMuonForwardTrackFinder::GetNDF(Int_t nClusters) {
996
997 // the same definition as in AliMUONTrack is implemented, since here we just add more clusters to the Muon track
998
999 Int_t ndf = 2 * nClusters - 5;
1000 return (ndf > 0) ? ndf : 0;
1001
1002}
1003
1004//============================================================================================================================================
1005
1006void AliMuonForwardTrackFinder::BookHistos() {
1007
1008 const Int_t nMaxNewTracks[] = {150, 200, 250, 600, 1000};
1009 const Double_t radiusPlane[] = {0.001, 0.010, 0.100, 5.0, 5.0};
1010
1011 fHistPtSpectrometer = new TH1D("hPtSpectrometer", "p_{T} as given by the Muon Spectrometer", 200, 0, 20.);
1012
1013 fHistPtMuonTrackWithGoodMatch = new TH1D("fHistPtMuonTrackWithGoodMatch", "p_{T} of muon track with good match", 200, 0, 20.);
1014 fHistPtMuonTrackWithBadMatch = new TH1D("fHistPtMuonTrackWithBadMatch", "p_{T} of muon track with bad match", 200, 0, 20.);
1015
1016 fHistRadiusEndOfAbsorber = new TH1D("hRadiusEndOfAbsorber", "Track radial distance at the end of the absorber", 1000, 0, 100.);
1017
1018 fHistNGoodClustersForFinalTracks = new TH1D("hNGoodClustersForFinalTracks", "Number of Good Clusters per Final Track", 20, -0.25, 9.75);
1019
1020 fHistDistanceGoodClusterFromTrackAtLastPlane = new TH1D("hDistanceGoodClusterFromTrackAtLastPlane",
1021 "Distance of MC Good Cluster from Track in last MFT plane", 200, 0., 2.);
1022
1023 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane =
1024 new TH1D("hDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane",
1025 "Good Cluster distance from track - Best Cluster distance from track in last MFT plane", 200, 0., 2.);
1026
1027 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1028
1029 fHistNTracksAfterExtrapolation[iPlane] = new TH1D(Form("hNTracksAfterExtrapolation_pl%02d", iPlane),
1030 Form("Number of Candidates after analysis of MFT plane %02d", iPlane),
1031 nMaxNewTracks[iPlane], -0.5, nMaxNewTracks[iPlane]-0.5);
1032
1033 fHistResearchRadius[iPlane] = new TH1D(Form("hResearchRadius_pl%02d", iPlane),
1034 Form("Research Radius for candidate clusters in MFT plane %02d", iPlane),
1035 1000, 0., radiusPlane[iPlane]);
1036
1037 fHistChi2Cluster_GoodCluster[iPlane] = new TH1D(Form("hChi2Cluster_GoodCluster_pl%02d", iPlane),
1038 Form("#chi^{2}_{clust} for Good clusters in MFT plane %02d", iPlane),
1039 100, 0., 15.);
1040
1041 fHistChi2Cluster_BadCluster[iPlane] = new TH1D(Form("hChi2Cluster_BadCluster_pl%02d", iPlane),
1042 Form("#chi^{2}_{clust} for Bad clusters in MFT plane %02d", iPlane),
1043 100, 0., 15.);
1044
1045 fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1046 Form("#chi^{2}/ndf at plane %d for GOOD candidates of trackable muons",iPlane),
1047 100, 0., 15.);
1048
1049 fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] = new TH1D(Form("fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons_pl%02d", iPlane),
1050 Form("#chi^{2}/ndf at plane %d for BAD candidates of trackable muons",iPlane),
1051 100, 0., 15.);
1052
1053 }
1054
1055 //------------------------------------------
1056
1057 fHistPtSpectrometer -> Sumw2();
1058 fHistPtMuonTrackWithGoodMatch -> Sumw2();
1059 fHistPtMuonTrackWithBadMatch -> Sumw2();
1060 fHistRadiusEndOfAbsorber -> Sumw2();
1061 fHistNGoodClustersForFinalTracks -> Sumw2();
1062
1063 fHistDistanceGoodClusterFromTrackAtLastPlane -> Sumw2();
1064 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Sumw2();
1065
1066 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1067
1068 fHistNTracksAfterExtrapolation[iPlane] -> Sumw2();
1069 fHistResearchRadius[iPlane] -> Sumw2();
1070
1071 fHistChi2Cluster_GoodCluster[iPlane] -> Sumw2();
1072 fHistChi2Cluster_BadCluster[iPlane] -> Sumw2();
1073
1074 fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1075 fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Sumw2();
1076
1077 }
1078
1079 fNtuFinalCandidates1 = new TNtuple("ntuFinalCandidates1", "Final Candidates (ALL)", "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
1080
1081 fNtuFinalBestCandidates1 = new TNtuple("ntuFinalBestCandidates1", "Final Best Candidates", "run:event:muonTrack:nFinalCandidates:nClustersMC:nGoodClusters:ptSpectrometer:thetaSpectrometer:etaSpectrometer");
1082
1083 fNtuFinalCandidates2 = new TNtuple("ntuFinalCandidates2", "Final Candidates (ALL)", "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
1084
1085 fNtuFinalBestCandidates2 = new TNtuple("ntuFinalBestCandidates2", "Final Best Candidates", "chi2AtPlane0:chi2AtPlane1:chi2AtPlane2:chi2AtPlane3:chi2AtPlane4:chi2AtPlane5:chi2AtPlane6:chi2AtPlane7:chi2AtPlane8:chi2AtPlane9:chi2AtPlane10:chi2AtPlane11:chi2AtPlane12:chi2AtPlane13:chi2AtPlane14");
1086
1087}
1088
1089//============================================================================================================================================
1090
1091void AliMuonForwardTrackFinder::SetTitleHistos() {
1092
1093 fHistPtSpectrometer -> SetXTitle("p_{T} [GeV/c]");
1094 fHistPtMuonTrackWithGoodMatch -> SetXTitle("p_{T} [GeV/c]");
1095 fHistPtMuonTrackWithBadMatch -> SetXTitle("p_{T} [GeV/c]");
1096 fHistRadiusEndOfAbsorber -> SetXTitle("R_{abs} [cm]");
1097 fHistNGoodClustersForFinalTracks -> SetXTitle("N_{GoodClusters}");
1098
1099 fHistDistanceGoodClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1100 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> SetXTitle("Distance [cm]");
1101
1102
1103 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1104
1105 fHistNTracksAfterExtrapolation[iPlane] -> SetXTitle("N_{tracks}");
1106 fHistResearchRadius[iPlane] -> SetXTitle("Research Radius [cm]");
1107
1108 fHistChi2Cluster_GoodCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1109 fHistChi2Cluster_BadCluster[iPlane] -> SetXTitle("#chi^{2}/ndf");
1110
1111 fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1112 fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> SetXTitle("#chi^{2}/ndf");
1113
1114 }
1115
1116}
1117
1118//===========================================================================================================================================
1119
1120void AliMuonForwardTrackFinder::BookPlanes() {
1121
1122 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1123 fGrMFTPlane[iPlane][kAllClusters] = new TGraph();
1124 fGrMFTPlane[iPlane][kAllClusters] -> SetName(Form("fGrMFTPlane_%02d_AllClusters",iPlane));
1125 fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerStyle(20);
1126 // fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.5);
1127 // fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.3);
1128 fGrMFTPlane[iPlane][kAllClusters] -> SetMarkerSize(0.2);
1129 }
1130
1131 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1132 fGrMFTPlane[iPlane][kClustersGoodChi2] = new TGraph();
1133 fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetName(Form("fGrMFTPlane_%02d_ClustersGoodChi2",iPlane));
1134 fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerStyle(20);
1135 // fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.8);
1136 // fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.4);
1137 fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerSize(0.3);
1138 fGrMFTPlane[iPlane][kClustersGoodChi2] -> SetMarkerColor(kBlue);
1139 }
1140
1141 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1142 fGrMFTPlane[iPlane][kClusterOfTrack] = new TGraph();
1143 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetName(Form("fGrMFTPlane_%02d_ClustersOfTrack",iPlane));
1144 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerStyle(25);
1145 // fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(1.2);
1146 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerSize(0.9);
1147 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetMarkerColor(kRed);
1148 fGrMFTPlane[iPlane][kClusterOfTrack] -> SetTitle(Form("Plane %d (%3.1f cm)", iPlane, fZPlane[iPlane]));
1149 }
1150
1151 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1152 fGrMFTPlane[iPlane][kClusterCorrectMC] = new TGraph();
1153 fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetName(Form("fGrMFTPlane_%02d_ClustersCorrectMC",iPlane));
1154 fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerStyle(20);
1155 // fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.8);
1156 fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerSize(0.5);
1157 fGrMFTPlane[iPlane][kClusterCorrectMC] -> SetMarkerColor(kGreen);
1158 }
1159
1160 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1161 fCircleExt[iPlane] = new TEllipse(0., 0., fRPlaneMax[iPlane], fRPlaneMax[iPlane]);
1162 fCircleInt[iPlane] = new TEllipse(0., 0., fRPlaneMin[iPlane], fRPlaneMin[iPlane]);
1163 }
1164
1165 fTxtDummy = new TLatex(0.10, 0.67, "Best Candidate:");
1166
1167 //---------------------------------------------------
1168
1169 fMrkAllClust = new TMarker(0.10, 0.32, 20);
1170 fMrkAllClust -> SetMarkerSize(0.5);
1171
1172 fMrkClustGoodChi2 = new TMarker(0.10, 0.26, 20);
1173 fMrkClustGoodChi2 -> SetMarkerSize(0.8);
1174 fMrkClustGoodChi2 -> SetMarkerColor(kBlue);
1175
1176 fMrkClustMC = new TMarker(0.10, 0.20, 20);
1177 fMrkClustMC -> SetMarkerSize(0.8);
1178 fMrkClustMC -> SetMarkerColor(kGreen);
1179
1180 fMrkClustOfTrack = new TMarker(0.10, 0.14, 25);
1181 fMrkClustOfTrack -> SetMarkerSize(1.2);
1182 fMrkClustOfTrack -> SetMarkerColor(kRed);
1183
1184 fTxtAllClust = new TLatex(0.15, 0.30, "All Clusters");
1185 fTxtAllClust -> SetTextSize(0.040);
1186
1187 fTxtClustGoodChi2 = new TLatex(0.15, 0.24, "Clusters involved in the research");
1188 fTxtClustGoodChi2 -> SetTextSize(0.040);
1189
1190 fTxtClustMC = new TLatex(0.15, 0.18, "MC good clusters");
1191 fTxtClustMC -> SetTextSize(0.040);
1192
1193 fTxtClustOfTrack = new TLatex(0.15, 0.12, "Clusters of the best candidate");
1194 fTxtClustOfTrack -> SetTextSize(0.040);
1195
1196}
1197
1198//===========================================================================================================================================
1199
1200void AliMuonForwardTrackFinder::ResetPlanes() {
1201
1202 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1203 for (Int_t iGr=0; iGr<4; iGr++) {
1204 Int_t nOldClusters = fGrMFTPlane[iPlane][iGr]->GetN();
1205 for (Int_t iPoint=nOldClusters-1; iPoint>=0; iPoint--) fGrMFTPlane[iPlane][iGr]->RemovePoint(iPoint);
1206 }
1207 }
1208
1209}
1210
1211//===========================================================================================================================================
1212
1213void AliMuonForwardTrackFinder::PrintParticleHistory() {
1214
1215 TString history = "";
1216
1217 TParticle *part = fStack->Particle(fLabelMC);
1218
1219 if (part->GetFirstMother() != -1) {
1220 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1221 if (partMother->GetFirstMother() != -1) history += "... #rightarrow ";
1222 Char_t newName[100];
1223 PDGNameConverter(partMother->GetName(), newName);
1224 history += Form("%s #rightarrow ", newName);
1225 }
1226 Char_t newName[100];
1227 PDGNameConverter(part->GetName(), newName);
1228 history += Form("%s at z = %5.1f cm", newName, part->Vz());
1229
1230 // printf("%s", history.Data());
1231
1232 fTxtMuonHistory = new TLatex(0.10, 0.86, history.Data());
1233
1234}
1235
1236//===========================================================================================================================================
1237
1238Bool_t AliMuonForwardTrackFinder::IsMother(Char_t *nameMother) {
1239
1240 Bool_t result = kFALSE;
1241
1242 TParticle *part = fStack->Particle(fLabelMC);
1243
1244 if (part->GetFirstMother() != -1) {
1245 TParticle *partMother = fStack->Particle(part->GetFirstMother());
1246 if (!strcmp(partMother->GetName(), nameMother)) result=kTRUE;
1247 }
1248
1249 return result;
1250
1251}
1252
1253//===========================================================================================================================================
1254
1255void AliMuonForwardTrackFinder::DrawPlanes() {
1256
1257 fCanvas -> Clear();
1258 if (fNPlanesMFT <= 5) fCanvas -> Divide(3,2);
1259 else if (fNPlanesMFT <= 11) fCanvas -> Divide(4,3);
1260 else if (fNPlanesMFT <= 19) fCanvas -> Divide(5,4);
1261
1262 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1263
1264 fCanvas->cd(fNPlanesMFT-iPlane+1);
1265
1266 fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetLimits(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1267 fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetRangeUser(-1.1*fRPlaneMax[fNPlanesMFT-1], +1.1*fRPlaneMax[fNPlanesMFT-1]);
1268 fGrMFTPlane[iPlane][kClusterOfTrack] -> GetXaxis() -> SetTitle("X [cm]");
1269 fGrMFTPlane[iPlane][kClusterOfTrack] -> GetYaxis() -> SetTitle("Y [cm]");
1270 fGrMFTPlane[iPlane][kClusterOfTrack] -> Draw("ap");
1271
1272 fCircleExt[iPlane] -> Draw("same");
1273 fCircleInt[iPlane] -> Draw("same");
1274
1275 if (fGrMFTPlane[iPlane][kAllClusters]->GetN()) fGrMFTPlane[iPlane][kAllClusters] -> Draw("psame");
1276 if (fGrMFTPlane[iPlane][kClustersGoodChi2]->GetN()) fGrMFTPlane[iPlane][kClustersGoodChi2] -> Draw("psame");
1277 if (fGrMFTPlane[iPlane][kClusterOfTrack]->GetN()) fGrMFTPlane[iPlane][kClusterOfTrack] -> Draw("psame");
1278 if (fGrMFTPlane[iPlane][kClusterCorrectMC]->GetN()) fGrMFTPlane[iPlane][kClusterCorrectMC] -> Draw("psame");
1279
1280 fTxtTrackChi2[iPlane] -> Draw("same");
1281
1282 }
1283
1284 fCanvas -> cd(1);
1285 fTxtMuonHistory -> Draw();
1286 fTxtDummy -> Draw("same");
1287 fTxtTrackGoodClusters -> Draw("same");
1288 fTxtTrackFinalChi2 -> Draw("same");
1289 fTxtFinalCandidates -> Draw("same");
1290
1291 fMrkAllClust -> Draw("same");
1292 fMrkClustGoodChi2 -> Draw("same");
1293 fMrkClustMC -> Draw("same");
1294 fMrkClustOfTrack -> Draw("same");
1295
1296 fTxtAllClust -> Draw("same");
1297 fTxtClustGoodChi2 -> Draw("same");
1298 fTxtClustMC -> Draw("same");
1299 fTxtClustOfTrack -> Draw("same");
1300
1301 // fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1302 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1303 if (IsMother("phi")) {
1304 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1305 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.phi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1306 }
1307 if (IsMother("J/psi")) {
1308 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.gif", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1309 fCanvas -> SaveAs(Form("%s/figures/eventDisplay/run%d_event%d_track%d.jPsi.eps", fOutDir.Data(), fRun, fEv, fCountRealTracksAnalyzedOfEvent));
1310 }
1311
1312}
1313
1314//===========================================================================================================================================
1315
1316void AliMuonForwardTrackFinder::Terminate() {
1317
1318 printf("\n");
1319 printf("---------------------------------------------------------------------------------------------------------------\n");
1320 printf("%8d tracks analyzed\n", fCountRealTracksAnalyzed);
1321 printf("%8d tracks with MC ref\n", fCountRealTracksWithRefMC);
1322 printf("%8d tracks with MC ref & trigger match\n", fCountRealTracksWithRefMC_andTrigger);
1323 printf("%8d tracks analyzed with final candidates\n", fCountRealTracksAnalyzedWithFinalCandidates);
1324// printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c\n", fCountRealTracksWithRefMC_andTrigger_andGoodPt, fLowPtCut);
1325// printf("%8d tracks with MC ref & trigger match & pt>%3.1f GeV/c & correct R_abs\n", fCountRealTracksWithRefMC_andTrigger_andGoodPt_andGoodTheta, fLowPtCut);
1326 printf("---------------------------------------------------------------------------------------------------------------\n");
1327
1328 WriteOutputTree();
1329 WriteHistos();
1330
1331}
1332
1333//==========================================================================================================================================
1334
1335void AliMuonForwardTrackFinder::FillOutputTree() {
1336
1337 if (!fMuonForwardTracks || !fOutputEventTree) return;
1338
1339 AliDebug(1, Form("Filling output tree %p with %p having %d entries whose 1st entry is %p",
1340 fOutputEventTree, fMuonForwardTracks, fMuonForwardTracks->GetEntries(), fMuonForwardTracks->At(0)));
1341
1342 fOutputEventTree->Fill();
1343 AliDebug(1, Form("\n\nFilled Tree: nEvents = %d!!!!\n\n", Int_t(fOutputEventTree->GetEntries())));
1344
1345}
1346
1347//==========================================================================================================================================
1348
1349void AliMuonForwardTrackFinder::WriteOutputTree() {
1350
1351 if (!fOutputEventTree || !fOutputTreeFile) return;
1352
1353 fOutputTreeFile -> cd();
1354
1355 fOutputEventTree -> Write();
1356 fOutputTreeFile -> Close();
1357
1358}
1359
1360//==========================================================================================================================================
1361
1362void AliMuonForwardTrackFinder::WriteHistos() {
1363
1364 TFile *fileOut = new TFile(Form("%s/MuonGlobalTracking.QA.run%d.root", fOutDir.Data(), fRun), "recreate");
1365
1366 fHistPtSpectrometer -> Write();
1367 fHistPtMuonTrackWithGoodMatch -> Write();
1368 fHistPtMuonTrackWithBadMatch -> Write();
1369 fHistRadiusEndOfAbsorber -> Write();
1370 fHistNGoodClustersForFinalTracks -> Write();
1371
1372 fHistDistanceGoodClusterFromTrackAtLastPlane -> Write();
1373 fHistDistanceGoodClusterFromTrackMinusDistanceBestClusterFromTrackAtLastPlane -> Write();
1374
1375 for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
1376
1377 fHistNTracksAfterExtrapolation[iPlane] -> Write();
1378 fHistResearchRadius[iPlane] -> Write();
1379
1380 fHistChi2Cluster_GoodCluster[iPlane] -> Write();
1381 fHistChi2Cluster_BadCluster[iPlane] -> Write();
1382
1383 fHistChi2AtPlaneFor_GOOD_CandidatesOfTrackableMuons[iPlane] -> Write();
1384 fHistChi2AtPlaneFor_BAD_CandidatesOfTrackableMuons[iPlane] -> Write();
1385
1386 }
1387
1388 fNtuFinalCandidates1 -> Write();
1389 fNtuFinalBestCandidates1 -> Write();
1390 fNtuFinalCandidates2 -> Write();
1391 fNtuFinalBestCandidates2 -> Write();
1392
1393 fileOut -> Close();
1394
1395}
1396
1397//===========================================================================================================================================
1398
1399void AliMuonForwardTrackFinder::PDGNameConverter(const Char_t *nameIn, Char_t *nameOut) {
1400
1401 if (!strcmp(nameIn, "mu+")) sprintf(nameOut, "#mu^{+}");
1402 else if (!strcmp(nameIn, "mu-")) sprintf(nameOut, "#mu^{-}");
1403 else if (!strcmp(nameIn, "pi+")) sprintf(nameOut, "#pi^{+}");
1404 else if (!strcmp(nameIn, "pi-")) sprintf(nameOut, "#pi^{-}");
1405 else if (!strcmp(nameIn, "K+")) sprintf(nameOut, "K^{+}");
1406 else if (!strcmp(nameIn, "K-")) sprintf(nameOut, "K^{-}");
1407 else if (!strcmp(nameIn, "K*+")) sprintf(nameOut, "K^{*+}");
1408 else if (!strcmp(nameIn, "K*-")) sprintf(nameOut, "K^{*-}");
1409 else if (!strcmp(nameIn, "K_S0")) sprintf(nameOut, "K_{S}^{0}");
1410 else if (!strcmp(nameIn, "K_L0")) sprintf(nameOut, "K_{L}^{0}");
1411 else if (!strcmp(nameIn, "K0")) sprintf(nameOut, "K^{0}");
1412 else if (!strcmp(nameIn, "K0_bar")) sprintf(nameOut, "#bar{K}^{0}");
1413 else if (!strcmp(nameIn, "K*0")) sprintf(nameOut, "K^{*0}");
1414 else if (!strcmp(nameIn, "K*0_bar")) sprintf(nameOut, "#bar{K}^{*0}");
1415 else if (!strcmp(nameIn, "rho0")) sprintf(nameOut, "#rho^{0}");
1416 else if (!strcmp(nameIn, "rho+")) sprintf(nameOut, "#rho^{+}");
1417 else if (!strcmp(nameIn, "rho-")) sprintf(nameOut, "#rho^{-}");
1418 else if (!strcmp(nameIn, "omega")) sprintf(nameOut, "#omega");
1419 else if (!strcmp(nameIn, "eta'")) sprintf(nameOut, "#eta'");
1420 else if (!strcmp(nameIn, "phi")) sprintf(nameOut, "#phi");
1421
1422 else if (!strcmp(nameIn, "D-")) sprintf(nameOut, "D^{-}");
1423 else if (!strcmp(nameIn, "D+")) sprintf(nameOut, "D^{+}");
1424 else if (!strcmp(nameIn, "D0")) sprintf(nameOut, "D^{0}");
1425 else if (!strcmp(nameIn, "D0_bar")) sprintf(nameOut, "#bar{D}^{0}");
1426 else if (!strcmp(nameIn, "D*-")) sprintf(nameOut, "D^{*-}");
1427 else if (!strcmp(nameIn, "D*+")) sprintf(nameOut, "D^{*+}");
1428 else if (!strcmp(nameIn, "D_s+")) sprintf(nameOut, "D_{s}^{+}");
1429 else if (!strcmp(nameIn, "D*_s+")) sprintf(nameOut, "D_{s}^{*+}");
1430
1431 else if (!strcmp(nameIn, "B-")) sprintf(nameOut, "B^{-}");
1432 else if (!strcmp(nameIn, "B+")) sprintf(nameOut, "B^{+}");
1433 else if (!strcmp(nameIn, "B_s0_bar")) sprintf(nameOut, "#bar{B}_{s}^{0}");
1434
1435 else if (!strcmp(nameIn, "antiproton")) sprintf(nameOut, "#bar{p}");
1436 else if (!strcmp(nameIn, "proton")) sprintf(nameOut, "p");
1437 else if (!strcmp(nameIn, "neutron")) sprintf(nameOut, "n");
1438 else if (!strcmp(nameIn, "Sigma+")) sprintf(nameOut, "#Sigma^{+}");
1439 else if (!strcmp(nameIn, "Delta+")) sprintf(nameOut, "#Delta{+}");
1440 else if (!strcmp(nameIn, "Delta--")) sprintf(nameOut, "#Delta{--}");
1441 else if (!strcmp(nameIn, "Lambda0")) sprintf(nameOut, "#Lambda_0");
1442 else if (!strcmp(nameIn, "Lambda0_bar")) sprintf(nameOut, "#bar{Lambda}_0");
1443
1444 else sprintf(nameOut, "%s", nameIn);
1445
1446}
1447
1448//===========================================================================================================================================
1449
1450void AliMuonForwardTrackFinder::SetDraw(Bool_t drawOption) {
1451
1452 fDrawOption = drawOption;
1453
1454 if (!fCanvas) {
1455 fCanvas = new TCanvas("tracking", "tracking", 1200, 800);
1456 fCanvas -> Divide(3,2);
1457 }
1458
1459}
1460
1461//===========================================================================================================================================
1462
1463Bool_t AliMuonForwardTrackFinder::InitGRP() {
1464
1465 //------------------------------------
1466 // Initialization of the GRP entry
1467 //------------------------------------
1468
1469 AliCDBEntry* entry = AliCDBManager::Instance()->Get("GRP/GRP/Data");
1470
1471 if (entry) {
1472
1473 TMap* m = dynamic_cast<TMap*>(entry->GetObject()); // old GRP entry
1474
1475 if (m) {
1476 AliInfo("Found a TMap in GRP/GRP/Data, converting it into an AliGRPObject");
1477 m->Print();
1478 fGRPData = new AliGRPObject();
1479 fGRPData->ReadValuesFromMap(m);
1480 }
1481
1482 else {
1483 AliInfo("Found an AliGRPObject in GRP/GRP/Data, reading it");
1484 fGRPData = dynamic_cast<AliGRPObject*>(entry->GetObject()); // new GRP entry
1485 entry->SetOwner(0);
1486 }
1487
1488 // FIX ME: The unloading of GRP entry is temporarily disabled
1489 // because ZDC and VZERO are using it in order to initialize
1490 // their reconstructor objects. In the future one has to think
1491 // of propagating AliRunInfo to the reconstructors.
1492 // AliCDBManager::Instance()->UnloadFromCache("GRP/GRP/Data");
1493 }
1494
1495 if (!fGRPData) {
1496 AliError("No GRP entry found in OCDB!");
1497 return kFALSE;
1498 }
1499
1500 TString lhcState = fGRPData->GetLHCState();
1501 if (lhcState==AliGRPObject::GetInvalidString()) {
1502 AliError("GRP/GRP/Data entry: missing value for the LHC state ! Using UNKNOWN");
1503 lhcState = "UNKNOWN";
1504 }
1505
1506 TString beamType = fGRPData->GetBeamType();
1507 if (beamType==AliGRPObject::GetInvalidString()) {
1508 AliError("GRP/GRP/Data entry: missing value for the beam type ! Using UNKNOWN");
1509 beamType = "UNKNOWN";
1510 }
1511
1512 Float_t beamEnergy = fGRPData->GetBeamEnergy();
1513 if (beamEnergy==AliGRPObject::GetInvalidFloat()) {
1514 AliError("GRP/GRP/Data entry: missing value for the beam energy ! Using 0");
1515 beamEnergy = 0;
1516 }
1517
1518 TString runType = fGRPData->GetRunType();
1519 if (runType==AliGRPObject::GetInvalidString()) {
1520 AliError("GRP/GRP/Data entry: missing value for the run type ! Using UNKNOWN");
1521 runType = "UNKNOWN";
1522 }
1523
1524 Int_t activeDetectors = fGRPData->GetDetectorMask();
1525 if (activeDetectors==AliGRPObject::GetInvalidUInt()) {
1526 AliError("GRP/GRP/Data entry: missing value for the detector mask ! Using 1074790399");
1527 activeDetectors = 1074790399;
1528 }
1529 AliDebug(1, Form("activeDetectors = %d", activeDetectors));
1530
1531 fRunInfo = new AliRunInfo(lhcState, beamType, beamEnergy, runType, activeDetectors);
1532 fRunInfo->Dump();
1533
1534 // *** Dealing with the magnetic field map
1535
1536 if ( TGeoGlobalMagField::Instance()->IsLocked() ) {
1537 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(AliMagF::kOverrideGRP)) {
1538 AliInfo("ExpertMode!!! GRP information will be ignored !");
1539 AliInfo("ExpertMode!!! Running with the externally locked B field !");
1540 }
1541 else {
1542 AliInfo("Destroying existing B field instance!");
1543 delete TGeoGlobalMagField::Instance();
1544 }
1545 }
1546 if ( !TGeoGlobalMagField::Instance()->IsLocked() ) {
1547 // Construct the field map out of the information retrieved from GRP.
1548 Bool_t ok = kTRUE;
1549 // L3
1550 Float_t l3Current = fGRPData->GetL3Current((AliGRPObject::Stats)0);
1551 if (l3Current == AliGRPObject::GetInvalidFloat()) {
1552 AliError("GRP/GRP/Data entry: missing value for the L3 current !");
1553 ok = kFALSE;
1554 }
1555
1556 Char_t l3Polarity = fGRPData->GetL3Polarity();
1557 if (l3Polarity == AliGRPObject::GetInvalidChar()) {
1558 AliError("GRP/GRP/Data entry: missing value for the L3 polarity !");
1559 ok = kFALSE;
1560 }
1561
1562 // Dipole
1563 Float_t diCurrent = fGRPData->GetDipoleCurrent((AliGRPObject::Stats)0);
1564 if (diCurrent == AliGRPObject::GetInvalidFloat()) {
1565 AliError("GRP/GRP/Data entry: missing value for the dipole current !");
1566 ok = kFALSE;
1567 }
1568
1569 Char_t diPolarity = fGRPData->GetDipolePolarity();
1570 if (diPolarity == AliGRPObject::GetInvalidChar()) {
1571 AliError("GRP/GRP/Data entry: missing value for the dipole polarity !");
1572 ok = kFALSE;
1573 }
1574
1575 // read special bits for the polarity convention and map type
1576 Int_t polConvention = fGRPData->IsPolarityConventionLHC() ? AliMagF::kConvLHC : AliMagF::kConvDCS2008;
1577 Bool_t uniformB = fGRPData->IsUniformBMap();
1578
1579 if (ok) {
1580 AliMagF* fld = AliMagF::CreateFieldMap(TMath::Abs(l3Current) * (l3Polarity ? -1:1),
1581 TMath::Abs(diCurrent) * (diPolarity ? -1:1),
1582 polConvention,uniformB,beamEnergy, beamType.Data());
1583 if (fld) {
1584 TGeoGlobalMagField::Instance()->SetField( fld );
1585 TGeoGlobalMagField::Instance()->Lock();
1586 AliInfo("Running with the B field constructed out of GRP !");
1587 }
1588 else AliFatal("Failed to create a B field map !");
1589 }
1590 else AliFatal("B field is neither set nor constructed from GRP ! Exitig...");
1591 }
1592
1593 return kTRUE;
1594}
1595
1596//====================================================================================================================================================
1597
1598Bool_t AliMuonForwardTrackFinder::SetRunNumber() {
1599
1600 AliCDBManager *man = AliCDBManager::Instance();
1601
1602 if (!fRunLoader) {
1603 AliError("No run loader found!");
1604 return kFALSE;
1605 }
1606 else {
1607 fRunLoader->LoadHeader();
1608 // read run number from gAlice
1609 if (fRunLoader->GetHeader()) {
1610 man->SetRun(fRunLoader->GetHeader()->GetRun());
1611 fRunLoader->UnloadHeader();
1612 }
1613 else {
1614 AliError("No run-loader header found!");
1615 return kFALSE;
1616 }
1617 }
1618
1619 return kTRUE;
1620
1621}
1622
1623//====================================================================================================================================================
1624