Fixes for cmake
[u/mrichter/AliRoot.git] / PWG1 / AliMaterialBudget.cxx
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
52 using namespace std;
53
54 ClassImp(AliMaterialBudget)
55
56 //________________________________________________________________________
57 AliMaterialBudget::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
77 AliMaterialBudget::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 //________________________________________________________________________
101 AliMaterialBudget::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
127 AliMaterialBudget::~AliMaterialBudget(){
128   //
129   //
130   //
131   if (fDebugLevel>0)  printf("AliMaterialBudget::~AliMaterialBudget\n");
132   if (fDebugStreamer) delete fDebugStreamer;
133   fDebugStreamer=0;
134 }
135
136
137 //________________________________________________________________________
138 void 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 //________________________________________________________________________
175 void 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 //________________________________________________________________________
194 void 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 //________________________________________________________________________
222 void 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
237 TTreeSRedirector *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
257 AliExternalTrackParam * 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()};
266   Double_t charge = TMath::Nint(part->GetPDG()->Charge()/3.);
267   if (ref->X()*ref->Px()+ref->Y()*ref->Py() <0){
268     pxyz[0]*=-1;
269     pxyz[1]*=-1;
270     pxyz[2]*=-1;
271     charge*=-1.;
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
279 Bool_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
298 void  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 = 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
516 void 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());
534     seed.Rotate(alpha-seed.GetAlpha());
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
563 void  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;
584   Double_t covar[14];
585   for (Int_t icov=0; icov<14; icov++) covar[icov]=0;
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
657 void 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 = 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 = 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
885 Bool_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
917 AliTrackReference * 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
943 AliTrackReference *  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
968 void 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
981 void 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
1007 Bool_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