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