Update of the HFE package
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskDCA.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15 //
16 // The analysis task:
17 // impact parameter resolution and pull study
18 // for tracks which survivied the particle cuts
19 // 
20 // 
21 // Authors:
22 //  Hongyan Yang <hongyan@physi.uni-heidelberg.de>
23 //  Carlo Bombonati <carlo.bombonati@cern.ch>
24 //
25 #include <Riostream.h>
26 #include <TChain.h>
27 #include <TFile.h>
28 #include <TH1F.h>
29 #include <TH1I.h>
30 #include <TList.h>
31 #include <TMath.h>
32 #include <TObjArray.h>
33 #include <TParticle.h>
34 #include <TString.h>
35
36 #include <TCanvas.h>
37
38 #include "AliAnalysisManager.h"
39
40 #include "AliCFManager.h"
41
42 #include "AliESDInputHandler.h"
43 #include "AliESDtrack.h"
44 #include "AliVertexerTracks.h"
45 #include "AliESDVertex.h"
46
47 #include "AliMCEventHandler.h"
48 #include "AliMCEvent.h"
49 #include "AliMCParticle.h"
50
51 #include "AliESDpid.h"
52 #include "AliHFEpid.h"
53 #include "AliHFEcuts.h"
54 #include "AliHFEdca.h"
55 #include "AliHFEtools.h"
56
57 #include "AliAnalysisTaskDCA.h"
58
59
60 //____________________________________________________________
61 AliAnalysisTaskDCA::AliAnalysisTaskDCA():
62   AliAnalysisTaskSE("Impact Parameter Resolution and Pull Analysis")
63   , fPlugins(0)
64   , fCuts(0x0)
65   , fDefaultPID(0x0)
66   , fHFEpid(0x0)
67   , fPIDdetectors("")
68   , fPIDstrategy(0)
69   , fCFM(0x0)
70   , fDCA(0x0)
71   , fNclustersITS(0x0)
72   , fMinNprimVtxContrbutor(0x0)
73   , fNEvents(0x0)
74   , fResidualList(0x0)
75   , fPullList(0x0)
76   , fDcaList(0x0)
77   , fKfDcaList(0x0)
78   , fMcVertexList(0x0)
79   , fDataDcaList(0x0)
80   , fDataVertexList(0x0)
81   , fDataPullList(0x0)
82   , fMcPidList(0x0) 
83   , fDataPidList(0x0)
84   , fHfeDcaList(0x0)
85   , fHfeDataDcaList(0x0)
86   , fOutput(0x0)
87 {
88   //
89   // Dummy constructor
90   //
91 }
92
93 //____________________________________________________________
94 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const char * name):
95   AliAnalysisTaskSE(name)
96   , fPlugins(0)
97   , fCuts(0x0)
98   , fDefaultPID(0x0)
99   , fHFEpid(0x0)
100   , fPIDdetectors("")
101   , fPIDstrategy(0)
102   , fCFM(0x0)
103   , fDCA(0x0)
104   , fNclustersITS(0x0)
105   , fMinNprimVtxContrbutor(0x0)
106   , fNEvents(0x0)
107   , fResidualList(0x0)
108   , fPullList(0x0)
109   , fDcaList(0x0) 
110   , fKfDcaList(0x0)
111   , fMcVertexList(0x0)
112   , fDataDcaList(0x0)
113   , fDataVertexList(0x0)
114   , fDataPullList(0x0)
115   , fMcPidList(0x0)
116   , fDataPidList(0x0)
117   , fHfeDcaList(0x0)
118   , fHfeDataDcaList(0x0)
119   , fOutput(0x0)
120 {
121   //
122   // Default constructor
123   //
124   DefineInput(0, TChain::Class());
125   DefineOutput(1, TH1I::Class());
126   DefineOutput(2, TList::Class());
127
128
129   //-CUTS SETTING-//
130   Int_t nMinTPCcluster = 100;
131   Float_t maxDcaXY = 0.5;
132   Float_t maxDcaZ = 1.0;
133   //--------------//
134   AliHFEcuts *hfecuts = new AliHFEcuts;
135   hfecuts->CreateStandardCuts();
136   hfecuts->SetMinNClustersTPC(nMinTPCcluster);
137   hfecuts->SetCutITSpixel(AliHFEextraCuts::kFirst);
138   hfecuts->SetCheckITSLayerStatus(kFALSE);
139   hfecuts->SetMaxImpactParam(maxDcaXY, maxDcaZ);
140   SetHFECuts(hfecuts);
141  
142   fDefaultPID = new AliESDpid();
143   fHFEpid = new AliHFEpid("PIDforDCAanalysis");
144
145 }
146
147 //____________________________________________________________
148 AliAnalysisTaskDCA::AliAnalysisTaskDCA(const AliAnalysisTaskDCA &ref):
149   AliAnalysisTaskSE(ref)
150   , fPlugins(ref.fPlugins)
151   , fCuts(ref.fCuts)  
152   , fDefaultPID(ref.fDefaultPID)
153   , fHFEpid(ref.fHFEpid)
154   , fPIDdetectors(ref.fPIDdetectors)
155   , fPIDstrategy(ref.fPIDstrategy)
156   , fCFM(ref.fCFM)
157   , fDCA(ref.fDCA)
158   , fNclustersITS(ref.fNclustersITS)
159   , fMinNprimVtxContrbutor(ref.fMinNprimVtxContrbutor)
160   , fNEvents(ref.fNEvents)
161   , fResidualList(ref.fResidualList)
162   , fPullList(ref.fPullList)
163   , fDcaList(ref.fDcaList)
164   , fKfDcaList(ref.fKfDcaList)
165   , fMcVertexList(ref.fMcVertexList)
166   , fDataDcaList(ref.fDataDcaList)
167   , fDataVertexList(ref.fDataVertexList)
168   , fDataPullList(ref.fDataPullList)
169   , fMcPidList(ref.fMcPidList)
170   , fDataPidList(ref.fDataPidList)
171   , fHfeDcaList(ref.fHfeDcaList)
172   , fHfeDataDcaList(ref.fHfeDataDcaList)
173   , fOutput(ref.fOutput)
174 {
175   //
176   // Copy Constructor
177   //
178   AliInfo("Copy Constructor");
179   ref.Copy(*this);
180 }
181
182 //____________________________________________________________
183 AliAnalysisTaskDCA &AliAnalysisTaskDCA::operator=(const AliAnalysisTaskDCA &ref){
184   //
185   // Assignment operator
186   //
187   if(this == &ref) return *this;
188   AliAnalysisTaskSE::operator=(ref);
189   fPlugins = ref.fPlugins;
190   fCuts = ref.fCuts;
191   fDefaultPID = ref.fDefaultPID;
192   fHFEpid = ref.fHFEpid;
193   fPIDdetectors = ref.fPIDdetectors;
194   fPIDstrategy = ref.fPIDstrategy;
195   fCFM = ref.fCFM;
196   fDCA = ref.fDCA;
197   fNclustersITS = ref.fNclustersITS;
198   fMinNprimVtxContrbutor = ref.fMinNprimVtxContrbutor;
199   fNEvents = ref.fNEvents;
200   fResidualList = ref.fResidualList;
201   fPullList = ref.fPullList;
202   fDcaList = ref.fDcaList;
203   fKfDcaList = ref.fKfDcaList;
204   fMcVertexList = ref.fMcVertexList;
205   fDataDcaList = ref.fDataDcaList;
206   fDataVertexList = ref.fDataVertexList;
207   fDataPullList = ref.fDataPullList;
208   fMcPidList = ref.fMcPidList;
209   fDataPidList = ref.fDataPidList;
210   fHfeDcaList = ref.fHfeDcaList;    
211   fHfeDataDcaList = ref.fHfeDataDcaList;
212   fOutput = ref.fOutput;
213
214   return *this;
215 }
216
217 //____________________________________________________________
218 AliAnalysisTaskDCA::~AliAnalysisTaskDCA(){
219   //
220   // Destructor
221   //
222
223   if(fDefaultPID) delete fDefaultPID;
224   if(fHFEpid) delete fHFEpid;
225   if(fCFM) delete fCFM;
226   if(fDCA) delete fDCA;  
227   if(fNEvents) delete fNEvents;
228   if(fResidualList){ 
229     fResidualList->Clear();
230     delete fResidualList;
231   }
232
233   if(fPullList){ 
234     fPullList->Clear();
235     delete fPullList;
236   }
237   
238   if(fDcaList){ 
239     fDcaList->Clear();
240     delete fDcaList;
241   }
242   if(fKfDcaList){ 
243     fKfDcaList->Clear();
244     delete fKfDcaList;
245   }
246
247   if(fMcVertexList){
248     fMcVertexList->Clear();
249     delete   fMcVertexList;
250   }
251
252   if(fDataDcaList){ 
253     fDataDcaList->Clear();
254     delete fDataDcaList;
255   }
256
257   if(fDataVertexList){
258     fDataVertexList->Clear();
259     delete   fDataVertexList;
260   }
261   if(fDataPullList){ 
262     fDataPullList->Clear();
263     delete fDataPullList;
264   }
265
266   if(fMcPidList){
267     fMcPidList -> Clear();
268     delete fMcPidList;
269   }
270   if(fDataPidList){
271     fDataPidList -> Clear();
272     delete fDataPidList;
273   }
274
275   if(fHfeDcaList) {
276     fHfeDcaList->Clear();
277     delete fHfeDcaList;
278   } 
279   
280   if(fHfeDataDcaList) {
281     fHfeDataDcaList->Clear();
282     delete fHfeDataDcaList;
283   } 
284   
285   if(fOutput){ 
286     fOutput->Clear();
287     delete fOutput;
288   }
289   
290 }
291
292 //____________________________________________________________
293 void AliAnalysisTaskDCA::UserCreateOutputObjects(){
294   // create output objects
295   // fNEvents
296   // residual and pull
297   //printf("\n=====UserCreateOutputObjects=====\n");
298   
299   // Automatic determination of the analysis mode
300   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
301   if(!inputHandler){
302     AliError("NoEvent Handler available");
303     return;
304   }
305
306   if(!TString(inputHandler->IsA()->GetName()).CompareTo("AliAODInputHandler")){
307     SetAODAnalysis();
308   } else {
309     SetESDAnalysis();
310     if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())
311       SetHasMCData();
312   }
313   
314
315   fNEvents = new TH1I("nEvents", "Number of Events in the Analysis", 5, -0.5, 4.5); // Number of Events neccessary for the analysis and not a QA histogram
316   if(!fOutput) fOutput = new TList;
317   // Initialize correction Framework and Cuts
318   fCFM = new AliCFManager;
319   MakeParticleContainer();
320   // Temporary fix: Initialize particle cuts with 0x0
321   for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++)
322     fCFM->SetParticleCutsList(istep, 0x0);
323   if(!fCuts){
324     AliWarning("Cuts not available. Default cuts will be used");
325     fCuts = new AliHFEcuts;
326     fCuts->CreateStandardCuts();
327   }
328   
329   fCuts->Initialize(fCFM);
330   
331   if(!fHFEpid) AliWarning("Hello, fHFEpid is not available");
332   cout<<"  Hello this is a cout "<<endl<<endl;
333
334   if(GetPlugin(kHFEpid)) {   
335     fHFEpid->SetHasMCData(HasMCData());
336     fHFEpid->AddDetector("TOF", 0);
337     fHFEpid->AddDetector("TPC", 1);
338     cout<<endl<<" ---> TPC and TOF added to the PID"<<endl;
339                 fHFEpid->ConfigureTOF();
340     fHFEpid->ConfigureTPCdefaultCut();
341     fHFEpid->InitializePID();
342   }
343
344   // dca study----------------------------------
345
346   
347   if(!fDCA) fDCA = new AliHFEdca;
348   if(!fResidualList) fResidualList = new TList();
349   if(!fPullList) fPullList = new TList();
350   if(!fDcaList) fDcaList = new TList();
351   if(!fKfDcaList) fKfDcaList = new TList();
352   if(!fMcVertexList) fMcVertexList = new TList();
353   if(!fDataDcaList) fDataDcaList = new TList();
354   if(!fDataVertexList) fDataVertexList = new TList();
355   if(!fDataPullList) fDataPullList = new TList();
356   if(!fMcPidList) fMcPidList = new TList();
357   if(!fDataPidList) fDataPidList = new TList();  
358   
359   if(!fHfeDcaList) fHfeDcaList = new TList();
360   if(!fHfeDataDcaList) fHfeDataDcaList = new TList();
361
362   if(HasMCData()) {    
363     if(GetPlugin(kImpactPar) ) {
364       fDCA->CreateHistogramsResidual(fResidualList);
365       fDCA->CreateHistogramsPull(fPullList);
366       fDCA->CreateHistogramsDca(fDcaList);
367       fOutput->AddAt(fResidualList,0);
368       fOutput->AddAt(fPullList,1);
369       fOutput->AddAt(fDcaList,2);
370     } 
371     if(GetPlugin(kKFdca)){
372       fDCA->CreateHistogramsKfDca(fKfDcaList);
373       fOutput->AddAt(fDcaList,3);
374     }
375     if(GetPlugin(kPrimVtx)){//<---
376       fDCA->CreateHistogramsVertex(fMcVertexList);
377       fOutput->AddAt(fMcVertexList,4);
378     }
379     if(GetPlugin(kCombinedPid)){//<---
380       fDCA->CreateHistogramsPid(fMcPidList);
381       fOutput->AddAt(fMcPidList, 5);
382     }
383     if(GetPlugin(kHFEpid)){//<---
384       fDCA->CreateHistogramsHfeDca(fHfeDcaList);
385       fOutput->AddAt(fHfeDcaList, 6);
386     }
387   } // mc case
388
389   if(!HasMCData())  { 
390     
391     if(GetPlugin(kPrimVtx)){
392       fDCA->CreateHistogramsDataVertex(fDataVertexList);  
393       fOutput->AddAt(fDataVertexList,0);
394     }    
395
396     if(GetPlugin(kCombinedPid)){
397       fDCA->CreateHistogramsDataDca(fDataDcaList);  
398       fDCA->CreateHistogramsDataPull(fDataPullList);  
399       fDCA->CreateHistogramsDataPid(fDataPidList);
400       fOutput->AddAt(fDataDcaList,1);
401       fOutput->AddAt(fDataPullList,2);
402       fOutput->AddAt(fDataPidList, 3);
403     }
404     if(GetPlugin(kHFEpid)){
405       fDCA->CreateHistogramsHfeDataDca(fHfeDataDcaList);
406       fOutput->AddAt(fHfeDataDcaList, 4);
407     }
408     
409
410
411   }  // data case
412   
413 }
414
415 //____________________________________________________________
416 void AliAnalysisTaskDCA::UserExec(Option_t *){
417   //
418   // Run the analysis
419   // 
420   //printf("\n=====UserExec=====\n");
421   if(HasMCData()) printf("WITH MC!\n");
422
423   AliDebug(3, "Processing ESD events");
424
425   if(!fInputEvent){
426     AliError("Reconstructed Event not available");
427     return;
428   }
429   if(HasMCData()){
430     AliDebug(4, Form("MC Event: %p", fMCEvent));
431     if(!fMCEvent){
432       AliError("No MC Event, but MC Data required");
433       return;
434     }
435   }
436
437   if(!fCuts){
438     AliError("HFE cuts not available");
439     return;
440   }
441  
442
443   // protection
444   if(IsESDanalysis() && HasMCData()){
445     // Protect against missing MC trees
446     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
447     if(!mcH){
448       AliError("No MC Event Handler available");
449       return;
450     }
451     if(!mcH->InitOk()) return;
452     if(!mcH->TreeK()) return;
453     if(!mcH->TreeTR()) return;
454   }
455
456   if(!IsAODanalysis()) {
457     AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
458     if(!inH){
459       AliError("No ESD Event Handler available");
460       return;
461     }
462     AliESDpid *workingPID = inH->GetESDpid();
463     if(workingPID){
464       AliDebug(1, "Using ESD PID from the input handler");
465       fHFEpid->SetESDpid(workingPID);
466     } else {
467       AliDebug(1, "Using default ESD PID");
468       fHFEpid->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
469     }
470     ProcessDcaAnalysis();
471   }
472
473   
474   PostData(1, fNEvents);
475   PostData(2, fOutput);
476 }
477 //____________________________________________________________
478 void AliAnalysisTaskDCA::ProcessDcaAnalysis(){
479
480   //printf("\n=====ProcessDcaAnalysis=====\n");
481
482   //
483   // Loop ESD
484   //
485   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
486   if(!fESD){
487     AliError("ESD Event required for ESD Analysis");
488       return;
489   }
490
491   AliMCEvent *fMC = 0x0;
492   if(HasMCData()){
493     fMC = dynamic_cast<AliMCEvent*>(fMCEvent);
494     if(!fMC){
495       AliError("MC Event required for Analysis");
496       return;
497     }
498   }
499
500   fNEvents->Fill(1);  // original event number before cut
501   fDCA->ApplyExtraCuts(fESD,fMinNprimVtxContrbutor);  // cut on primVtx contributors
502   fNEvents->Fill(3);  // events number after cut
503   fCFM->SetRecEventInfo(fESD);
504  
505   // event cut level
506   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
507
508   AliESDtrack *track = 0x0;  
509   AliMCParticle *mctrack = 0x0;
510   AliESDVertex *vtxESDSkip = 0x0;
511   AliHFEpidObject hfetrack;
512
513   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
514     
515     track = fESD->GetTrack(itrack);
516     if(HasMCData()) mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel())));
517
518     // RecPrim: primary cuts
519     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
520     // RecKine: ITSTPC cuts  
521     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
522     // HFEcuts: ITS layers cuts
523     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
524     
525     if(track->GetITSclusters(0)<=fNclustersITS) continue;  // require number of ITS clusters
526     
527     // track accepted, do PID
528     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
529     hfetrack.SetRecTrack(track);
530     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
531
532     //printf("Track %d passed all the cuts!\n",itrack);
533
534     if(HasMCData()){
535       if(GetPlugin(kPrimVtx))
536         fDCA->FillHistogramsVtx(fESD, fMC);
537       if(GetPlugin(kImpactPar)) 
538         fDCA->FillHistogramsDca(fESD, track, fMC);
539       if(GetPlugin(kKFdca)) 
540         fDCA->FillHistogramsKfDca(fESD, track, fMC);
541       if(GetPlugin(kCombinedPid)) 
542         fDCA->FillHistogramsPid(track, fMC);
543       if(GetPlugin(kHFEpid)) { // data-like
544         if(fHFEpid->IsSelected(&hfetrack)){ 
545           
546           //      printf("Found an electron in p+p collision! from HFE pid \n");
547           if(!vtxESDSkip){
548             // method from Andrea D 28.05.2010
549             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
550             vertexer->SetITSMode();
551             vertexer->SetMinClusters(fNclustersITS);
552             Int_t skipped[2];
553             skipped[0] = (Int_t)track->GetID();
554             vertexer->SetSkipTracks(1,skipped);
555             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
556             delete vertexer; vertexer = NULL;
557             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
558           }
559           //printf("\n[ABOUT TO FILL HFE DCA: MC!]\n");
560           fDCA->FillHistogramsHfeDataDca(fESD, track, vtxESDSkip); 
561         }
562       } // plugin for hfepid 
563     }  // MC
564
565     if(!HasMCData()){
566       if(GetPlugin(kPrimVtx))
567         fDCA->FillHistogramsDataVtx(fESD);
568       if(GetPlugin(kCombinedPid)) {
569
570         // method from Andrea D 28.05.2010
571         AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
572         vertexer->SetITSMode();
573         vertexer->SetMinClusters(fNclustersITS);
574         Int_t skipped[2];
575         skipped[0] = (Int_t)track->GetID();
576         vertexer->SetSkipTracks(1,skipped);
577         vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
578         delete vertexer; vertexer = NULL;
579         if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
580
581         fDCA->FillHistogramsDataDca(fESD, track, vtxESDSkip);
582         fDCA->FillHistogramsDataPid(track);
583       }
584       if(GetPlugin(kHFEpid)) {
585         if(fHFEpid->IsSelected(&hfetrack)) {
586           //      printf("Found an electron in p+p collision! from HFE pid \n");
587           if(!vtxESDSkip){
588             // method from Andrea D 28.05.2010
589             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
590             vertexer->SetITSMode();
591             vertexer->SetMinClusters(fNclustersITS);
592             Int_t skipped[2];
593             skipped[0] = (Int_t)track->GetID();
594             vertexer->SetSkipTracks(1,skipped);
595             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
596             delete vertexer; vertexer = NULL;
597             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
598           }
599         printf("\n[ABOUT TO FILL HFE DCA: DATA!]\n");
600           fDCA->FillHistogramsHfeDataDca(fESD, track,vtxESDSkip);    
601         } 
602       } // plugin for hfepid
603     }  // data case
604
605   } // track loop
606
607 }
608
609
610 //____________________________________________________________
611 void AliAnalysisTaskDCA::Terminate(Option_t *){
612   //
613   // Terminate not implemented at the moment
614   //  
615   //printf("\n=====Terminate=====\n");
616   
617   if(GetPlugin(kPostProcess)){
618     fOutput = dynamic_cast<TList *>(GetOutputData(1));
619     if(!fOutput){
620       AliError("Results not available");
621       return;
622     }
623     PostProcess();
624   }
625   
626 }
627
628
629 //____________________________________________________________
630 void AliAnalysisTaskDCA::Load(TString filename){
631
632   //printf("\n=====Load=====\n");
633
634   // no need for postprocessing for the moment
635   TFile *input = TFile::Open(filename.Data());
636   if(!input || input->IsZombie()){
637     AliError("Cannot read file");
638     return;
639   }
640
641   input->Close();
642   delete input;
643   
644   
645 }
646
647 //____________________________________________________________
648 void AliAnalysisTaskDCA::PostProcess(){
649   // do post processing
650   // should do fitting here for dca resolution
651   // moved to an external macro to do the job
652   
653   //printf("\n=====PostProcess=====\n");
654   Load("HFEdca.root");
655   TCanvas *c1 = new TCanvas("c1", "number of analyzed events", 300, 400);
656   fNEvents->Draw();
657   c1->SaveAs("temp.png");
658  
659 }
660
661
662
663
664 //____________________________________________________________
665 void AliAnalysisTaskDCA::PrintStatus() const {
666   
667   //
668   // Print Analysis status
669   //
670   printf("\n\tAnalysis Settings\n\t========================================\n");
671   printf("\t Running on %s\n", !HasMCData()?"p+p collision data":"MC sample");
672   printf("\t Cuts: %s\n", (fCuts != NULL) ? "YES" : "NO");
673   printf("\t Impact parameter analysis is %s\n", GetPlugin(kImpactPar)?"ON":"OFF");
674   printf("\t Using AliKFParticle for analysis? %s\n", GetPlugin(kKFdca)?"ON":"OFF");
675   printf("\t Primary vertex analysis is %s\n", GetPlugin(kPrimVtx)?"ON":"OFF");
676   printf("\t Combined pid analysis is %s\n", GetPlugin(kCombinedPid)?"ON":"OFF");
677   printf("\t HFE pid analysis is %s\n", GetPlugin(kHFEpid)?"ON":"OFF");
678   printf("\t Post process analysis is %s\n", GetPlugin(kPostProcess)?"ON":"OFF");
679   printf("\t ");
680   printf("\n");
681 }
682
683 //__________________________________________                                                  
684 void AliAnalysisTaskDCA::SwitchOnPlugin(Int_t plug){
685   //                                            
686   // Switch on Plugin          
687   // Available:                                  
688   //  - analyze impact parameter
689   //  - Post Processing  
690
691   AliDebug(2,Form("SwitchOnPlugin %d",plug));  
692
693   switch(plug){
694   case kPostProcess: 
695     SETBIT(fPlugins, plug); 
696     break;
697   case kImpactPar: 
698     SETBIT(fPlugins, plug); 
699     break;
700   case kPrimVtx: 
701     SETBIT(fPlugins, plug); 
702     break;
703   case kCombinedPid:
704     SETBIT(fPlugins, plug); 
705     break;
706   case kHFEpid:
707     SETBIT(fPlugins, plug); 
708     break;
709   case kKFdca:
710     SETBIT(fPlugins, plug); 
711     break;
712   default: 
713     AliError("Unknown Plugin");
714   };
715 }
716
717
718 //____________________________________________________________
719 void AliAnalysisTaskDCA::MakeParticleContainer(){
720
721   //printf("\n=====MakeParticleContainer=====\n");
722   //
723   // Create the particle container (borrowed from AliAnalysisTaskHFE)
724   //
725   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
726   const Double_t kPtmin = 0.1, kPtmax = 10.;
727   const Double_t kEtamin = -0.9, kEtamax = 0.9;
728   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
729
730   //arrays for the number of bins in each dimension
731   Int_t iBin[kNvar];
732   iBin[0] = 40; //bins in pt
733   iBin[1] =  8; //bins in eta 
734   iBin[2] = 18; // bins in phi
735
736   //arrays for lower bounds :
737   Double_t* binEdges[kNvar];
738   for(Int_t ivar = 0; ivar < kNvar; ivar++)
739     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
740
741   //values for bin lower bounds
742   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);  
743   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
744   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
745
746   //one "container" for MC
747   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
748   const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
749   AliCFContainer* container = new AliCFContainer("container","container for tracks", (kNcutStepsTrack + 1 + 2*(kNcutStepsESDtrack + 1)), kNvar, iBin);
750
751   //setting the bin limits
752   for(Int_t ivar = 0; ivar < kNvar; ivar++)
753     container -> SetBinLimits(ivar, binEdges[ivar]);
754   fCFM->SetParticleContainer(container);
755
756   //create correlation matrix for unfolding
757   Int_t thnDim[2*kNvar];
758   for (int k=0; k<kNvar; k++) {
759     //first half  : reconstructed 
760     //second half : MC
761     thnDim[k]      = iBin[k];
762     thnDim[k+kNvar] = iBin[k];
763   }
764
765
766 }
767
768 //____________________________________________________________
769 void AliAnalysisTaskDCA::AddPIDdetector(TString detector){
770   
771   //
772   // Adding PID detector to the task
773   //
774   //printf("\n=====AddPIDdetector=====\n");
775   
776   if(!fPIDdetectors.Length()) 
777     fPIDdetectors = detector;
778   else
779     fPIDdetectors += ":" + detector;
780 }
781