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