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