Line data Source code
1 : /*************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 :
17 : #include "TTree.h"
18 : #include "TSystem.h"
19 : #include "TMath.h"
20 : #include "TArrayF.h"
21 : #include "TGrid.h"
22 :
23 : #include "AliCodeTimer.h"
24 : #include "AliLog.h"
25 : #include "AliGeomManager.h"
26 : #include "AliESDEvent.h"
27 : #include "AliESDMuonTrack.h"
28 : #include "AliESDMuonGlobalTrack.h"
29 : #include "AliMFTTracker.h"
30 : #include "AliMFTTrack.h"
31 : #include "AliRun.h"
32 : #include "AliRunLoader.h"
33 : #include "AliHeader.h"
34 : #include "AliGenEventHeader.h"
35 : #include "AliMFT.h"
36 : #include "AliMUONTrackExtrap.h"
37 : #include "AliMUONTrack.h"
38 : #include "AliMUONESDInterface.h"
39 : #include "AliMuonForwardTrack.h"
40 : #include "AliMUONConstants.h"
41 : #include "AliMUONRawClusterV2.h"
42 : #include "AliMFTTrack.h"
43 : #include "AliMFTTrackFinder.h"
44 : #include "AliMFTCATrack.h"
45 : #include "AliMFTCACell.h"
46 : #include "AliMFTGeometry.h"
47 : #include "AliMFTHalfSegmentation.h"
48 : #include "AliMFTHalfDiskSegmentation.h"
49 : #include "AliMFTTrackReconstructor.h"
50 : #include "AliMUONTrackParam.h"
51 :
52 : /// \cond CLASSIMP
53 12 : ClassImp(AliMFTTracker);
54 : /// \endcond
55 :
56 12 : const Double_t AliMFTTracker::fRadLengthSi = AliMFTConstants::fRadLengthSi;
57 :
58 : ///
59 : /// AliMFTTracker constructor
60 : ///
61 : AliMFTTracker::AliMFTTracker() :
62 0 : AliTracker(),
63 0 : fESD(0),
64 0 : fMFT(0),
65 0 : fSegmentation(NULL),
66 0 : fTrackFinder(0),
67 0 : fNPlanesMFT(0),
68 0 : fNPlanesMFTAnalyzed(0),
69 0 : fSigmaClusterCut(2),
70 0 : fScaleSigmaClusterCut(1.),
71 0 : fNMaxMissingMFTClusters(0),
72 0 : fGlobalTrackingDiverged(kFALSE),
73 0 : fCandidateTracks(0),
74 0 : fMFTTracks(0),
75 0 : fMUONTrack(0),
76 0 : fCurrentTrack(0),
77 0 : fFinalBestCandidate(0),
78 0 : fXExtrapVertex(0),
79 0 : fYExtrapVertex(0),
80 0 : fZExtrapVertex(0),
81 0 : fXExtrapVertexError(0),
82 0 : fYExtrapVertexError(0),
83 0 : fXVertexMC(0),
84 0 : fYVertexMC(0),
85 0 : fZVertexMC(0),
86 0 : fBransonCorrection(kFALSE)
87 0 : {
88 0 : AliMFTGeometry *mftGeo = AliMFTGeometry::Instance();
89 0 : AliMFTSegmentation * fSegmentation = mftGeo->GetSegmentation();
90 0 : if (!fSegmentation) AliFatal("No segmentation available");
91 :
92 0 : fMFT = (AliMFT*) gAlice->GetDetector("MFT");
93 0 : fTrackFinder = new AliMFTTrackFinder();
94 0 : fTrackFinder->Init(gSystem->ExpandPathName("$(ALICE_ROOT)/ITSMFT/MFT/data/param_10ch.txt" ));
95 0 : SetNPlanesMFT(AliMFTConstants::kNDisks);
96 0 : AliMUONTrackExtrap::SetField(); // set the magnetic field for track extrapolations
97 :
98 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
99 0 : fMFTClusterArray[iPlane] = new TClonesArray("AliMFTCluster");
100 0 : fMFTClusterArrayFront[iPlane] = new TClonesArray("AliMFTCluster");
101 0 : fMFTClusterArrayBack[iPlane] = new TClonesArray("AliMFTCluster");
102 0 : fMFTClusterArray[iPlane] -> SetOwner(kTRUE);
103 0 : fMFTClusterArrayFront[iPlane] -> SetOwner(kTRUE);
104 0 : fMFTClusterArrayBack[iPlane] -> SetOwner(kTRUE);
105 0 : fMinResearchRadiusAtPlane[iPlane] = 0.;
106 : }
107 :
108 0 : fCandidateTracks = new TClonesArray("AliMuonForwardTrack",50000);
109 0 : fMFTTracks = new TClonesArray("AliMFTTrack",100);
110 0 : fMFTTracks -> SetOwner(kTRUE);
111 :
112 0 : }
113 :
114 : //====================================================================================================================================================
115 :
116 0 : AliMFTTracker::~AliMFTTracker() {
117 :
118 : // destructor
119 :
120 0 : delete fTrackFinder;
121 :
122 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
123 0 : fMFTClusterArray[iPlane] -> Delete();
124 0 : delete fMFTClusterArray[iPlane];
125 0 : delete fMFTClusterArrayFront[iPlane];
126 0 : delete fMFTClusterArrayBack[iPlane];
127 : }
128 :
129 0 : delete fCandidateTracks;
130 0 : delete fMFTTracks;
131 :
132 0 : }
133 :
134 : //==============================================================================================
135 : ///
136 : /// Loads the MFT clusters
137 : ///
138 : /// \param cTree TTree containing the ALiMFTCluster objects
139 : ///
140 : /// \return 0
141 : Int_t AliMFTTracker::LoadClusters(TTree *cTree) {
142 :
143 0 : AliCodeTimerAuto("",0);
144 :
145 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
146 0 : AliDebug(1, Form("Setting Address for Branch Plane_%02d", iPlane));
147 0 : cTree->SetBranchAddress(Form("Plane_%02d",iPlane), &fMFTClusterArray[iPlane]);
148 : }
149 :
150 0 : if (!cTree->GetEvent()) return kFALSE;
151 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
152 0 : AliInfo(Form("plane %02d: nClusters = %d", iPlane, fMFTClusterArray[iPlane]->GetEntries()));
153 : }
154 :
155 : //AddClustersFromUnderlyingEvent();
156 : //AddClustersFromPileUpEvents();
157 :
158 0 : SeparateFrontBackClusters();
159 :
160 0 : fTrackFinder->LoadClusters(fMFTClusterArrayFront, fMFTClusterArrayBack);
161 :
162 0 : return 0;
163 :
164 0 : }
165 : //==============================================================================================
166 : ///
167 : /// Loads the tracks found by the track finder
168 : ///
169 :
170 : void AliMFTTracker::LoadTracks() {
171 :
172 0 : AliCodeTimerAuto("",0);
173 :
174 0 : Int_t nTracks = fTrackFinder->GetNtracks();
175 : AliMFTCATrack * catrack = NULL;
176 0 : for (Int_t i = 0 ; i < nTracks; i++) {
177 0 : catrack = fTrackFinder->GetTrack(i);
178 0 : new ((*fMFTTracks)[i]) AliMFTTrack(catrack);
179 0 : LinearFit((AliMFTTrack*)fMFTTracks->At(i));
180 : }
181 :
182 0 : }
183 :
184 : //==============================================================================================
185 :
186 : void AliMFTTracker::UnloadClusters() {
187 :
188 : //--------------------------------------------------------------------
189 : // This function unloads MFT clusters
190 : //--------------------------------------------------------------------
191 :
192 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
193 0 : fMFTClusterArray[iPlane] -> Clear("C");
194 0 : fMFTClusterArrayFront[iPlane] -> Clear("C");
195 0 : fMFTClusterArrayBack[iPlane] -> Clear("C");
196 : }
197 :
198 0 : }
199 :
200 : //==============================================================================================
201 :
202 : Int_t AliMFTTracker::Clusters2Tracks(AliESDEvent *event) {
203 0 : AliCodeTimerAuto("",0);
204 :
205 : // Int_t nGlobalTrack = 0;
206 : //--------------------------------------------------------------------
207 : // This functions reconstructs the Muon Forward Tracks
208 : // The clusters must be already loaded !
209 : //--------------------------------------------------------------------
210 :
211 : // Tree of AliMuonForwardTrack objects. Created outside the ESD framework for cross-check purposes
212 :
213 : Double_t xVertex = 0., yVertex = 0., zVertex = 0.;
214 0 : Double_t zvr[2] = {0.,0.};
215 :
216 0 : const AliESDVertex* esdVert = event->GetVertex();
217 0 : if (esdVert->GetNContributors() > 0 || !strcmp(esdVert->GetTitle(),"vertexer: smearMC")) {
218 0 : xVertex = esdVert->GetX();
219 0 : yVertex = esdVert->GetY();
220 0 : zVertex = esdVert->GetZ();
221 0 : zvr[0] = zVertex - 3.*AliMFTConstants::fPrimaryVertexResZ; // 5 microns
222 0 : zvr[1] = zVertex + 3.*AliMFTConstants::fPrimaryVertexResZ;
223 0 : GetVertexFromMC();
224 : } else {
225 0 : printf("No vertex in ESD! Get it from MC.\n");
226 0 : GetVertexFromMC();
227 0 : xVertex = fXExtrapVertex;
228 0 : yVertex = fYExtrapVertex;
229 0 : zVertex = fZExtrapVertex;
230 0 : zvr[0] = zVertex - 3.*AliMFTConstants::fPrimaryVertexResZ; // 5 microns
231 0 : zvr[1] = zVertex + 3.*AliMFTConstants::fPrimaryVertexResZ;
232 : }
233 :
234 : //printf("Vertex: %f %f %f \n",zvr[0],zvr[1],zVertex);
235 0 : fTrackFinder->SetZVertRange(zvr,zVertex);
236 :
237 0 : fTrackFinder->FindTracks();
238 0 : fTrackFinder->BuildRoads();
239 0 : fTrackFinder->FilterTracks();
240 :
241 0 : LoadTracks();
242 :
243 : //AliInfo(" ------ Print TrackFinder Out -------");
244 : //
245 : //fTrackFinder->PrintAll();
246 :
247 0 : AliInfo("Track Finder Done");
248 : // Reconstruction of the MFT tracks from the Tracks found by the Cellular Automaton
249 :
250 0 : AliMFTTrackReconstructor * trackReco = new AliMFTTrackReconstructor();
251 0 : trackReco->EventReconstruct(fMFTTracks);
252 : // ----> Standalone track reconstruction done
253 :
254 0 : Int_t nTracksMUON = event->GetNumberOfMuonTracks();
255 0 : Int_t nTracksMFT = fTrackFinder->GetNtracks();
256 :
257 : AliMFTCATrack * caTrack = NULL;
258 : AliMUONTrack * muonTrack = NULL;
259 : AliMFTCACell * caCell = NULL;
260 : AliMuonForwardTrack * mfwdTrack = NULL;
261 :
262 : Double_t equivalentSilicon = 0.0028;
263 : Double_t equivalentSiliconBeforeFront = 0.0028;
264 : Double_t equivalentSiliconBeforeBack = 0.0050;
265 :
266 0 : Double_t zEndOfMFTTrack,
267 : xTr[AliMFTConstants::fNMaxPlanes],
268 : yTr[AliMFTConstants::fNMaxPlanes],
269 : zTr[AliMFTConstants::fNMaxPlanes];
270 0 : Int_t planeID[AliMFTConstants::fNMaxPlanes];
271 :
272 : Double_t phic, phicp, phis;
273 : Short_t trackCharge;
274 : Double_t caTrackPhi, caTrackThe, caTrackOrg[3], caTrackBegOfAbs[3];
275 : Double_t Ux, Uy, Uz, Tx, Ty, Tz;
276 : Double_t addChi2TrackAtCluster = 0.;
277 0 : AliMUONTrackParam trackParamMM, trackParamMM0;
278 : AliMUONRawCluster *muonCluster = 0x0;
279 :
280 0 : AliMFTCluster *mftCluster[AliMFTConstants::fNMaxPlanes];
281 0 : for (Int_t i = 0; i < AliMFTConstants::fNMaxPlanes; i++) {
282 0 : mftCluster[i] = 0x0;
283 : }
284 :
285 : Bool_t saveAllMatch = kTRUE;
286 :
287 : Double_t chi2cut = 2.0;
288 :
289 0 : AliInfo(Form("Number of ESD MUON tracks: %d\n", nTracksMUON));
290 :
291 0 : TFile *outputFileMFTTracks = new TFile("MFT.Tracks.root", "update");
292 :
293 : Int_t myEventID = 0;
294 0 : while (outputFileMFTTracks->cd(Form("Event%d",myEventID))) myEventID++;
295 0 : outputFileMFTTracks->mkdir(Form("Event%d",myEventID));
296 0 : outputFileMFTTracks->cd(Form("Event%d",myEventID));
297 :
298 0 : TTree *outputTreeMFTTracks = new TTree("MFTTracks", "Tree of MFT tracks");
299 0 : TClonesArray *mftTracks = new TClonesArray("AliMFTCATrack");
300 0 : outputTreeMFTTracks->Branch("tracks", &mftTracks);
301 :
302 0 : TTree *outputTreeMuonForwardTracks = new TTree("MuonForwardTracks", "Tree of muon forward racks");
303 0 : TClonesArray *mfwdTracks = new TClonesArray("AliMuonForwardTrack");
304 0 : outputTreeMuonForwardTracks->Branch("tracks", &mfwdTracks);
305 :
306 0 : TTree *outputTreeEvent = new TTree("Events", "Tree of events");
307 0 : outputTreeEvent->Branch("fXVertexMC", &fXVertexMC);
308 0 : outputTreeEvent->Branch("fYVertexMC", &fYVertexMC);
309 0 : outputTreeEvent->Branch("fZVertexMC", &fZVertexMC);
310 :
311 0 : TTree *outputTreeMFTCells = new TTree("MFTCells", "Tree of MFT CA cells");
312 0 : TClonesArray *mftCells = new TClonesArray("AliMFTCACell");
313 0 : outputTreeMFTCells->Branch("cells", &mftCells);
314 :
315 : Int_t iTrack=0, iTrackMatchA=0, iTrackMatchB=0, iTotalCells=0;
316 0 : while (iTrack < nTracksMUON) {
317 :
318 0 : const AliESDMuonTrack *esdTrack = event->GetMuonTrack(iTrack);
319 :
320 0 : trackCharge = esdTrack->Charge();
321 :
322 0 : if (muonTrack) delete muonTrack;
323 0 : muonTrack = new AliMUONTrack();
324 0 : AliMUONESDInterface::ESDToMUON(*esdTrack, *muonTrack, kFALSE);
325 :
326 0 : if (!muonTrack->GetTrackParamAtCluster()->First()) {
327 0 : AliInfo("Skipping track, no parameters available!!!");
328 0 : iTrack++;
329 0 : continue;
330 : }
331 : //printf("Muon track %d start chi2: %f , %f , %d \n",iTrack,muonTrack->GetGlobalChi2(),muonTrack->GetGlobalChi2()/muonTrack->GetNDF(),muonTrack->GetNDF());
332 :
333 0 : if (mfwdTrack) delete mfwdTrack;
334 0 : mfwdTrack = new AliMuonForwardTrack(*muonTrack);
335 :
336 : // go with the track param to the vertex x,y,z (with Branson) or
337 : // to vertex z (without Branson)
338 :
339 0 : trackParamMM = (*((AliMUONTrackParam*)(mfwdTrack->GetTrackParamAtCluster()->First())));
340 :
341 0 : if (fBransonCorrection) {
342 0 : AliMUONTrackExtrap::ExtrapToVertex(&trackParamMM,fXExtrapVertex,fYExtrapVertex,fZExtrapVertex,fXExtrapVertexError,fYExtrapVertexError);
343 : } else {
344 0 : AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(&trackParamMM,fZExtrapVertex);
345 : }
346 : // copy to save this track param
347 0 : trackParamMM0 = trackParamMM;
348 :
349 0 : for (Int_t iTrackMFT = 0 ; iTrackMFT < nTracksMFT; iTrackMFT++) {
350 :
351 0 : caTrack = fTrackFinder->GetTrack(iTrackMFT);
352 :
353 0 : caTrackPhi = caTrack->GetPhi();
354 0 : caTrackThe = caTrack->GetTheta();
355 0 : caTrackOrg[0] = caTrack->GetVertX();
356 0 : caTrackOrg[1] = caTrack->GetVertY();
357 0 : caTrackOrg[2] = caTrack->GetVertZ();
358 :
359 0 : caTrackPhi *= TMath::DegToRad();
360 0 : caTrackThe *= TMath::DegToRad();
361 :
362 0 : Ux = TMath::Sin(caTrackThe)*TMath::Cos(caTrackPhi);
363 0 : Uy = TMath::Sin(caTrackThe)*TMath::Sin(caTrackPhi);
364 0 : Uz = TMath::Cos(caTrackThe);
365 : /*
366 : caTrackBegOfAbs[2] = zBegOfAbsorber;
367 :
368 : Tz = caTrackBegOfAbs[2] - caTrackOrg[2];
369 : caTrackBegOfAbs[0] = caTrackOrg[0] + Ux*Tz;
370 : caTrackBegOfAbs[1] = caTrackOrg[1] + Uy*Tz;
371 :
372 : printf("CATrack: %3d %10.4f %10.4f %10.4f %d \n",iTrackMFT,caTrackBegOfAbs[0],caTrackBegOfAbs[1],caTrackBegOfAbs[2],caTrack->GetMCindex());
373 : */
374 : // extract hit x,y,z from the cells
375 : Int_t mftClsId1, mftClsId2, layer1, layer2;
376 : Int_t nptr = 0;
377 : AliMFTCluster *cls1, *cls2;
378 0 : for (Int_t iCell = 0; iCell < caTrack->GetNcells(); iCell++) {
379 :
380 0 : caCell = caTrack->GetCell(iCell);
381 0 : caTrack->SetCellGID(iCell,iTotalCells);
382 :
383 0 : mftClsId1 = caCell->GetMFTClsId()[0];
384 0 : mftClsId2 = caCell->GetMFTClsId()[1];
385 0 : layer1 = caCell->GetLayers()[0];
386 0 : layer2 = caCell->GetLayers()[1];
387 0 : if (layer1%2 == 0) { // FRONT
388 0 : cls1 = (AliMFTCluster*)fMFTClusterArrayFront[layer1/2]->At(mftClsId1);
389 0 : } else { // BACK
390 0 : cls1 = (AliMFTCluster*)fMFTClusterArrayBack[layer1/2]->At(mftClsId1);
391 : }
392 0 : if (layer2%2 == 0) { // FRONT
393 0 : cls2 = (AliMFTCluster*)fMFTClusterArrayFront[layer2/2]->At(mftClsId2);
394 0 : } else { // BACK
395 0 : cls2 = (AliMFTCluster*)fMFTClusterArrayBack[layer2/2]->At(mftClsId2);
396 : }
397 :
398 : //printf("Cell %5d MFTClsId %5d %5d \n",iCell,mftClsId1,mftClsId2);
399 : //printf("Cls1: %10.4f %10.4f %10.4f \n",cls1->GetX(),cls1->GetY(),cls1->GetZ());
400 : //printf("Cls2: %10.4f %10.4f %10.4f \n",cls2->GetX(),cls2->GetY(),cls2->GetZ());
401 :
402 0 : new ((*mftCells)[iTotalCells++]) AliMFTCACell(*caCell);
403 :
404 0 : if (nptr == 0) {
405 0 : xTr[nptr] = caCell->GetHit2()[0];
406 0 : yTr[nptr] = caCell->GetHit2()[1];
407 0 : zTr[nptr] = caCell->GetHit2()[2];
408 0 : planeID[nptr] = caCell->GetLayers()[1];
409 0 : mftCluster[nptr] = cls2;
410 0 : nptr++;
411 0 : xTr[nptr] = caCell->GetHit1()[0];
412 0 : yTr[nptr] = caCell->GetHit1()[1];
413 0 : zTr[nptr] = caCell->GetHit1()[2];
414 0 : planeID[nptr] = caCell->GetLayers()[0];
415 0 : mftCluster[nptr] = cls1;
416 0 : nptr++;
417 0 : } else {
418 0 : xTr[nptr] = caCell->GetHit1()[0];
419 0 : yTr[nptr] = caCell->GetHit1()[1];
420 0 : zTr[nptr] = caCell->GetHit1()[2];
421 0 : planeID[nptr] = caCell->GetLayers()[0];
422 0 : mftCluster[nptr] = cls1;
423 0 : nptr++;
424 : }
425 : } // end loop over cells
426 :
427 : // estimate the charge sign
428 : phis = 0.;
429 0 : for (Int_t iptr = 0; iptr < nptr; iptr++) {
430 0 : phic = TMath::ATan(yTr[iptr]/xTr[iptr])*TMath::RadToDeg();
431 0 : phic += 90.;
432 0 : if (iptr > 0) {
433 0 : phis += phic-phicp;
434 0 : }
435 : phicp = phic;
436 : }
437 :
438 0 : caCell = caTrack->GetCell(0);
439 :
440 0 : caTrack->SetChargeSign(-(Short_t)(TMath::Sign(1.,phis)));
441 :
442 : // match by MC label
443 0 : if (saveAllMatch || caTrack->GetMCindex() == muonTrack->GetMCLabel()) {
444 :
445 0 : if (saveAllMatch) {
446 0 : if (muonTrack) delete muonTrack;
447 0 : muonTrack = new AliMUONTrack();
448 0 : AliMUONESDInterface::ESDToMUON(*esdTrack, *muonTrack, kFALSE);
449 0 : if (mfwdTrack) delete mfwdTrack;
450 0 : mfwdTrack = new AliMuonForwardTrack(*muonTrack);
451 0 : trackParamMM = trackParamMM0;
452 : }
453 :
454 0 : for (Int_t iptr = 0; iptr < nptr; iptr++) {
455 :
456 0 : muonCluster = mftCluster[iptr]->CreateMUONCluster();
457 : //printf("MFTCluster: %3d %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f \n",iptr,mftCluster[iptr]->GetX(),mftCluster[iptr]->GetY(),mftCluster[iptr]->GetZ(),xTr[iptr],yTr[iptr],zTr[iptr],muonCluster->GetX(),muonCluster->GetY(),muonCluster->GetZ());
458 :
459 : // extrapolation to z
460 : //printf("Extrap to (b): %10.4f \n",muonCluster->GetZ());
461 0 : AliMUONTrackExtrap::ExtrapToZCov(&trackParamMM,muonCluster->GetZ());
462 :
463 : // add MCS (Multiple Coulomb Scattering)
464 : // front/back correct ?
465 0 : if (iptr > 0) {
466 0 : if (planeID[iptr]%2 == 0) {
467 : // back
468 0 : AliMUONTrackExtrap::AddMCSEffect(&trackParamMM,
469 0 : (equivalentSilicon+equivalentSiliconBeforeBack)/fRadLengthSi,-1.);
470 : } else {
471 : // front
472 : // ... this is zero, anyway ...
473 0 : AliMUONTrackExtrap::AddMCSEffect(&trackParamMM,
474 0 : (equivalentSilicon+equivalentSiliconBeforeFront)/fRadLengthSi,-1.);
475 : }
476 : }
477 :
478 : //printf("1b (%2d): %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %f\n",iptr,trackParamMM.GetNonBendingCoor(), trackParamMM.GetBendingCoor(),trackParamMM.GetZ(),muonCluster->GetX(),muonCluster->GetY(),muonCluster->GetZ(),trackParamMM.GetTrackChi2());
479 0 : addChi2TrackAtCluster = RunKalmanFilter(trackParamMM,*muonCluster);
480 0 : trackParamMM.SetTrackChi2(trackParamMM.GetTrackChi2()+addChi2TrackAtCluster);
481 :
482 0 : mfwdTrack->AddTrackParamAtMFTCluster(trackParamMM,*muonCluster,iptr);
483 0 : mfwdTrack->SetGlobalChi2(trackParamMM.GetTrackChi2());
484 0 : mfwdTrack->SetTrackMCId(caTrack->GetMCindex());
485 : //printf("2b: %10.4f %10.4f %10.4f %f\n",trackParamMM.GetNonBendingCoor(), trackParamMM.GetBendingCoor(),trackParamMM.GetZ(),trackParamMM.GetTrackChi2());
486 0 : delete muonCluster;
487 :
488 : } // end MFT cluster loop
489 :
490 0 : if (saveAllMatch) {
491 0 : if (mfwdTrack->GetNormalizedChi2() < chi2cut) {
492 : //printf("Muon forward track %2d:%2d end chi2: %f , %f , %d \n",iTrack,iTrackMFT,mfwdTrack->GetGlobalChi2(),mfwdTrack->GetGlobalChi2()/mfwdTrack->GetNDF(),mfwdTrack->GetNDF());
493 0 : new ((*mfwdTracks)[iTrackMatchB++]) AliMuonForwardTrack(*mfwdTrack);
494 : } else {
495 : //printf("... not matched (b) ... %2d:%2d end chi2: %f , %f , %d \n",iTrack,iTrackMFT,mfwdTrack->GetGlobalChi2(),mfwdTrack->GetGlobalChi2()/mfwdTrack->GetNDF(),mfwdTrack->GetNDF());
496 : }
497 : }
498 :
499 : } // end save all || match by MC label
500 : } // end MFT track loop
501 :
502 0 : if (!saveAllMatch) {
503 :
504 : //printf("Muon forward track %d end chi2: %f , %f , %d \n",iTrack,mfwdTrack->GetGlobalChi2(),mfwdTrack->GetGlobalChi2()/mfwdTrack->GetNDF(),mfwdTrack->GetNDF());
505 0 : if (mfwdTrack->GetNormalizedChi2() < chi2cut) {
506 0 : new ((*mfwdTracks)[iTrackMatchB++]) AliMuonForwardTrack(*mfwdTrack);
507 : }
508 :
509 : }
510 : /*
511 : for (Int_t i = 0; i < muonTrack->GetTrackParamAtCluster()->GetEntries(); i++) {
512 : AliMUONTrackParam trackParamMM0(*((AliMUONTrackParam*)(muonTrack->GetTrackParamAtCluster()->At(i))));
513 : printf("%d %10.4f %10.4f \n",i,trackParamMM0.P(),trackParamMM0.GetZ());
514 : }
515 : */
516 0 : iTrack++;
517 :
518 0 : } // end MUON track loop
519 :
520 0 : for (Int_t iTrackMFT = 0 ; iTrackMFT < nTracksMFT; iTrackMFT++) {
521 :
522 0 : caTrack = fTrackFinder->GetTrack(iTrackMFT);
523 0 : new ((*mftTracks)[iTrackMFT]) AliMFTCATrack(*caTrack);
524 :
525 : }
526 :
527 0 : outputTreeMFTTracks->Fill();
528 0 : outputTreeMFTTracks->Write();
529 0 : outputTreeMuonForwardTracks->Fill();
530 0 : outputTreeMuonForwardTracks->Write();
531 0 : outputTreeMFTCells->Fill();
532 0 : outputTreeMFTCells->Write();
533 0 : outputTreeEvent->Fill();
534 0 : outputTreeEvent->Write();
535 :
536 0 : outputFileMFTTracks->Close();
537 :
538 0 : mftTracks->Delete();
539 0 : delete mftTracks;
540 :
541 0 : mftCells->Delete();
542 0 : delete mftCells;
543 :
544 0 : mfwdTracks->Delete();
545 0 : delete mfwdTracks;
546 :
547 0 : fTrackFinder->Clear("");
548 :
549 : return 0;
550 :
551 0 : }
552 :
553 : //=========================================================================================================================================
554 :
555 : Double_t AliMFTTracker::RunKalmanFilter(AliMUONTrackParam &trackParamAtCluster, AliMUONVCluster &cluster)
556 : {
557 : /// Compute new track parameters and their covariances including new cluster using kalman filter
558 : /// return the additional track chi2
559 : /// copied from AliMUONTrackReconstructorK::RunKalmanFilter
560 0 : AliDebug(1,"Enter RunKalmanFilter");
561 :
562 : // Get actual track parameters (p)
563 0 : TMatrixD param(trackParamAtCluster.GetParameters());
564 :
565 : // Get new cluster parameters (m)
566 0 : TMatrixD clusterParam(5,1);
567 0 : clusterParam.Zero();
568 0 : clusterParam(0,0) = cluster.GetX();
569 0 : clusterParam(2,0) = cluster.GetY();
570 :
571 : // Compute the actual parameter weight (W)
572 0 : TMatrixD paramWeight(trackParamAtCluster.GetCovariances());
573 0 : if (paramWeight.Determinant() != 0) {
574 0 : paramWeight.Invert();
575 : } else {
576 0 : AliWarning(" Determinant = 0");
577 0 : return 2.*AliMUONTrack::MaxChi2();
578 : }
579 :
580 : // Compute the new cluster weight (U)
581 0 : TMatrixD clusterWeight(5,5);
582 0 : clusterWeight.Zero();
583 0 : clusterWeight(0,0) = 1. / cluster.GetErrX2();
584 0 : clusterWeight(2,2) = 1. / cluster.GetErrY2();
585 :
586 : // Compute the new parameters covariance matrix ( (W+U)^-1 )
587 0 : TMatrixD newParamCov(paramWeight,TMatrixD::kPlus,clusterWeight);
588 0 : if (newParamCov.Determinant() != 0) {
589 0 : newParamCov.Invert();
590 : } else {
591 0 : AliWarning(" Determinant = 0");
592 0 : return 2.*AliMUONTrack::MaxChi2();
593 : }
594 :
595 : // Save the new parameters covariance matrix
596 0 : trackParamAtCluster.SetCovariances(newParamCov);
597 :
598 : // Compute the new parameters (p' = ((W+U)^-1)U(m-p) + p)
599 0 : TMatrixD tmp(clusterParam,TMatrixD::kMinus,param);
600 0 : TMatrixD tmp2(clusterWeight,TMatrixD::kMult,tmp); // U(m-p)
601 0 : TMatrixD newParam(newParamCov,TMatrixD::kMult,tmp2); // ((W+U)^-1)U(m-p)
602 0 : newParam += param; // ((W+U)^-1)U(m-p) + p
603 :
604 : // Save the new parameters
605 0 : trackParamAtCluster.SetParameters(newParam);
606 :
607 : // Compute the additional chi2 (= ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m))
608 0 : tmp = newParam; // p'
609 0 : tmp -= param; // (p'-p)
610 0 : TMatrixD tmp3(paramWeight,TMatrixD::kMult,tmp); // W(p'-p)
611 0 : TMatrixD addChi2Track(tmp,TMatrixD::kTransposeMult,tmp3); // ((p'-p)^-1)W(p'-p)
612 0 : tmp = newParam; // p'
613 0 : tmp -= clusterParam; // (p'-m)
614 0 : TMatrixD tmp4(clusterWeight,TMatrixD::kMult,tmp); // U(p'-m)
615 0 : addChi2Track += TMatrixD(tmp,TMatrixD::kTransposeMult,tmp4); // ((p'-p)^-1)W(p'-p) + ((p'-m)^-1)U(p'-m)
616 :
617 0 : return addChi2Track(0,0);
618 :
619 0 : }
620 :
621 : //=========================================================================================================================================
622 :
623 : void AliMFTTracker::SeparateFrontBackClusters() {
624 0 : AliCodeTimerAuto("",0);
625 :
626 0 : AliMFTGeometry *mftGeo = AliMFTGeometry::Instance();
627 0 : AliMFTSegmentation * seg = mftGeo->GetSegmentation();
628 :
629 0 : AliMFTHalfSegmentation *halfSeg[2];
630 0 : for (int i=0; i<2; i++) halfSeg[i] = seg->GetHalf(i);
631 :
632 0 : AliMFTHalfDiskSegmentation *halfDiskSeg[2];
633 :
634 0 : for (Int_t iPlane=0; iPlane<AliMFTConstants::kNDisks; iPlane++) {
635 0 : fMFTClusterArrayFront[iPlane]->Delete();
636 0 : fMFTClusterArrayBack[iPlane] ->Delete();
637 0 : for (int i=0; i<2; i++) halfDiskSeg[i] = halfSeg[i]->GetHalfDisk(iPlane);
638 :
639 0 : for (Int_t iCluster=0; iCluster<fMFTClusterArray[iPlane]->GetEntries(); iCluster++) {
640 0 : AliMFTCluster *cluster = (AliMFTCluster*) fMFTClusterArray[iPlane]->At(iCluster);
641 0 : if (mftGeo->GetLadderID(cluster->GetDetElemID())<(halfDiskSeg[mftGeo->GetHalfMFTID(cluster->GetDetElemID())]->GetNLadders())/2) {
642 0 : new ((*fMFTClusterArrayFront[iPlane])[fMFTClusterArrayFront[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
643 : }
644 : else {
645 0 : new ((*fMFTClusterArrayBack[iPlane])[fMFTClusterArrayBack[iPlane]->GetEntries()]) AliMFTCluster(*cluster);
646 : }
647 : }
648 : }
649 :
650 0 : }
651 :
652 : //======================================================================================================================================
653 :
654 : void AliMFTTracker::GetVertexFromMC() {
655 :
656 0 : AliRunLoader *runLoader = AliRunLoader::Open("galice.root");
657 0 : if (!runLoader) {
658 0 : AliError("no run loader found in file galice.root");
659 0 : return;
660 : }
661 :
662 0 : runLoader->CdGAFile();
663 0 : runLoader->LoadgAlice();
664 0 : runLoader->LoadHeader();
665 0 : runLoader->GetEvent(gAlice->GetEvNumber());
666 :
667 0 : TArrayF vtx(3);
668 0 : runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vtx);
669 0 : AliInfo(Form("Primary vertex from MC found in (%f, %f, %f)\n",vtx[0], vtx[1], vtx[2]));
670 :
671 0 : fXVertexMC = vtx[0];
672 0 : fYVertexMC = vtx[1];
673 0 : fZVertexMC = vtx[2];
674 :
675 0 : fXExtrapVertex = gRandom->Gaus(vtx[0], AliMFTConstants::fPrimaryVertexResX);
676 0 : fYExtrapVertex = gRandom->Gaus(vtx[1], AliMFTConstants::fPrimaryVertexResY);
677 0 : fZExtrapVertex = gRandom->Gaus(vtx[2], AliMFTConstants::fPrimaryVertexResZ);
678 0 : fXExtrapVertexError = AliMFTConstants::fXVertexTolerance;
679 0 : fYExtrapVertexError = AliMFTConstants::fYVertexTolerance;
680 0 : AliInfo(Form("Set ESD vertex from MC (%f +/- %f, %f +/- %f, %f)",
681 : fXExtrapVertex,fXExtrapVertexError,fYExtrapVertex,fYExtrapVertexError,fZExtrapVertex));
682 :
683 0 : }
684 :
685 : //======================================================================================================================================
686 :
687 : void AliMFTTracker::AddClustersFromUnderlyingEvent() {
688 :
689 0 : AliInfo("Adding clusters from underlying event");
690 :
691 0 : if (!fMFT) return;
692 :
693 0 : TGrid::Connect("alien://");
694 :
695 0 : TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForUnderlyingEvent());
696 0 : if (!fileWithClustersToAdd) return;
697 0 : if (!(fileWithClustersToAdd->IsOpen())) return;
698 0 : if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetUnderlyingEventID())))) return;
699 :
700 0 : TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
701 : TTree *treeIn = 0;
702 :
703 0 : for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
704 :
705 0 : treeIn = (TTree*) gDirectory->Get("TreeR");
706 :
707 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
708 0 : if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
709 0 : for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
710 0 : return;
711 : }
712 0 : else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
713 : }
714 :
715 0 : treeIn -> GetEntry(0);
716 :
717 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
718 0 : printf("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries());
719 0 : Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
720 0 : for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
721 0 : AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
722 0 : for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
723 0 : new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
724 : }
725 0 : printf("after = %d\n",fMFTClusterArray[iPlane]->GetEntries());
726 : }
727 :
728 0 : for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
729 :
730 0 : }
731 :
732 : //======================================================================================================================================
733 :
734 : void AliMFTTracker::AddClustersFromPileUpEvents() {
735 :
736 0 : AliInfo("Adding clusters from pile-up event(s)");
737 :
738 0 : if (!fMFT) return;
739 :
740 0 : TGrid::Connect("alien://");
741 :
742 0 : TFile* fileWithClustersToAdd = TFile::Open(fMFT->GetFileNameForPileUpEvents());
743 0 : if (!fileWithClustersToAdd) return;
744 0 : if (!(fileWithClustersToAdd->IsOpen())) return;
745 :
746 0 : TClonesArray *recPointsPerPlaneToAdd[AliMFTConstants::fNMaxPlanes] = {0};
747 : TTree *treeIn = 0;
748 :
749 0 : for (Int_t iPileUp=0; iPileUp<AliMFTConstants::fNMaxPileUpEvents; iPileUp++) {
750 :
751 0 : if (!(fileWithClustersToAdd->cd(Form("Event%d",fMFT->GetPileUpEventID(iPileUp))))) continue;
752 :
753 0 : for (Int_t iPlane=0; iPlane<AliMFTConstants::fNMaxPlanes; iPlane++) recPointsPerPlaneToAdd[iPlane] = new TClonesArray("AliMFTCluster");
754 :
755 0 : treeIn = (TTree*) gDirectory->Get("TreeR");
756 :
757 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
758 0 : if (!(treeIn->GetBranch(Form("Plane_%02d",iPlane)))) {
759 0 : for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
760 0 : return;
761 : }
762 0 : else treeIn->SetBranchAddress(Form("Plane_%02d",iPlane), &(recPointsPerPlaneToAdd[iPlane]));
763 : }
764 :
765 0 : treeIn -> GetEntry(0);
766 :
767 0 : for (Int_t iPlane=0; iPlane<fNPlanesMFT; iPlane++) {
768 0 : AliInfo(Form("plane %d -> before = %d ",iPlane,fMFTClusterArray[iPlane]->GetEntries()));
769 0 : Int_t nClusters = recPointsPerPlaneToAdd[iPlane]->GetEntries();
770 0 : for (Int_t iCluster=0; iCluster<nClusters; iCluster++) {
771 0 : AliMFTCluster *newCluster = (AliMFTCluster*) recPointsPerPlaneToAdd[iPlane]->At(iCluster);
772 0 : for (Int_t iTrack=0; iTrack<newCluster->GetNMCTracks(); iTrack++) newCluster->SetMCLabel(iTrack, newCluster->GetMCLabel(iTrack)+AliMFTConstants::fLabelOffsetMC);
773 0 : new ((*fMFTClusterArray[iPlane])[fMFTClusterArray[iPlane]->GetEntries()]) AliMFTCluster(*newCluster);
774 : }
775 0 : AliInfo(Form("after = %d\n",fMFTClusterArray[iPlane]->GetEntries()));
776 : }
777 :
778 0 : for (Int_t jPlane=0; jPlane<AliMFTConstants::fNMaxPlanes; jPlane++) delete recPointsPerPlaneToAdd[jPlane];
779 :
780 0 : }
781 :
782 0 : }
783 :
784 : //======================================================================================================================================
785 : //___________________________________________________________________________
786 : Bool_t AliMFTTracker::LinearFit(AliMFTTrack * track) {
787 :
788 0 : if (!track) return 0;
789 0 : if (!track->GetCATrack()) return 0;
790 :
791 0 : AliMFTCATrack * caTrack = track->GetCATrack();
792 0 : Int_t nCells = caTrack->GetNcells();
793 0 : AliDebug(1,Form("NCell = %d ",nCells));
794 :
795 0 : Double_t xTrErrDet = 0.0025/TMath::Sqrt(12.);
796 0 : Double_t yTrErrDet = 0.0025/TMath::Sqrt(12.);
797 : Double_t xTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
798 : Double_t yTrErrMS = 0.00055; // estimated at p = 5.5 GeV/c
799 :
800 : Int_t nDet=0;
801 : const Int_t nMaxCell = 10;
802 :
803 0 : if (nCells>nMaxCell) AliError(Form("Number of Cell = %d; Bigger than allowed value = %d", nCells,nMaxCell));
804 :
805 0 : Double_t xcl[nMaxCell];
806 0 : Double_t ycl[nMaxCell];
807 0 : Double_t zcl[nMaxCell];
808 0 : Double_t xerr[nMaxCell];
809 0 : Double_t yerr[nMaxCell];
810 : // Double_t &a; Double_t &ae; Double_t &b; Double_t &be;
811 : // Int_t skip;
812 :
813 : AliMFTCACell *caCell = NULL;
814 0 : for (int iCell=0; iCell<nCells; iCell++) {
815 0 : caCell = caTrack->GetCell(iCell);
816 0 : if(iCell==0){
817 0 : xcl[nDet] = caCell->GetHit1()[0];
818 0 : ycl[nDet] = caCell->GetHit1()[1];
819 0 : zcl[nDet] = caCell->GetHit1()[2];
820 0 : xerr[nDet] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
821 0 : yerr[nDet] = TMath::Sqrt(yTrErrDet*yTrErrDet+yTrErrMS*yTrErrMS);
822 0 : nDet++;
823 0 : }
824 0 : xcl[nDet] = caCell->GetHit2()[0];
825 0 : ycl[nDet] = caCell->GetHit2()[1];
826 0 : zcl[nDet] = caCell->GetHit2()[2];
827 0 : xerr[nDet] = TMath::Sqrt(xTrErrDet*xTrErrDet+xTrErrMS*xTrErrMS);
828 0 : yerr[nDet] = TMath::Sqrt(yTrErrDet*yTrErrDet+yTrErrMS*yTrErrMS);
829 0 : nDet++;
830 :
831 : }
832 :
833 :
834 :
835 :
836 : // y=a*x+b
837 : //
838 : // const Int_t nMaxh = 100;
839 : // Double_t xCl[nMaxh], yCl[nMaxh], yErr[nMaxh];
840 : // Int_t idet = 0;
841 : // for (Int_t i = 0; i < nDet; i++) {
842 : // if (i == skip) continue;
843 : // xCl[idet] = xcl[i];
844 : // yCl[idet] = ycl[i];
845 : // yErr[idet] = yerr[i];
846 : // idet++;
847 : // }
848 : //
849 : // Double_t S1, SXY, SX, SY, SXX, SsXY, SsXX, SsYY, Xm, Ym, s, delta, difx;
850 : //
851 : // S1 = SXY = SX = SY = SXX = 0.0;
852 : // SsXX = SsYY = SsXY = Xm = Ym = 0.;
853 : // difx = 0.;
854 : // for (Int_t i = 0; i < idet; i++) {
855 : // S1 += 1.0/(yErr[i]*yErr[i]);
856 : // SXY += xCl[i]*yCl[i]/(yErr[i]*yErr[i]);
857 : // SX += xCl[i]/(yErr[i]*yErr[i]);
858 : // SY += yCl[i]/(yErr[i]*yErr[i]);
859 : // SXX += xCl[i]*xCl[i]/(yErr[i]*yErr[i]);
860 : // if (i > 0) difx += TMath::Abs(xCl[i]-xCl[i-1]);
861 : // Xm += xCl[i];
862 : // Ym += yCl[i];
863 : // SsXX += xCl[i]*xCl[i];
864 : // SsYY += yCl[i]*yCl[i];
865 : // SsXY += xCl[i]*yCl[i];
866 : // }
867 : // delta = SXX*S1 - SX*SX;
868 : // if (delta == 0.) {
869 : // return kFALSE;
870 : // }
871 : // a = (SXY*S1 - SX*SY)/delta;
872 : // b = (SY*SXX - SX*SXY)/delta;
873 : //
874 : // Ym /= (Double_t)idet;
875 : // Xm /= (Double_t)idet;
876 : // SsYY -= (Double_t)idet*(Ym*Ym);
877 : // SsXX -= (Double_t)idet*(Xm*Xm);
878 : // SsXY -= (Double_t)idet*(Ym*Xm);
879 : // Double_t eps = 1.E-24;
880 : // if ((idet > 2) && (TMath::Abs(difx) > eps) && ((SsYY-(SsXY*SsXY)/SsXX) > 0.)) {
881 : // s = TMath::Sqrt((SsYY-(SsXY*SsXY)/SsXX)/(idet-2));
882 : // be = s*TMath::Sqrt(1./(Double_t)idet+(Xm*Xm)/SsXX);
883 : // ae = s/TMath::Sqrt(SsXX);
884 : // } else {
885 : // be = 0.;
886 : // ae = 0.;
887 : // }
888 : //
889 : return kTRUE;
890 :
891 0 : }
892 :
893 :
|
