remove obselete clean QA macros command
[u/mrichter/AliRoot.git] / PWG1 / AliMCTrackingTestTask.cxx
1 //
2 // This class is the task to check the 
3 // Propagation method used in the 
4 //               1. AliExternalTrackParam 
5 //               2. AliTracker
6 //
7 // Principle - Creates AliExternalTrackParam form 1 Track Refernece - 
8 //             Propagate it to other
9 // Magnetic filed and the geomtry has to bec 
10
11 //
12 // ROOT includes
13 #include <TChain.h>
14 #include <TMath.h>
15 #include <TVectorD.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18
19 // ALIROOT includes
20 #include <TTreeStream.h>
21 #include <AliAnalysisManager.h>
22 #include <AliESDInputHandler.h>
23 #include "AliStack.h"
24 #include "AliMCEvent.h"
25 #include "AliMCEventHandler.h"
26
27 #include <AliESD.h>
28 #include "AliMCTrackingTestTask.h"
29 #include "AliGenInfoMaker.h"
30 #include "AliHelix.h"
31
32 //
33 #include "AliMCInfo.h"
34 #include "AliComparisonObject.h"
35 #include "AliESDRecInfo.h"
36 #include "AliTPCParamSR.h"
37 #include "AliTracker.h"
38 #include "AliTPCseed.h"
39
40 // STL includes
41 #include <iostream>
42
43 using namespace std;
44
45 ClassImp(AliMCTrackingTestTask)
46
47 //________________________________________________________________________
48 AliMCTrackingTestTask::AliMCTrackingTestTask() : 
49   AliAnalysisTask(), 
50   fMCinfo(0),     //! MC event handler
51   fESD(0),
52   fDebugStreamer(0),
53   fStreamLevel(0),
54   fDebugLevel(0),
55   fDebugOutputPath()
56 {
57   //
58   // Default constructor (should not be used)
59   //
60 }
61
62 AliMCTrackingTestTask::AliMCTrackingTestTask(const AliMCTrackingTestTask& /*info*/) : 
63   AliAnalysisTask(), 
64   fMCinfo(0),     //! MC event handler
65   fESD(0),
66   //
67   fDebugStreamer(0),
68   fStreamLevel(0),
69   fDebugLevel(),
70   fDebugOutputPath()
71 {
72   //
73   // Default constructor 
74   //
75 }
76
77
78
79 //________________________________________________________________________
80 AliMCTrackingTestTask::AliMCTrackingTestTask(const char *name) : 
81   AliAnalysisTask(name, "AliMCTrackingTestTask"), 
82   fMCinfo(0),     //! MC event handler
83   fESD(0),
84   fDebugStreamer(0),
85   fStreamLevel(0),
86   fDebugLevel(0),
87   fDebugOutputPath()
88 {
89   //
90   // Normal constructor
91   //
92   // Input slot #0 works with a TChain
93   DefineInput(0, TChain::Class());
94   // Output slot #0 writes into a TList
95   DefineOutput(0, TObjArray::Class());
96   //
97   //
98 }
99
100 AliMCTrackingTestTask::~AliMCTrackingTestTask(){
101   //
102   //
103   //
104   if (fDebugLevel>0)  printf("AliMCTrackingTestTask::~AliMCTrackingTestTask\n");
105   if (fDebugStreamer) delete fDebugStreamer;
106   fDebugStreamer=0;
107 }
108
109
110 //________________________________________________________________________
111 void AliMCTrackingTestTask::ConnectInputData(Option_t *) 
112 {
113   //
114   // Connect the input data
115   //
116   if(fDebugLevel>3)
117     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
118
119   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
120   if (!tree) {
121     //Printf("ERROR: Could not read chain from input slot 0");
122   }
123   else {
124     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
125     if (!esdH) {
126       //Printf("ERROR: Could not get ESDInputHandler");
127     }
128     else {
129       fESD = esdH->GetEvent();
130       //Printf("*** CONNECTED NEW EVENT ****");
131     }  
132   }
133   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());  
134   mcinfo->SetReadTR(kTRUE);
135   
136   fMCinfo = mcinfo->MCEvent();
137
138
139 }
140
141
142
143
144
145
146 //________________________________________________________________________
147 void AliMCTrackingTestTask::CreateOutputObjects() 
148 {
149   //
150   // Connect the output objects
151   //
152   if(fDebugLevel>3)
153     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
154   
155 }
156
157
158 //________________________________________________________________________
159 void AliMCTrackingTestTask::Exec(Option_t *) {
160   //
161   // Execute analysis for current event 
162   //
163
164   if(fDebugLevel>3)
165     cout << "AliMCTrackingTestTask::Exec()" << endl;
166     
167
168   // If MC has been connected   
169
170   if (!fMCinfo){
171     cout << "Not MC info\n" << endl;
172   }else{
173     ProcessMCInfo();
174     //mcinfo->Print();
175     //DumpInfo();
176   }
177   //
178 }      
179
180
181
182
183 //________________________________________________________________________
184 void AliMCTrackingTestTask::Terminate(Option_t *) {
185     //
186     // Terminate loop
187     //
188   if(fDebugLevel>3)
189     printf("AliMCTrackingTestTask: Terminate() \n");  
190   //
191   if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
192   if (fDebugStreamer) delete fDebugStreamer;
193   fDebugStreamer = 0;
194   return;
195 }
196
197
198
199 TTreeSRedirector *AliMCTrackingTestTask::GetDebugStreamer(){
200   //
201   // Get Debug streamer
202   // In case debug streamer not yet initialized and StreamLevel>0 create new one
203   //
204   if (fStreamLevel==0) return 0;
205   if (fDebugStreamer) return fDebugStreamer;
206   TString dsName;
207   dsName=GetName();
208   dsName+="Debug.root";
209   dsName.ReplaceAll(" ","");
210   fDebugStreamer = new TTreeSRedirector(dsName.Data());
211   return fDebugStreamer;
212 }
213
214
215
216
217
218
219 AliExternalTrackParam * AliMCTrackingTestTask::MakeTrack(const AliTrackReference* ref, TParticle*part)
220 {
221   //
222   // Make track out of the track ref
223   // part - TParticle used to determine chargr
224   // the covariance matrix - equal 0 - starting from ideal MC position
225   Double_t xyz[3]={ref->X(),ref->Y(),ref->Z()};
226   Double_t pxyz[3]={ref->Px(),ref->Py(),ref->Pz()};
227   Double_t cv[21];
228   for (Int_t i=0; i<21;i++) cv[i]=0;
229   if (!part->GetPDG()) return 0;
230   AliExternalTrackParam * param = new AliExternalTrackParam(xyz,pxyz,cv,TMath::Nint(part->GetPDG()->Charge()/3.));
231   return param;
232 }
233
234 Bool_t  AliMCTrackingTestTask::PropagateToPoint(AliExternalTrackParam *param, Double_t *xyz, Double_t mass, Float_t step){
235   // 
236   // Propagate track to point xyz using 
237   // AliTracker::PropagateTo functionality
238   //
239   //  param - track parameters
240   //  xyz   - position to propagate
241   //  mass  - particle mass
242   //  step  - step to be used
243   Double_t radius=TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
244   AliTracker::PropagateTrackTo(param, radius+step, mass, step, kTRUE,0.99);
245   AliTracker::PropagateTrackTo(param, radius+0.5, mass, step*0.1, kTRUE,0.99);
246   Double_t sxyz[3]={0,0,0};
247   AliESDVertex vertex(xyz,sxyz);
248   Bool_t isOK = param->PropagateToDCA(&vertex,AliTracker::GetBz(),10);
249   return isOK;
250 }
251
252
253 void  AliMCTrackingTestTask::ProcessMCInfo(){
254   //
255   //
256   //
257    //
258   TParticle * particle= new TParticle;
259   TClonesArray * trefs = new TClonesArray("AliTrackReference");
260   const Double_t kPcut=0.1;
261   //
262   //
263   // Process tracks
264   //
265   Int_t npart = fMCinfo->GetNumberOfTracks();
266   if (npart==0) return;
267   Double_t vertex[4]={0,0,0,0};
268   fMCinfo->GetParticleAndTR(0, particle, trefs);
269   if (particle){
270     vertex[0]=particle->Vx();
271     vertex[1]=particle->Vy();
272     vertex[2]=particle->Vz();
273     vertex[3]=particle->R();
274   }
275   //
276   //
277
278   for (Int_t ipart=0;ipart<npart;ipart++){
279     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
280     if (status<0) continue;
281     if (!particle) continue;
282     if (!trefs) continue;
283     Int_t nref = trefs->GetEntries();
284     if (nref<5) continue;
285     AliTrackReference * tpcIn=0;
286     AliTrackReference * tpcOut=0;
287     AliTrackReference * trdIn=0;
288     AliTrackReference * trdOut=0;
289     AliTrackReference * itsIn=0;
290     AliTrackReference * itsOut=0;
291     Double_t rmax=0;
292     Double_t rmin=1000;
293     for (Int_t iref=0;iref<nref;iref++){
294       AliTrackReference * ref = (AliTrackReference*)trefs->At(iref);
295       if (!ref) continue;
296       
297       Float_t dir = ref->X()*ref->Px()+ref->Y()*ref->Py();
298       
299       if (dir<0) break; // oposite direction - looping track - return back
300       if (ref->P()<kPcut) continue;
301       if (ref->R()<rmax) break;
302       //if (ref->R()<rmin)  break; 
303       //TPC
304       if (ref->DetectorId()==AliTrackReference::kTPC){
305         if (!tpcIn) {
306           tpcIn  = ref;
307         }else{
308           if (ref->R()>tpcIn->R()) tpcOut = ref;
309         }       
310       }
311       //ITS
312       if (ref->DetectorId()==AliTrackReference::kITS){
313         if (!itsIn) {
314           itsIn  = ref;
315         }else{
316           if (ref->R()>itsIn->R()) itsOut = ref;
317         }       
318       }
319       //TRD
320       if (ref->DetectorId()==AliTrackReference::kTRD){
321         if (!trdIn) {
322           trdIn  = ref;
323         }else{
324           if (ref->R()>trdIn->R()) trdOut = ref;
325         }       
326       }      
327       if (ref->R()<rmin) rmin=ref->R();
328       if (ref->R()>rmax) rmax=ref->R();
329     }
330     if (tpcIn && tpcOut) {
331       ProcessRefTracker(tpcIn,tpcOut,particle,1);
332       ProcessRefTracker(tpcIn,tpcOut,particle,3);
333     }
334     if (itsIn && itsOut) ProcessRefTracker(itsIn,itsOut,particle,0);
335     if (trdIn && trdOut) ProcessRefTracker(trdIn,trdOut,particle,2);
336   }
337 }
338
339
340
341 void AliMCTrackingTestTask::ProcessRefTracker(AliTrackReference* refIn,  AliTrackReference* refOut, TParticle*part,Int_t type){
342   //
343   // Test propagation from In to out
344   //
345   AliExternalTrackParam *param = 0;
346   AliExternalTrackParam *paramMC = 0;
347   Double_t xyzIn[3]={refIn->X(),refIn->Y(), refIn->Z()};
348   Double_t mass = part->GetMass();
349   Double_t step=1;
350   //
351   param=MakeTrack(refOut,part);
352   paramMC=MakeTrack(refOut,part);
353   if (!param) return;
354   if (type<3) PropagateToPoint(param,xyzIn, mass, step);
355   if (type==3) {
356     AliTPCseed seed;
357     seed.Set(param->GetX(),param->GetAlpha(),param->GetParameter(),param->GetCovariance());
358     Float_t alpha= TMath::ATan2(refIn->Y(),refIn->X());
359     seed.Rotate(alpha-seed.GetAlpha());
360     seed.SetMass(mass);
361     for (Float_t xlayer= seed.GetX(); xlayer>refIn->R(); xlayer-=step){
362       seed.PropagateTo(xlayer);
363     }
364     seed.PropagateTo(refIn->R());
365     param->Set(seed.GetX(),seed.GetAlpha(),seed.GetParameter(),seed.GetCovariance());
366   }
367   TTreeSRedirector *pcstream = GetDebugStreamer();
368   TVectorD gpos(3);
369   TVectorD gmom(3);
370   param->GetXYZ(gpos.GetMatrixArray());
371   param->GetPxPyPz(gmom.GetMatrixArray());
372   if (pcstream){
373     (*pcstream)<<"MC"<<
374       "type="<<type<<
375       "step="<<step<<
376       "refIn.="<<refIn<<
377       "refOut.="<<refOut<<
378       "p.="<<part<<
379       "par.="<<param<<   
380       "parMC.="<<paramMC<<   
381       "gpos.="<<&gpos<<
382       "gmom.="<<&gmom<<
383       "\n";
384   }
385 }
386
387
388
389 void AliMCTrackingTestTask::FinishTaskOutput()
390 {
391   //
392   // According description in AliAnalisysTask this method is call
393   // on the slaves before sending data
394   //
395   Terminate("slave");
396   gSystem->Exec("pwd");
397   RegisterDebugOutput();
398
399 }
400
401
402 void AliMCTrackingTestTask::RegisterDebugOutput(){
403   //
404   //
405   //
406   //
407   // store  - copy debug output to the destination position
408   // currently ONLY for local copy
409   TString dsName;
410   dsName=GetName();
411   dsName+="Debug.root";
412   dsName.ReplaceAll(" ","");
413   TString dsName2=fDebugOutputPath.Data();
414   gSystem->MakeDirectory(dsName2.Data());
415   dsName2+=gSystem->HostName();
416   gSystem->MakeDirectory(dsName2.Data());
417   dsName2+="/";
418   dsName2+=gSystem->BaseName(gSystem->pwd());
419   dsName2+="/";
420   gSystem->MakeDirectory(dsName2.Data());
421   dsName2+=dsName;
422   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
423   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
424   TFile::Cp(dsName.Data(),dsName2.Data());
425 }
426
427
428 /*
429   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
430   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
431   AliXRDPROOFtoolkit tool; 
432   TChain * chain = tool.MakeChain("mctracking.txt","MC",0,100);
433   chain->Lookup();
434   //
435   //
436   chain->SetAlias("pdg","(p.fPdgCode)");
437   chain->SetAlias("dPRec","(refOut.P()-par.P())/refIn.P()");
438   chain->SetAlias("dPMC","(refOut.P()-refIn.P())/refIn->P()");
439   chain->SetAlias("dPtRec","(refOut.Pt()-par.Pt())/refIn.Pt()");
440   chain->SetAlias("dPtMC","(refOut.Pt()-refIn.Pt())/refIn->Pt()");
441
442
443   // ITS
444   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==0&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
445   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
446   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
447   gPad->SaveAs("picLoss/dPcorr_ITS_step1.gif");
448   gPad->SaveAs("picLoss/dPcorr_ITS_step1.eps");
449   // TPC
450   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==1&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
451   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
452   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
453   gPad->SaveAs("picLoss/dPcorr_TPC_step1.gif");
454   gPad->SaveAs("picLoss/dPcorr_TPC_step1.eps");
455   //
456    // TPC
457   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
458   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
459   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
460   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.gif");
461   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1.eps");
462
463
464   // TRD
465   chain->Draw("-sqrt(dPRec):-sqrt(dPMC)","abs(pdg)!=11&&type==2&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
466   htemp->SetYTitle("#sqrt{#DeltaP_{rec}/P}");
467   htemp->SetXTitle("#sqrt{#DeltaP_{mc}/P}");
468   gPad->SaveAs("picLoss/dPcorr_TRD_step1.gif");
469   gPad->SaveAs("picLoss/dPcorr_TRD_step1.eps");
470
471   //
472   //
473   //
474   chain->Draw("(par.Pt()-refIn.Pt())/refIn.Pt()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
475   his->SetXTitle("(P_{trec}-P_{tmc})/P_{tmc}");
476   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.eps");
477   gPad->SaveAs("picLoss/dPtcorr_TPCseed_step1_1D.gif");
478
479   chain->Draw("(par.P()-refIn.P())/refIn.P()>>his(100,-0.02,0.02)","abs(pdg)!=11&&type==3&&p.Pt()<0.5&&abs(p.R())<1&&abs(refOut.fZ)<220");
480   his->SetXTitle("(P_{rec}-P_{mc})/P_{mc}");
481   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.eps");
482   gPad->SaveAs("picLoss/dPcorr_TPCseed_step1_1D.gif");
483 */