2480a44ed30e3a02d123f8a509f299291cbec18a
[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->ConfigureTPCrejection();
340     fHFEpid->InitializePID();
341   }
342
343   // dca study----------------------------------
344
345   
346   if(!fDCA) fDCA = new AliHFEdca;
347   if(!fResidualList) fResidualList = new TList();
348   if(!fPullList) fPullList = new TList();
349   if(!fDcaList) fDcaList = new TList();
350   if(!fKfDcaList) fKfDcaList = new TList();
351   if(!fMcVertexList) fMcVertexList = new TList();
352   if(!fDataDcaList) fDataDcaList = new TList();
353   if(!fDataVertexList) fDataVertexList = new TList();
354   if(!fDataPullList) fDataPullList = new TList();
355   if(!fMcPidList) fMcPidList = new TList();
356   if(!fDataPidList) fDataPidList = new TList();  
357   
358   if(!fHfeDcaList) fHfeDcaList = new TList();
359   if(!fHfeDataDcaList) fHfeDataDcaList = new TList();
360
361   if(HasMCData()) {    
362     if(GetPlugin(kImpactPar) ) {
363       fDCA->CreateHistogramsResidual(fResidualList);
364       fDCA->CreateHistogramsPull(fPullList);
365       fDCA->CreateHistogramsDca(fDcaList);
366       fOutput->AddAt(fResidualList,0);
367       fOutput->AddAt(fPullList,1);
368       fOutput->AddAt(fDcaList,2);
369     } 
370     if(GetPlugin(kKFdca)){
371       fDCA->CreateHistogramsKfDca(fKfDcaList);
372       fOutput->AddAt(fDcaList,3);
373     }
374     if(GetPlugin(kPrimVtx)){//<---
375       fDCA->CreateHistogramsVertex(fMcVertexList);
376       fOutput->AddAt(fMcVertexList,4);
377     }
378     if(GetPlugin(kCombinedPid)){//<---
379       fDCA->CreateHistogramsPid(fMcPidList);
380       fOutput->AddAt(fMcPidList, 5);
381     }
382     if(GetPlugin(kHFEpid)){//<---
383       fDCA->CreateHistogramsHfeDca(fHfeDcaList);
384       fOutput->AddAt(fHfeDcaList, 6);
385     }
386   } // mc case
387
388   if(!HasMCData())  { 
389     
390     if(GetPlugin(kPrimVtx)){
391       fDCA->CreateHistogramsDataVertex(fDataVertexList);  
392       fOutput->AddAt(fDataVertexList,0);
393     }    
394
395     if(GetPlugin(kCombinedPid)){
396       fDCA->CreateHistogramsDataDca(fDataDcaList);  
397       fDCA->CreateHistogramsDataPull(fDataPullList);  
398       fDCA->CreateHistogramsDataPid(fDataPidList);
399       fOutput->AddAt(fDataDcaList,1);
400       fOutput->AddAt(fDataPullList,2);
401       fOutput->AddAt(fDataPidList, 3);
402     }
403     if(GetPlugin(kHFEpid)){
404       fDCA->CreateHistogramsHfeDataDca(fHfeDataDcaList);
405       fOutput->AddAt(fHfeDataDcaList, 4);
406     }
407     
408
409
410   }  // data case
411   
412 }
413
414 //____________________________________________________________
415 void AliAnalysisTaskDCA::UserExec(Option_t *){
416   //
417   // Run the analysis
418   // 
419   //printf("\n=====UserExec=====\n");
420   if(HasMCData()) printf("WITH MC!\n");
421
422   AliDebug(3, "Processing ESD events");
423
424   if(!fInputEvent){
425     AliError("Reconstructed Event not available");
426     return;
427   }
428   if(HasMCData()){
429     AliDebug(4, Form("MC Event: %p", fMCEvent));
430     if(!fMCEvent){
431       AliError("No MC Event, but MC Data required");
432       return;
433     }
434   }
435
436   if(!fCuts){
437     AliError("HFE cuts not available");
438     return;
439   }
440  
441
442   // protection
443   if(IsESDanalysis() && HasMCData()){
444     // Protect against missing MC trees
445     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
446     if(!mcH){
447       AliError("No MC Event Handler available");
448       return;
449     }
450     if(!mcH->InitOk()) return;
451     if(!mcH->TreeK()) return;
452     if(!mcH->TreeTR()) return;
453   }
454
455   if(!IsAODanalysis()) {
456     AliESDInputHandler *inH = dynamic_cast<AliESDInputHandler *>(fInputHandler);
457     if(!inH){
458       AliError("No ESD Event Handler available");
459       return;
460     }
461     AliESDpid *workingPID = inH->GetESDpid();
462     if(workingPID){
463       AliDebug(1, "Using ESD PID from the input handler");
464       fHFEpid->SetESDpid(workingPID);
465     } else {
466       AliDebug(1, "Using default ESD PID");
467       fHFEpid->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
468     }
469     ProcessDcaAnalysis();
470   }
471
472   
473   PostData(1, fNEvents);
474   PostData(2, fOutput);
475 }
476 //____________________________________________________________
477 void AliAnalysisTaskDCA::ProcessDcaAnalysis(){
478
479   //printf("\n=====ProcessDcaAnalysis=====\n");
480
481   //
482   // Loop ESD
483   //
484   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
485   if(!fESD){
486     AliError("ESD Event required for ESD Analysis");
487       return;
488   }
489
490   AliMCEvent *fMC = 0x0;
491   if(HasMCData()){
492     fMC = dynamic_cast<AliMCEvent*>(fMCEvent);
493     if(!fMC){
494       AliError("MC Event required for Analysis");
495       return;
496     }
497   }
498
499   fNEvents->Fill(1);  // original event number before cut
500   fDCA->ApplyExtraCuts(fESD,fMinNprimVtxContrbutor);  // cut on primVtx contributors
501   fNEvents->Fill(3);  // events number after cut
502   fCFM->SetRecEventInfo(fESD);
503  
504   // event cut level
505   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fESD)) return;
506
507   AliESDtrack *track = 0x0;  
508   AliMCParticle *mctrack = 0x0;
509   AliESDVertex *vtxESDSkip = 0x0;
510   AliHFEpidObject hfetrack;
511
512   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
513     
514     track = fESD->GetTrack(itrack);
515     if(HasMCData()) mctrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(TMath::Abs(track->GetLabel())));
516
517     // RecPrim: primary cuts
518     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecPrim + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
519     // RecKine: ITSTPC cuts  
520     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
521     // HFEcuts: ITS layers cuts
522     if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS + AliHFEcuts::kNcutStepsMCTrack, track)) continue;
523     
524     if(track->GetITSclusters(0)<=fNclustersITS) continue;  // require number of ITS clusters
525     
526     // track accepted, do PID
527     hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis);
528     hfetrack.SetRecTrack(track);
529     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
530
531     //printf("Track %d passed all the cuts!\n",itrack);
532
533     if(HasMCData()){
534       if(GetPlugin(kPrimVtx))
535         fDCA->FillHistogramsVtx(fESD, fMC);
536       if(GetPlugin(kImpactPar)) 
537         fDCA->FillHistogramsDca(fESD, track, fMC);
538       if(GetPlugin(kKFdca)) 
539         fDCA->FillHistogramsKfDca(fESD, track, fMC);
540       if(GetPlugin(kCombinedPid)) 
541         fDCA->FillHistogramsPid(track, fMC);
542       if(GetPlugin(kHFEpid)) { // data-like
543         if(fHFEpid->IsSelected(&hfetrack)){ 
544           
545           //      printf("Found an electron in p+p collision! from HFE pid \n");
546           if(!vtxESDSkip){
547             // method from Andrea D 28.05.2010
548             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
549             vertexer->SetITSMode();
550             vertexer->SetMinClusters(fNclustersITS);
551             Int_t skipped[2];
552             skipped[0] = (Int_t)track->GetID();
553             vertexer->SetSkipTracks(1,skipped);
554             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
555             delete vertexer; vertexer = NULL;
556             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
557           }
558           //printf("\n[ABOUT TO FILL HFE DCA: MC!]\n");
559           fDCA->FillHistogramsHfeDataDca(fESD, track, vtxESDSkip); 
560         }
561       } // plugin for hfepid 
562     }  // MC
563
564     if(!HasMCData()){
565       if(GetPlugin(kPrimVtx))
566         fDCA->FillHistogramsDataVtx(fESD);
567       if(GetPlugin(kCombinedPid)) {
568
569         // method from Andrea D 28.05.2010
570         AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
571         vertexer->SetITSMode();
572         vertexer->SetMinClusters(fNclustersITS);
573         Int_t skipped[2];
574         skipped[0] = (Int_t)track->GetID();
575         vertexer->SetSkipTracks(1,skipped);
576         vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
577         delete vertexer; vertexer = NULL;
578         if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
579
580         fDCA->FillHistogramsDataDca(fESD, track, vtxESDSkip);
581         fDCA->FillHistogramsDataPid(track);
582       }
583       if(GetPlugin(kHFEpid)) {
584         if(fHFEpid->IsSelected(&hfetrack)) {
585           //      printf("Found an electron in p+p collision! from HFE pid \n");
586           if(!vtxESDSkip){
587             // method from Andrea D 28.05.2010
588             AliVertexerTracks *vertexer = new AliVertexerTracks(fESD->GetMagneticField());
589             vertexer->SetITSMode();
590             vertexer->SetMinClusters(fNclustersITS);
591             Int_t skipped[2];
592             skipped[0] = (Int_t)track->GetID();
593             vertexer->SetSkipTracks(1,skipped);
594             vtxESDSkip = (AliESDVertex*)vertexer->FindPrimaryVertex(fESD);
595             delete vertexer; vertexer = NULL;
596             if(vtxESDSkip->GetNContributors()<fMinNprimVtxContrbutor) continue;
597           }
598         printf("\n[ABOUT TO FILL HFE DCA: DATA!]\n");
599           fDCA->FillHistogramsHfeDataDca(fESD, track,vtxESDSkip);    
600         } 
601       } // plugin for hfepid
602     }  // data case
603
604   } // track loop
605
606 }
607
608
609 //____________________________________________________________
610 void AliAnalysisTaskDCA::Terminate(Option_t *){
611   //
612   // Terminate not implemented at the moment
613   //  
614   //printf("\n=====Terminate=====\n");
615   
616   if(GetPlugin(kPostProcess)){
617     fOutput = dynamic_cast<TList *>(GetOutputData(1));
618     if(!fOutput){
619       AliError("Results not available");
620       return;
621     }
622     PostProcess();
623   }
624   
625 }
626
627
628 //____________________________________________________________
629 void AliAnalysisTaskDCA::Load(TString filename){
630
631   //printf("\n=====Load=====\n");
632
633   // no need for postprocessing for the moment
634   TFile *input = TFile::Open(filename.Data());
635   if(!input || input->IsZombie()){
636     AliError("Cannot read file");
637     return;
638   }
639
640   input->Close();
641   delete input;
642   
643   
644 }
645
646 //____________________________________________________________
647 void AliAnalysisTaskDCA::PostProcess(){
648   // do post processing
649   // should do fitting here for dca resolution
650   // moved to an external macro to do the job
651   
652   //printf("\n=====PostProcess=====\n");
653   Load("HFEdca.root");
654   TCanvas *c1 = new TCanvas("c1", "number of analyzed events", 300, 400);
655   fNEvents->Draw();
656   c1->SaveAs("temp.png");
657  
658 }
659
660
661
662
663 //____________________________________________________________
664 void AliAnalysisTaskDCA::PrintStatus() const {
665   
666   //
667   // Print Analysis status
668   //
669   printf("\n\tAnalysis Settings\n\t========================================\n");
670   printf("\t Running on %s\n", !HasMCData()?"p+p collision data":"MC sample");
671   printf("\t Cuts: %s\n", (fCuts != NULL) ? "YES" : "NO");
672   printf("\t Impact parameter analysis is %s\n", GetPlugin(kImpactPar)?"ON":"OFF");
673   printf("\t Using AliKFParticle for analysis? %s\n", GetPlugin(kKFdca)?"ON":"OFF");
674   printf("\t Primary vertex analysis is %s\n", GetPlugin(kPrimVtx)?"ON":"OFF");
675   printf("\t Combined pid analysis is %s\n", GetPlugin(kCombinedPid)?"ON":"OFF");
676   printf("\t HFE pid analysis is %s\n", GetPlugin(kHFEpid)?"ON":"OFF");
677   printf("\t Post process analysis is %s\n", GetPlugin(kPostProcess)?"ON":"OFF");
678   printf("\t ");
679   printf("\n");
680 }
681
682 //__________________________________________                                                  
683 void AliAnalysisTaskDCA::SwitchOnPlugin(Int_t plug){
684   //                                            
685   // Switch on Plugin          
686   // Available:                                  
687   //  - analyze impact parameter
688   //  - Post Processing  
689
690   AliDebug(2,Form("SwitchOnPlugin %d",plug));  
691
692   switch(plug){
693   case kPostProcess: 
694     SETBIT(fPlugins, plug); 
695     break;
696   case kImpactPar: 
697     SETBIT(fPlugins, plug); 
698     break;
699   case kPrimVtx: 
700     SETBIT(fPlugins, plug); 
701     break;
702   case kCombinedPid:
703     SETBIT(fPlugins, plug); 
704     break;
705   case kHFEpid:
706     SETBIT(fPlugins, plug); 
707     break;
708   case kKFdca:
709     SETBIT(fPlugins, plug); 
710     break;
711   default: 
712     AliError("Unknown Plugin");
713   };
714 }
715
716
717 //____________________________________________________________
718 void AliAnalysisTaskDCA::MakeParticleContainer(){
719
720   //printf("\n=====MakeParticleContainer=====\n");
721   //
722   // Create the particle container (borrowed from AliAnalysisTaskHFE)
723   //
724   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
725   const Double_t kPtmin = 0.1, kPtmax = 10.;
726   const Double_t kEtamin = -0.9, kEtamax = 0.9;
727   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
728
729   //arrays for the number of bins in each dimension
730   Int_t iBin[kNvar];
731   iBin[0] = 40; //bins in pt
732   iBin[1] =  8; //bins in eta 
733   iBin[2] = 18; // bins in phi
734
735   //arrays for lower bounds :
736   Double_t* binEdges[kNvar];
737   for(Int_t ivar = 0; ivar < kNvar; ivar++)
738     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
739
740   //values for bin lower bounds
741   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);  
742   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
743   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
744
745   //one "container" for MC
746   const Int_t kNcutStepsESDtrack = AliHFEcuts::kNcutStepsRecTrack + 1;
747   const Int_t kNcutStepsTrack = AliHFEcuts::kNcutStepsMCTrack + kNcutStepsESDtrack;
748   AliCFContainer* container = new AliCFContainer("container","container for tracks", (kNcutStepsTrack + 1 + 2*(kNcutStepsESDtrack + 1)), kNvar, iBin);
749
750   //setting the bin limits
751   for(Int_t ivar = 0; ivar < kNvar; ivar++)
752     container -> SetBinLimits(ivar, binEdges[ivar]);
753   fCFM->SetParticleContainer(container);
754
755   //create correlation matrix for unfolding
756   Int_t thnDim[2*kNvar];
757   for (int k=0; k<kNvar; k++) {
758     //first half  : reconstructed 
759     //second half : MC
760     thnDim[k]      = iBin[k];
761     thnDim[k+kNvar] = iBin[k];
762   }
763
764
765 }
766
767 //____________________________________________________________
768 void AliAnalysisTaskDCA::AddPIDdetector(TString detector){
769   
770   //
771   // Adding PID detector to the task
772   //
773   //printf("\n=====AddPIDdetector=====\n");
774   
775   if(!fPIDdetectors.Length()) 
776     fPIDdetectors = detector;
777   else
778     fPIDdetectors += ":" + detector;
779 }
780