]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/TPC/AliMaterialBudget.cxx
-- remove obsolate classes
[u/mrichter/AliRoot.git] / PWG1 / TPC / AliMaterialBudget.cxx
CommitLineData
7cc34f08 1
2
3
4
5//
6//
7
8// This class estiamtes the material budget of the inner detectors in ALICE based
9// on the "upper/lower track matching"-method of the ALICE TPC.
10
11//
12//
13//
14
15//
16
17
18
19
20// ROOT includes
21#include <TChain.h>
22#include <TMath.h>
23#include <TVectorD.h>
24#include <TSystem.h>
25#include <TFile.h>
26#include "TGeoGlobalMagField.h"
27
28// ALIROOT includes
29#include <TTreeStream.h>
30#include <AliAnalysisManager.h>
31#include <AliESDInputHandler.h>
32#include "AliStack.h"
33#include "AliMCEvent.h"
34#include "AliMCEventHandler.h"
35
36#include <AliESD.h>
37#include "AliMaterialBudget.h"
38#include "AliGenInfoMaker.h"
39#include "AliHelix.h"
40
41//
42#include "AliMCInfo.h"
43#include "AliComparisonObject.h"
44#include "AliESDRecInfo.h"
45#include "AliTPCParamSR.h"
46#include "AliTracker.h"
47#include "AliTPCseed.h"
48
49// STL includes
50#include <iostream>
51
52using namespace std;
53
54ClassImp(AliMaterialBudget)
55
56//________________________________________________________________________
57AliMaterialBudget::AliMaterialBudget() :
58 AliAnalysisTask(),
59 fMCinfo(0), //! MC event handler
60 fESD(0),
61 fDebugStreamer(0),
62 fStreamLevel(0),
63 fDebugLevel(0),
64 fDebugOutputPath(),
65 fListHist(0),
66 fHistMult(0),
67 fCutMaxD(5), // maximal distance in rfi ditection
68 fCutMaxDz(40), // maximal distance in z ditection
69 fCutTheta(0.03), // maximal distan theta
70 fCutMinDir(-0.99) // direction vector products
71{
72 //
73 // Default constructor (should not be used)
74 //
75}
76
77AliMaterialBudget::AliMaterialBudget(const AliMaterialBudget& /*info*/) :
78 AliAnalysisTask(),
79 fMCinfo(0), //! MC event handler
80 fESD(0),
81 //
82 fDebugStreamer(0),
83 fStreamLevel(0),
84 fDebugLevel(),
85 fDebugOutputPath(),
86 fListHist(0),
87 fHistMult(0),
88 fCutMaxD(5), // maximal distance in rfi ditection
89 fCutMaxDz(40), // maximal distance in z ditection
90 fCutTheta(0.03), // maximal distan theta
91 fCutMinDir(-0.99) // direction vector products
92{
93 //
94 // Default constructor
95 //
96}
97
98
99
100//________________________________________________________________________
101AliMaterialBudget::AliMaterialBudget(const char *name) :
102 AliAnalysisTask(name, "AliMaterialBudget"),
103 fMCinfo(0), //! MC event handler
104 fESD(0),
105 fDebugStreamer(0),
106 fStreamLevel(0),
107 fDebugLevel(0),
108 fDebugOutputPath(),
109 fListHist(0),
110 fHistMult(0),
111 fCutMaxD(5), // maximal distance in rfi ditection
112 fCutMaxDz(40), // maximal distance in z ditection
113 fCutTheta(0.03), // maximal distan theta
114 fCutMinDir(-0.99) // direction vector products
115{
116 //
117 // Normal constructor
118 //
119 // Input slot #0 works with a TChain
120 DefineInput(0, TChain::Class());
121 // Output slot #0 writes into a TList
122 DefineOutput(0, TList::Class());
123 //
124 //
125}
126
127AliMaterialBudget::~AliMaterialBudget(){
128 //
129 //
130 //
131 if (fDebugLevel>0) printf("AliMaterialBudget::~AliMaterialBudget\n");
132 if (fDebugStreamer) delete fDebugStreamer;
133 fDebugStreamer=0;
134}
135
136
137//________________________________________________________________________
138void AliMaterialBudget::ConnectInputData(Option_t *)
139{
140 //
141 // Connect the input data
142 //
143 if(fDebugLevel>3)
144 cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
145
146 TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
147 if (!tree) {
148 //Printf("ERROR: Could not read chain from input slot 0");
149 }
150 else {
151 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
152 if (!esdH) {
153 //Printf("ERROR: Could not get ESDInputHandler");
154 }
155 else {
156 esdH->SetActiveBranches("ESDfriend");
157 fESD = esdH->GetEvent();
158 //Printf("*** CONNECTED NEW EVENT ****");
159 }
160 }
161 AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
162 mcinfo->SetReadTR(kTRUE);
163
164 fMCinfo = mcinfo->MCEvent();
165
166
167}
168
169
170
171
172
173
174//________________________________________________________________________
175void AliMaterialBudget::CreateOutputObjects()
176{
177 //
178 // Connect the output objects
179 //
180 if(fDebugLevel>3)
181 cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
182 //
183 fListHist = new TList();
184 fListHist->SetOwner();
185 //
186 fHistMult = new TH1F("HistMult", "Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
187 fListHist->Add(fHistMult);
188
189
190}
191
192
193//________________________________________________________________________
194void AliMaterialBudget::Exec(Option_t *) {
195 //
196 // Execute analysis for current event
197 //
198
199 if(fDebugLevel>3)
200 cout << "AliMaterialBudget::Exec()" << endl;
201
202 fHistMult->Fill(fESD->GetNumberOfTracks());
203 //FindPairs(fESD); // nearly everything takes place in find pairs...
204
205 // If MC has been connected
206
207 if (!fMCinfo){
208 cout << "Not MC info\n" << endl;
209 }else{
210 ProcessMCInfo();
211 //mcinfo->Print();
212 //DumpInfo();
213 }
214 //
215 PostData(0, fListHist);
216}
217
218
219
220
221//________________________________________________________________________
222void AliMaterialBudget::Terminate(Option_t *) {
223 //
224 // Terminate loop
225 //
226 if(fDebugLevel>3)
227 printf("AliMaterialBudget: Terminate() \n");
228 //
229 if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
230 if (fDebugStreamer) delete fDebugStreamer;
231 fDebugStreamer = 0;
232 return;
233}
234
235
236
237TTreeSRedirector *AliMaterialBudget::GetDebugStreamer(){
238 //
239 // Get Debug streamer
240 // In case debug streamer not yet initialized and StreamLevel>0 create new one
241 //
242 if (fStreamLevel==0) return 0;
243 if (fDebugStreamer) return fDebugStreamer;
244 TString dsName;
245 dsName=GetName();
246 dsName+="Debug.root";
247 dsName.ReplaceAll(" ","");
248 fDebugStreamer = new TTreeSRedirector(dsName.Data());
249 return fDebugStreamer;
250}
251
252
253
254
255
256
257AliExternalTrackParam * AliMaterialBudget::MakeTrack(const AliTrackReference* ref, TParticle*part)
258{
259 //
260 // Make track out of the track ref
261 // part - TParticle used to determine chargr
262 // the covariance matrix - equal 0 - starting from ideal MC position
263 if (!part->GetPDG()) return 0;
264 Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
265 Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
c6d45264 266 Int_t charge = TMath::Nint(part->GetPDG()->Charge()/3.);
7cc34f08 267 if (ref->X()*ref->Px()+ref->Y()*ref->Py() <0){
268 pxyz[0]*=-1;
269 pxyz[1]*=-1;
270 pxyz[2]*=-1;
c6d45264 271 charge*=-1;
7cc34f08 272 }
273 Double_t cv[21];
274 for (Int_t i=0; i<21;i++) cv[i]=0;
275 AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,charge);
276 return param;
277}
278
279Bool_t AliMaterialBudget::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
280 //
281 // Propagate track to point xyz using
282 // AliTracker::PropagateTo functionality
283 //
284 // param - track parameters
285 // xyz - position to propagate
286 // mass - particle mass
287 // step - step to be used
288 Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
289 AliTracker::PropagateTrackToBxByBz(param, radius+step, mass, step, kTRUE,0.99);
290 AliTracker::PropagateTrackToBxByBz(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
291 Double_t sxyz[3]={0,0,0};
292 AliESDVertex vertex(xyz,sxyz);
293 Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
294 return isOK;
295}
296
297
298void AliMaterialBudget::ProcessMCInfo(){
299 //
300 //
301 //
302 //
303 TParticle * particle= new TParticle;
304 TClonesArray * trefs = new TClonesArray("AliTrackReference");
305 const Double_t kPcut=0.05;
306 const Double_t kMinDrITS = 2.; // minimal distance between references
307 const Double_t kMinDrTRD = 8.; // minimal distance between references
308 const Double_t kMinDrTOF = 10.; // minimal distance between references
309 //
310 //
311 // Process tracks
312 //
313 Int_t npart = fMCinfo->GetNumberOfTracks();
314 if (npart==0) return;
315 Double_t vertex[4]={0,0,0,0};
316 fMCinfo->GetParticleAndTR(0, particle, trefs);
317 if (particle){
318 vertex[0]=particle->Vx();
319 vertex[1]=particle->Vy();
320 vertex[2]=particle->Vz();
321 vertex[3]=particle->R();
322 }
323 //
324 //
325 AliTrackReference dummy,*pdummy= &dummy;
326 AliExternalTrackParam epdummy,*pepdummy= &epdummy;
327 Int_t nRefITS =0;
328 Int_t nRefTPC =0;
329 Int_t nRefTRD =0;
330 Int_t nRefTOF =0;
331 AliTrackReference * refITS0, *refITS1;
332 AliTrackReference * refTPC0, *refTPC1;
333 AliTrackReference * refTPCIn0, *refTPCIn1;
334 AliTrackReference * refTRD0, *refTRD1;
335 AliTrackReference * refTOF0, *refTOF1;
336 AliTrackReference *refMinR;
337 //
338 for (Int_t ipart=0;ipart<npart;ipart++){
339 //Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
340 AliMCParticle * pp = (AliMCParticle*) fMCinfo->GetTrack(ipart);
341 if (!pp) continue;
342 if (particle->P()<kPcut) continue;
343 Double_t mass = particle->GetMass();
344
345 // RESET
346 nRefITS =0;
347 nRefTPC =0;
348 nRefTRD =0;
349 nRefTOF =0;
350 refITS0=pdummy; refITS1=pdummy;
351 refTPC0=pdummy; refTPC1=pdummy;
352 refTPCIn0=pdummy; refTPCIn1=pdummy;
353 refTRD0=pdummy; refTRD1=pdummy;
354 refTOF0=pdummy; refTOF1=pdummy;
355 refMinR = pdummy;
356 //
357 Int_t nref = pp->GetNumberOfTrackReferences();
358 if (nref==0) continue;
359 for (Int_t iref = 0; iref < nref; iref++) {
360 AliTrackReference *ref = pp->GetTrackReference(iref);
361 if (ref->DetectorId()==AliTrackReference::kDisappeared) continue;
362 // if (ref.Px()*particle.Px()+ref.Py()*particle.Py()<0) break; // returning track
363 //
364 if (ref->DetectorId()==AliTrackReference::kITS){
365 if (TMath::Abs(ref->R()-refITS1->R())>kMinDrITS) {
366 refITS1 = ref;
367 nRefITS++;
368 }
369 if (refITS0==pdummy) refITS0=ref;
370 }
371 if (ref->DetectorId()==AliTrackReference::kTPC){
372 nRefTPC++;
373 refTPC1 = ref;
374 if (refTPC0==pdummy) refTPC0=ref;
375 }
376 if (ref->DetectorId()==AliTrackReference::kTRD){
377 if (TMath::Abs(ref->R()-refTRD1->R())>kMinDrTRD) {
378 refTRD1 = ref;
379 nRefTRD++;
380 }
381 if (refTRD0==pdummy) refTRD0=ref;
382 }
383 if (ref->DetectorId()==AliTrackReference::kTOF){
384 if (TMath::Abs(ref->X()-refTOF1->X()) + TMath::Abs(ref->Y()-refTOF1->Y()) + TMath::Abs(ref->Z()-refTOF1->Z())>kMinDrTOF) {
385 refTOF1 = ref;
386 nRefTOF++;
387 }
388 if (refTOF0==pdummy) refTOF0=ref;
389 }
390 //
391 // "find inner track ref"
392 if (ref->DetectorId()==AliTrackReference::kTPC){
393 if (ref->Px()*ref->X()+ref->Py()*ref->Y()<0){
394 // track in
395 if (refTPCIn0 == pdummy) refTPCIn0=ref;
396 if (refTPCIn0 != pdummy && refTPCIn0->R()>ref->R())
397 refTPCIn0=ref;
398 }
399 if (ref->Px()*ref->X()+ref->Py()*ref->Y()>0){
400 // track in
401 if (refTPCIn1 == pdummy) refTPCIn1=ref;
402 if (refTPCIn1 != pdummy && refTPCIn1->R()>ref->R())
403 refTPCIn1=ref;
404 }
405 }
406
407
408 if (refMinR==pdummy && ref->P()>0 ){
409 refMinR=ref;
410 }
411 if (refMinR->R()>ref->R() && ref->P()>0 ) refMinR=ref;
412
413 }
414 //
415 AliExternalTrackParam * trackMC = pepdummy;
416 //track0->GetDZ(0,0,0,bz,dvertex0)
417 Float_t dist[2]={0,0};
418 AliMagF* field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
419 Double_t esdfield= fESD->GetMagneticField();
420 Double_t xyz[3]={0,0,0};
421 Double_t bxyz[3]={0,0,0};
422 field->Field(xyz,bxyz);
423 if (refMinR->P()>0) {
424 trackMC = MakeTrack(refMinR,particle);
425 trackMC->GetDZ(0,0,0,bxyz[2],dist);
426 }
427 Double_t alphaTOF0 = TMath::ATan2(refTOF0->Y(),refTOF0->X());
428 Double_t alphaTOF1 = TMath::ATan2(refTOF1->Y(),refTOF1->X());
429 Int_t dsecTOF = TMath::Nint(180*(alphaTOF0-alphaTOF1)/(TMath::Pi()*20.)-0.5);
430 //
431 // make the two different TPC tracks and propagate them to their DCA
432 //
433 Double_t dP = 0;
434 Bool_t statusProp = kFALSE;
435 Double_t dY = 0;
436 Double_t dZ = 0;
437 AliExternalTrackParam * track0 = pepdummy;
438 AliExternalTrackParam * track1 = pepdummy;
439 AliExternalTrackParam * otrack0 = pepdummy;
440 AliExternalTrackParam * otrack1 = pepdummy;
441 if (refTPCIn0!=pdummy && refTPCIn1!=pdummy) {
442 track0 = MakeTrack(refTPCIn0,particle);
443 track1 = MakeTrack(refTPCIn1,particle);
444 otrack0 = MakeTrack(refTPCIn0,particle);
445 otrack1 = MakeTrack(refTPCIn1,particle);
446 dP = track0->P() - track1->P(); // momentum loss
447 statusProp = AliMaterialBudget::PropagateCosmicToDCA(track0,track1,mass);
448 if (statusProp) {
449 dY = track0->GetY() - track1->GetY();
450 dZ = track0->GetZ() - track1->GetZ();
451 }
452 }
453 //
454 TTreeSRedirector *pcstream = GetDebugStreamer();
455 if (pcstream){
456 char name[100];
457 for (Int_t id=0; id<3;id++){
458
459 // if (id==0) sprintf(name,"mcAll"); // all tracks: inconvenient to cut if on is only interest in tracks which reach the TPC
460 if (id==0) continue; // require TPC
461 if (id==1) sprintf(name,"mcITS");
462 if (id==2) sprintf(name,"mcTPC");
463 if (id==1&& nRefITS==0) continue;
464 if (id==2&& nRefTPC==0) continue;
465
466 (*pcstream)<<name<<
467 "ipart="<<ipart<<
468 "p.="<<particle<<
469 "mass="<<mass<<
470 "tbfield="<<bxyz[2]<<
471 "esdbfield="<<esdfield<<
472 // counter
473 "nref="<<nref<<
474 "nRefITS="<<nRefITS<<
475 "nRefTPC="<<nRefTPC<<
476 "nRefTRD="<<nRefTRD<<
477 "nRefTOF="<<nRefTOF<<
478 //references
479 "refMinR.="<<refMinR<<
480 "refITS0.="<<refITS0<<
481 "refITS1.="<<refITS1<<
482 "refTPC0.="<<refTPC0<<
483 "refTPC1.="<<refTPC1<<
484 "refTPCIn0.="<<refTPCIn0<<
485 "refTPCIn1.="<<refTPCIn1<<
486 "refTRD0.="<<refTRD0<<
487 "refTRD1.="<<refTRD1<<
488 "refTOF0.="<<refTOF0<<
489 "refTOF1.="<<refTOF1<<
490 //trigger variables
491 "dsecTOF="<<dsecTOF<< // delta TOF sectors
492 "aTOF0="<<alphaTOF0<<
493 "aTOF1="<<alphaTOF1<<
494 //
495 // track
496 "dr="<<dist[0]<<
497 "dz="<<dist[1]<<
498 //
499 // "two" TPC tracks
500 "status="<<statusProp<<
501 "dP="<<dP<<
502 "otrack0.="<<otrack0<<
503 "otrack1.="<<otrack1<<
504 "track0.="<<track0<<
505 "track1.="<<track1<<
506 "dY="<<dY<<
507 "dZ="<<dZ<<
508 "\n";
509 }
510 }
511 }
512}
513
514
515
516void AliMaterialBudget::ProcessRefTracker(AliTrackReference* refIn, AliTrackReference* refOut, TParticle*part,Int_t type){
517 //
518 // Test propagation from In to out
519 //
520 AliExternalTrackParam *param = 0;
521 AliExternalTrackParam *paramMC = 0;
522 Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
523 Double_t mass = part->GetMass();
524 Double_t step=1;
525 //
526 param=MakeTrack(refOut,part);
527 paramMC=MakeTrack(refOut,part);
528 if (!param) return;
529 if (type<3) PropagateToPoint(param,xyzIn, mass, step);
530 if (type==3) {
531 AliTPCseed seed;
532 seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
533 Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
26a35193 534 if(seed.Rotate(alpha-seed.GetAlpha()) == kFALSE) return;
7cc34f08 535 seed.SetMass(mass);
536 for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
537 seed.PropagateTo(xlayer);
538 }
539 seed.PropagateTo(refIn->R());
540 param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
541 }
542 TTreeSRedirector *pcstream = GetDebugStreamer();
543 TVectorD gpos(3);
544 TVectorD gmom(3);
545 param->GetXYZ(gpos.GetMatrixArray());
546 param->GetPxPyPz(gmom.GetMatrixArray());
547 if (pcstream){
548 (*pcstream)<<"MC"<<
549 "type="<<type<<
550 "step="<<step<<
551 "refIn.="<<refIn<<
552 "refOut.="<<refOut<<
553 "p.="<<part<<
554 "par.="<<param<<
555 "parMC.="<<paramMC<<
556 "gpos.="<<&gpos<<
557 "gmom.="<<&gmom<<
558 "\n";
559 }
560}
561
562
563void AliMaterialBudget::FitTrackRefs(TParticle * part, TClonesArray * trefs){
564 //
565 //
566 //
567 //
568 const Int_t kMinRefs=6;
569 Int_t nrefs = trefs->GetEntries();
570 if (nrefs<kMinRefs) return; // we should have enough references
571 Int_t iref0 =-1;
572 Int_t iref1 =-1;
573
574 for (Int_t iref=0; iref<nrefs; iref++){
575 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
576 if (!ref) continue;
577 Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
578 if (dir<0) break;
579 if (ref->DetectorId()!=AliTrackReference::kTPC) continue;
580 if (iref0<0) iref0 = iref;
581 iref1 = iref;
582 }
583 if (iref1-iref0<kMinRefs) return;
c6d45264 584 Double_t covar[15];
585 for (Int_t icov=0; icov<15; icov++) covar[icov]=0;
7cc34f08 586 covar[0]=1;
587 covar[2]=1;
588 covar[5]=1;
589 covar[9]=1;
590 covar[14]=1;
591
592 AliTrackReference * refIn = (AliTrackReference*)trefs->At(iref0);
593 AliTrackReference * refOut = (AliTrackReference*)trefs->At(iref1);
594 AliExternalTrackParam *paramPropagate= MakeTrack(refIn,part);
595 AliExternalTrackParam *paramUpdate = MakeTrack(refIn,part);
596 paramUpdate->AddCovariance(covar);
597 Double_t mass = part->GetMass();
598 Double_t charge = part->GetPDG()->Charge()/3.;
599/*
600 Float_t alphaIn= TMath::ATan2(refIn->Y(),refIn->X());
601 Float_t radiusIn= refIn->R();
602 Float_t alphaOut= TMath::ATan2(refOut->Y(),refOut->X());
603 Float_t radiusOut= refOut->R();
604*/
605 Bool_t isOKP=kTRUE;
606 Bool_t isOKU=kTRUE;
607 AliMagF * field = (AliMagF*) TGeoGlobalMagField::Instance()->GetField();
608 for (Int_t iref = iref0; iref<=iref1; iref++){
609 AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
610 Float_t alphaC= TMath::ATan2(ref->Y(),ref->X());
611 Double_t pos[3] = {ref->X(), ref->Y(), ref->Z()};
612 Double_t mag[3];
613 field->Field(pos,mag);
614 isOKP&=paramPropagate->Rotate(alphaC);
615 isOKU&=paramUpdate->Rotate(alphaC);
616 for (Float_t xref= paramPropagate->GetX(); xref<ref->R(); xref++){
617 isOKP&=paramPropagate->PropagateTo(xref, mag[2]);
618 isOKU&=paramUpdate->PropagateTo(xref, mag[2]);
619 }
620 isOKP&=paramPropagate->PropagateTo(ref->R(), mag[2]);
621 isOKU&=paramUpdate->PropagateTo(ref->R(), mag[2]);
622 Double_t clpos[2] = {0, ref->Z()};
623 Double_t clcov[3] = { 0.005,0,0.005};
624 isOKU&= paramUpdate->Update(clpos, clcov);
625 }
626 TTreeSRedirector *pcstream = GetDebugStreamer();
627 if (pcstream){
628 TVectorD gposU(3);
629 TVectorD gmomU(3);
630 TVectorD gposP(3);
631 TVectorD gmomP(3);
632 paramUpdate->GetXYZ(gposU.GetMatrixArray());
633 paramUpdate->GetPxPyPz(gmomU.GetMatrixArray());
634 paramPropagate->GetXYZ(gposP.GetMatrixArray());
635 paramPropagate->GetPxPyPz(gmomP.GetMatrixArray());
636
637 (*pcstream)<<"MCupdate"<<
638 "isOKU="<<isOKU<<
639 "isOKP="<<isOKP<<
640 "m="<<mass<<
641 "q="<<charge<<
642 "part.="<<part<<
643 "refIn.="<<refIn<<
644 "refOut.="<<refOut<<
645 "pP.="<<paramPropagate<<
646 "pU.="<<paramUpdate<<
647 "gposU.="<<&gposU<<
648 "gmomU.="<<&gmomU<<
649 "gposP.="<<&gposP<<
650 "gmomP.="<<&gmomP<<
651 "\n";
652 }
653}
654
655
656
657void AliMaterialBudget::FindPairs(AliESDEvent * event) {
658 //
659 // This function matches the cosmic tracks and calculates the energy loss in the material.
660 // If accessible the "true" energy loss is determined with the MC track references.
661 //
662 //
663 // Find cosmic pairs
664 //
665 // Track0 is choosen in upper TPC part
666 // Track1 is choosen in lower TPC part
667 //
668 if(fDebugLevel>3)
669 cout << "AliMaterialBudget::FindPairs()" << endl;
670
671
672 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
673 Int_t ntracks=event->GetNumberOfTracks();
674 TObjArray tpcSeeds(ntracks);
675 if (ntracks==0) return;
676 Double_t vtxx[3]={0,0,0};
677 Double_t svtxx[3]={0.000001,0.000001,100.};
678 AliESDVertex vtx(vtxx,svtxx);
679 //
680 //track loop
681 //
682 for (Int_t i=0;i<ntracks;++i) {
683 AliESDtrack *track = event->GetTrack(i);
684 const AliExternalTrackParam * trackIn = track->GetInnerParam();
685 const AliExternalTrackParam * trackOut = track->GetOuterParam();
686 if (!trackIn) continue;
687 if (!trackOut) continue;
688 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
689 TObject *calibObject;
690 AliTPCseed *seed = 0;
691 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
692 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
693 }
694 if (seed) tpcSeeds.AddAt(seed,i);
695 }
696
697 if (ntracks<2) return;
698 //
699 // Find pairs
700 //
701 for (Int_t i=0;i<ntracks;++i) {
702 AliESDtrack *track0 = event->GetTrack(i);
703 // track0 - choosen upper part
704 if (!track0) continue;
705 if (!track0->GetOuterParam()) continue;
706 if (track0->GetOuterParam()->GetAlpha()<0) continue;
707 Double_t dir0[3];
708 track0->GetDirection(dir0);
709 for (Int_t j=0;j<ntracks;++j) {
710 if (i==j) continue;
711 AliESDtrack *track1 = event->GetTrack(j);
712 //track 1 lower part
713 if (!track1) continue;
714 if (!track1->GetOuterParam()) continue;
715 if (track1->GetOuterParam()->GetAlpha()>0) continue;
716 //
717 Double_t dir1[3];
718 track1->GetDirection(dir1);
719
720 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
721 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
722 if (! seed0) continue;
723 if (! seed1) continue;
724 //
725 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
726 Float_t d0 = track0->GetLinearD(0,0);
727 Float_t d1 = track1->GetLinearD(0,0);
728 //
729 // conservative cuts - convergence to be guarantied
730 // applying before track propagation
731 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
732 if (dir>fCutMinDir) continue; // direction vector product
733 Float_t bz = AliTracker::GetBz();
734 Float_t dvertex0[2]; //distance to 0,0
735 Float_t dvertex1[2]; //distance to 0,0
736 track0->GetDZ(0,0,0,bz,dvertex0);
737 track1->GetDZ(0,0,0,bz,dvertex1);
738 if (TMath::Abs(dvertex0[1])>250) continue;
739 if (TMath::Abs(dvertex1[1])>250) continue;
740 //
741 //
742 //
743 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
744 AliExternalTrackParam param0(*track0);
745 AliExternalTrackParam param1(*track1);
746 //
747 // Propagate using Magnetic field and correct fo material budget
748 //
749 AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,0.0005,3,kTRUE);
750 AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,0.0005,3,kTRUE);
751 //
752 // Propagate rest to the 0,0 DCA - z should be ignored
753 //
754 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
755 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
756 //
757 param0.GetDZ(0,0,0,bz,dvertex0);
758 param1.GetDZ(0,0,0,bz,dvertex1);
759 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
760 //
761 Double_t xyz0[3];//,pxyz0[3];
762 Double_t xyz1[3];//,pxyz1[3];
763 param0.GetXYZ(xyz0);
764 param1.GetXYZ(xyz1);
765 Bool_t isPair = IsPair(&param0,&param1);
766 //
767 // HERE WE WILL PUT THE ACCESS TO THE MC TRACKS AND MATCH THESE !!!!
768 //
769 Int_t label0 = TMath::Abs(track0->GetLabel());
770 AliMCParticle *mcParticle0 = (AliMCParticle*) fMCinfo->GetTrack(label0);
771 TParticle *particle0 = mcParticle0->Particle();
772 AliTrackReference *ref0 = GetFirstTPCTrackRef(mcParticle0); // get the first TPC track reference
773 if (!ref0) continue;
774 AliExternalTrackParam *paramMC0 = 0;
775 paramMC0 = MakeTrack(ref0, particle0);
776 //
777 Int_t label1 = TMath::Abs(track1->GetLabel());
778 AliMCParticle *mcParticle1 = (AliMCParticle*) fMCinfo->GetTrack(label1);
779 TParticle *particle1 = mcParticle1->Particle();
780 AliTrackReference *ref1 = GetFirstTPCTrackRef(mcParticle1); // get the first TPC track reference
781 if (!ref1) continue;
782 AliExternalTrackParam *paramMC1 = 0;
783 paramMC1 = MakeTrack(ref1, particle1);
784 //
785 // ACCESS TOF INFORMATION
786 Int_t nTrackRefTOF0 = 0;
787 Int_t nTrackRefITS0 = 0;
788 AliTrackReference * refLastTOF0 = 0;
789 AliTrackReference * refFirstTOF0 = GetAllTOFinfo(mcParticle0, nTrackRefTOF0, nTrackRefITS0);
790 Float_t alphaTOF0 = 0;
791 if (refFirstTOF0) alphaTOF0 = refFirstTOF0->Alpha();
792 //
793 Int_t nTrackRefTOF1 = 0;
794 Int_t nTrackRefITS1 = 0;
795 AliTrackReference * refLastTOF1 = 0;
796 AliTrackReference * refFirstTOF1 =GetAllTOFinfo(mcParticle1, nTrackRefTOF1, nTrackRefITS1);
797 Float_t alphaTOF1 = 0;
798 if (refFirstTOF1) alphaTOF1 = refFirstTOF1->Alpha();
799 //cout << " STATUS "<<nTrackRefTOF0<<" "<<refFirstTOF0<<" "<<refLastTOF0<<" " <<refFirstTOF1<<" "<<refLastTOF1<<endl;
800 //
801 //
802 //
803 if (fStreamLevel>0){
804 TTreeSRedirector * cstream = GetDebugStreamer();
805 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
806 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
807 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
808 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
809 //
810 //
811 //
812 if (cstream) {
813 (*cstream) << "Track0" <<
814 "dir="<<dir<< // direction
815 "OK="<<isPair<< // will be accepted
816 "b0="<<b0<< // propagate status
817 "b1="<<b1<< // propagate status
818 //
819 "Particle.="<<particle0<< // TParticle with generated momentum
820 "NTrackRefTOF0="<<nTrackRefTOF0<< // Number of TOF track references upper
821 "NTrackRefTOF1="<<nTrackRefTOF1<< // Number of TOF track references lower
822 "NTrackRefITS0="<<nTrackRefITS0<< // Number of ITS track references upper
823 "NTrackRefITS1="<<nTrackRefITS1<< // Number of ITS track references lower
824 "Alpha0="<<alphaTOF0<< // alpha upper
825 "Alpha1="<<alphaTOF1<< // alpha lower
826 "RefFirstTOF0.="<<refFirstTOF0<< // first tof reference upper
827 "RefLastTOF0.="<<refLastTOF0<< // last tof reference upper
828 "RefFirstTOF1.="<<refFirstTOF1<< // first tof reference lower
829 "RefLastTOF1.="<<refLastTOF1<< // last tof reference lower
830 //
831 "Orig0.=" << track0 << // original track 0
832 "Orig1.=" << track1 << // original track 1
833 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
834 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
835 "Ip0.="<<ip0<< // inner param - upper
836 "Ip1.="<<ip1<< // inner param - lower
837 "Op0.="<<op0<< // outer param - upper
838 "Op1.="<<op1<< // outer param - lower
839 //
840 "paramTrackRef0.="<<paramMC0<< // "ideal" MC true track parameters from track references - upper
841 "paramTrackRef1.="<<paramMC1<< // "ideal" MC true track parameters from track references - lower
842 //
843 "v00="<<dvertex0[0]<< // distance using kalman
844 "v01="<<dvertex0[1]<< //
845 "v10="<<dvertex1[0]<< //
846 "v11="<<dvertex1[1]<< //
847 "d0="<<d0<< // linear distance to 0,0
848 "d1="<<d1<< // linear distance to 0,0
849 //
850 //
851 //
852 "x00="<<xyz0[0]<< // global position close to vertex
853 "x01="<<xyz0[1]<<
854 "x02="<<xyz0[2]<<
855 //
856 "x10="<<xyz1[0]<< // global position close to vertex
857 "x11="<<xyz1[1]<<
858 "x12="<<xyz1[2]<<
859 //
860 "dir00="<<dir0[0]<< // direction upper
861 "dir01="<<dir0[1]<<
862 "dir02="<<dir0[2]<<
863 //
864 "dir10="<<dir1[0]<< // direction lower
865 "dir11="<<dir1[1]<<
866 "dir12="<<dir1[2]<<
867 //
868 //
869 "Seed0.=" << seed0 << // original seed 0
870 "Seed1.=" << seed1 << // original seed 1
871 //
872 "\n";
873 }
874 }
875 delete paramMC0;
876 delete paramMC1;
877 }
878 }
879
880 return;
881
882
883}
884
885Bool_t AliMaterialBudget::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
886 //
887 //
888 /*
889 // 0. Same direction - OPOSITE - cutDir +cutT
890 TCut cutDir("cutDir","dir<-0.99")
891 // 1.
892 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
893 //
894 // 2. The same rphi
895 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
896 //
897 //
898 //
899 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
900 // 1/Pt diff cut
901 */
902 const Double_t *p0 = tr0->GetParameter();
903 const Double_t *p1 = tr1->GetParameter();
904 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
905 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
906 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
907
908 Double_t d0[3], d1[3];
909 tr0->GetDirection(d0);
910 tr1->GetDirection(d1);
911 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
912 //
913 return kTRUE;
914}
915
916
917AliTrackReference * AliMaterialBudget::GetFirstTPCTrackRef(AliMCParticle *mcParticle)
918{
919 // return first TPC track reference
920 if(!mcParticle) return 0;
921
922 // find first track reference
923 // check direction to select proper reference point for looping tracks
924 Int_t nTrackRef = mcParticle->GetNumberOfTrackReferences();
925 AliTrackReference *ref = 0;
926 AliTrackReference *refIn = 0;
927 for (Int_t iref = 0; iref < nTrackRef; iref++) {
928 ref = mcParticle->GetTrackReference(iref);
929 if(ref && (ref->DetectorId()==AliTrackReference::kTPC))
930 {
931 //Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
932 //if(dir < 0.) break;
933
934 refIn = ref;
935 break;
936 }
937 }
938
939 return refIn;
940}
941
942
943AliTrackReference * AliMaterialBudget::GetAllTOFinfo(AliMCParticle *mcParticle, Int_t &nTrackRef, Int_t &nTrackRefITS, Int_t retValue) {
944 //
945 //
946 //
947
948 if(!mcParticle) return 0;
949 Int_t counter = 0;
950 nTrackRef = 0;
951 nTrackRefITS = 0;
952 AliTrackReference *ref = 0;
953 for (Int_t iref = 0; iref < mcParticle->GetNumberOfTrackReferences(); iref++) {
954 ref = mcParticle->GetTrackReference(iref);
955 if(ref && (ref->DetectorId()==AliTrackReference::kTOF)) {
956 counter = iref;
957 nTrackRef++;
958 }
959 if(ref && (ref->DetectorId()==AliTrackReference::kITS)) nTrackRefITS++;
960 }
961 if (nTrackRef ==0) return 0;
962 if (retValue == 0) return mcParticle->GetTrackReference(counter - nTrackRef +1);
963 return mcParticle->GetTrackReference(counter);
964
965}
966
967
968void AliMaterialBudget::FinishTaskOutput()
969{
970 //
971 // According description in AliAnalisysTask this method is call
972 // on the slaves before sending data
973 //
974 Terminate("slave");
975 gSystem->Exec("pwd");
976 RegisterDebugOutput();
977
978}
979
980
981void AliMaterialBudget::RegisterDebugOutput(){
982 //
983 //
984 //
985 //
986 // store - copy debug output to the destination position
987 // currently ONLY for local copy
988 TString dsName;
989 dsName=GetName();
990 dsName+="Debug.root";
991 dsName.ReplaceAll(" ","");
992 TString dsName2=fDebugOutputPath.Data();
993 gSystem->MakeDirectory(dsName2.Data());
994 dsName2+=gSystem->HostName();
995 gSystem->MakeDirectory(dsName2.Data());
996 dsName2+="/";
997 dsName2+=gSystem->BaseName(gSystem->pwd());
998 dsName2+="/";
999 gSystem->MakeDirectory(dsName2.Data());
1000 dsName2+=dsName;
1001 AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
1002 printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
1003 TFile::Cp(dsName.Data(),dsName2.Data());
1004}
1005
1006
1007Bool_t AliMaterialBudget::PropagateCosmicToDCA(AliExternalTrackParam *param0, AliExternalTrackParam *param1, Double_t mass){
1008 //
1009 // param0 - upper part of cosmic track
1010 // param1 - lower part of cosmic track
1011 //
1012 // 0. Propagate both tracks to DCA to (0,0,0)
1013 // 1. After propagation to DCA rotate track param1 to corrdinate system of track1 <-> rotate param0 to coordinate system of param 1 ????
1014 // 2. Propagate track 1 to refernce x from track0
1015 //
1016
1017 // step 0.
1018
1019 Float_t d0 = param0->GetLinearD(0,0);
1020 Float_t d1 = param1->GetLinearD(0,0);
1021 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
1022 //
1023 // propagate in the beginning taking all material into account
1024 //
1025 AliTracker::PropagateTrackToBxByBz(param0,dmax+1.,mass,0.5,kTRUE,0.99,-1.);
1026 AliTracker::PropagateTrackToBxByBz(param1,dmax+1.,mass,0.5,kTRUE,0.99,1.);
1027 //
1028 Double_t vtxx[3]={0,0,0};
1029 Double_t svtxx[3]={0.000001,0.000001,100.};
1030 AliESDVertex vtx(vtxx,svtxx);
1031 //
1032 Bool_t b0 = param0->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1033 Bool_t b1 = param1->PropagateToDCA(&vtx,AliTracker::GetBz(),1000);
1034
1035 if (!(b0 && b1)) return 0;
1036
1037 // step 1.
1038
1039 Float_t dAlpha = param0->GetAlpha();
1040 param1->Rotate(dAlpha);
1041
1042 // step 2.
1043
1044 Float_t refX = param0->GetX();
1045 param1->PropagateTo(refX,AliTracker::GetBz());
1046
1047 return kTRUE;
1048
1049
1050}
1051
1052