PWGHF/vertexingHF/AliCFTaskForDStarAnalysis.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, 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 /* $Id$ */
  17
  18 //-----------------------------------------------------------------------
  19 //   Class for DStar corrections:
  20 //
  21 //   The D0 cutting varibles position in the container and container
  22 //   binning method from a C.Zampolli example
  23 //   In this way a simple comparison between D0 and D0 from D* is possible.
  24 //
  25 //-----------------------------------------------------------------------
  26 // Author : A.Grelli, Utrecht University
  27 //
  28 //         a.grelli@uu.nl
  29 //-----------------------------------------------------------------------
  30
  31 #include <TCanvas.h>
  32 #include <TParticle.h>
  33 #include <TDatabasePDG.h>
  34 #include <TH1I.h>
  35 #include <TStyle.h>
  36 #include <TFile.h>
  37
  38 #include "AliCFTaskForDStarAnalysis.h"
  39 #include "AliStack.h"
  40 #include "AliMCEvent.h"
  41 #include "AliCFManager.h"
  42 #include "AliCFContainer.h"
  43 #include "AliLog.h"
  44 #include "AliAnalysisManager.h"
  45 #include "AliAODHandler.h"
  46 #include "AliAODEvent.h"
  47 #include "AliAODRecoDecay.h"
  48 #include "AliAODRecoDecayHF.h"
  49 #include "AliAODRecoCascadeHF.h"
  50 #include "AliAODRecoDecayHF2Prong.h"
  51 #include "AliAODMCParticle.h"
  52 #include "AliAODMCHeader.h"
  53 #include "AliESDtrack.h"
  54 #include "TChain.h"
  55 #include "THnSparse.h"
  56 #include "TH2D.h"
  57
  58 //__________________________________________________________________________
  59 AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis() :
  60   AliAnalysisTaskSE(),
  61   fCFManager(0x0),
  62   fHistEventsProcessed(0x0),
  63   fCorrelation(0x0),
  64   fCountRecoDStarSel(0),
  65   fEvents(0),
  66   fMinITSClusters(5),
  67   fMinITSClustersSoft(4),
  68   fAcceptanceUnf(kTRUE)
  69 {
  70   //
  71   //Default ctor
  72   //
  73 }
  74 //___________________________________________________________________________
  75 AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const Char_t* name) :
  76   AliAnalysisTaskSE(name),
  77   fCFManager(0x0),
  78   fHistEventsProcessed(0x0),
  79   fCorrelation(0x0),
  80   fCountRecoDStarSel(0),
  81   fEvents(0),
  82   fMinITSClusters(5),
  83   fMinITSClustersSoft(4),
  84   fAcceptanceUnf(kTRUE)
  85 {
  86   //
  87   // Constructor. Initialization of Inputs and Outputs
  88   //
  89   Info("AliCFTaskForDStarAnalysis","Calling Constructor");
  90
  91   DefineOutput(1,TH1I::Class());
  92   DefineOutput(2,AliCFContainer::Class());
  93   DefineOutput(3,THnSparseD::Class());
  94 }
  95
  96 //___________________________________________________________________________
  97 AliCFTaskForDStarAnalysis& AliCFTaskForDStarAnalysis::operator=(const AliCFTaskForDStarAnalysis& c)
  98 {
  99   //
 100   // Assignment operator
 101   //
 102   if (this!=&c) {
 103     AliAnalysisTaskSE::operator=(c) ;
 104     fCFManager  = c.fCFManager;
 105     fHistEventsProcessed = c.fHistEventsProcessed;
 106   }
 107   return *this;
 108 }
 109
 110 //___________________________________________________________________________
 111 AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const AliCFTaskForDStarAnalysis& c) :
 112   AliAnalysisTaskSE(c),
 113   fCFManager(c.fCFManager),
 114   fHistEventsProcessed(c.fHistEventsProcessed),
 115   fCorrelation(c.fCorrelation),
 116   fCountRecoDStarSel(c.fCountRecoDStarSel),
 117   fEvents(c.fEvents),
 118   fMinITSClusters(c.fMinITSClusters),
 119   fMinITSClustersSoft(c.fMinITSClustersSoft),
 120   fAcceptanceUnf(c.fAcceptanceUnf)
 121 {
 122   //
 123   // Copy Constructor
 124   //
 125 }
 126
 127 //___________________________________________________________________________
 128 AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
 129   //
 130   //destructor
 131   //
 132   if (fCFManager)           delete fCFManager ;
 133   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 134   if (fCorrelation)         delete fCorrelation ;
 135 }
 136
 137 //_________________________________________________
 138 void AliCFTaskForDStarAnalysis::UserExec(Option_t *)
 139 {
 140   //
 141   // Main loop function
 142   //
 143
 144   if (!fInputEvent) {
 145     Error("UserExec","NO EVENT FOUND!");
 146     return;
 147   }
 148
 149   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 150
 151   TClonesArray *arrayDStartoD0pi=0;
 152
 153   if(!aodEvent && AODEvent() && IsStandardAOD()) {
 154     // In case there is an AOD handler writing a standard AOD, use the AOD
 155     // event in memory rather than the input (ESD) event.
 156     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 157     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 158     // have to taken from the AOD event hold by the AliAODExtension
 159     AliAODHandler* aodHandler = (AliAODHandler*)
 160       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 161     if(aodHandler->GetExtensions()) {
 162       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 163       AliAODEvent *aodFromExt = ext->GetAOD();
 164       arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 165     }
 166   } else {
 167     arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 168   }
 169
 170   if (!arrayDStartoD0pi) {
 171     AliError("Could not find array of HF vertices");
 172     return;
 173   }
 174
 175   // fix for temporary bug in ESDfilter
 176   // the AODs with null vertex pointer didn't pass the PhysSel
 177   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 178
 179   fEvents++;
 180   if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
 181
 182   fCFManager->SetRecEventInfo(aodEvent);
 183   fCFManager->SetMCEventInfo(aodEvent);
 184
 185   // event selection
 186   Double_t containerInput[15] ;
 187   Double_t containerInputMC[15] ;
 188
 189   //loop on the MC event
 190
 191   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 192   if (!mcArray) {
 193     AliError("Could not find Monte-Carlo in AOD");
 194     return;
 195   }
 196
 197   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 198   if (!mcHeader) {
 199     AliError("Could not find MC Header in AOD");
 200     return;
 201   }
 202
 203   // AOD primary vertex
 204   AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 205   Double_t zPrimVertex = vtx1->GetZ();
 206   Double_t zMCVertex = mcHeader->GetVtxZ();
 207   Bool_t vtxFlag = kTRUE;
 208   TString title=vtx1->GetTitle();
 209   if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
 210
 211   for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 212     AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 213     if (!mcPart) {
 214       AliWarning("MC Particle not found in tree, skipping");
 215       continue;
 216     }
 217
 218     // check the MC-level cuts
 219     if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
 220
 221     // fill the container for Gen-level selection
 222     Double_t vectorMC[9] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
 223
 224     if (GetDStarMCParticle(mcPart, mcArray, vectorMC)){
 225
 226       containerInputMC[0] = vectorMC[0];
 227       containerInputMC[1] = vectorMC[1] ;
 228       containerInputMC[2] = vectorMC[2] ;
 229       containerInputMC[3] = vectorMC[3] ;
 230       containerInputMC[4] = vectorMC[4] ;
 231       containerInputMC[5] = vectorMC[5] ;
 232       containerInputMC[6] = 0.;
 233       containerInputMC[7] = 0.;
 234       containerInputMC[8] = 0.;
 235       containerInputMC[9] = -100000.;
 236       containerInputMC[10] = 1.01;
 237       containerInputMC[11] = vectorMC[6];
 238       containerInputMC[12] = zMCVertex;     // z vertex
 239       containerInputMC[13] = vectorMC[7];   // pt D0 pion
 240       containerInputMC[14] = vectorMC[8];   // pt D0 kaon
 241
 242       fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
 243
 244       // check the MC-Acceptance level cuts
 245       // since standard CF functions are not applicable, using Kine Cuts on daughters
 246
 247       Int_t daughter0 = mcPart->GetDaughter(0);
 248       Int_t daughter1 = mcPart->GetDaughter(1);
 249
 250       AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
 251
 252       //D0
 253       AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 254       if (!mcPartDaughter0) {
 255         AliError("Could not find Monte-Carlo in AOD");
 256         return;
 257       }
 258
 259       // Soft Pion
 260       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 261       if (!mcPartDaughter1) {
 262         AliError("Could not find Monte-Carlo in AOD");
 263         return;
 264       }
 265
 266       // Acceptance variables for the soft pion
 267       Double_t eta1 = mcPartDaughter1->Eta();
 268       Double_t pt1  = mcPartDaughter1->Pt();
 269
 270       Int_t daughD00 = 0;
 271       Int_t daughD01 = 0;
 272
 273       // Just to be sure to take the right particles
 274       if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
 275         daughD00 = mcPartDaughter0->GetDaughter(0);
 276         daughD01 = mcPartDaughter0->GetDaughter(1);
 277       }else{
 278         daughD00 = mcPartDaughter1->GetDaughter(0);
 279         daughD01 = mcPartDaughter1->GetDaughter(1);
 280       }
 281
 282       // the D0 daughters
 283       AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
 284       AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
 285
 286       if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
 287         AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
 288         return;
 289       }
 290
 291       // D0 daughters - needed for acceptance
 292       Double_t theD0pt0 =  mcD0PartDaughter0->Pt();
 293       Double_t theD0pt1 =  mcD0PartDaughter1->Pt();
 294       Double_t theD0eta0 = mcD0PartDaughter0->Eta();
 295       Double_t theD0eta1 = mcD0PartDaughter1->Eta();
 296
 297       // ACCEPTANCE REQUESTS ---------
 298
 299       // soft pion
 300       Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
 301       // Do daughters
 302       Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1);
 303       Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
 304
 305       if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
 306
 307         AliDebug(2, "D* Daughter particles in acceptance");
 308
 309         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
 310
 311         // check on the vertex
 312         if (vtxFlag){
 313           // filling the container if the vertex is ok
 314           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
 315
 316           Bool_t refitFlag = kTRUE;
 317           for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
 318             AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
 319
 320             // refit only for D0 daughters
 321             if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
 322               if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
 323                 refitFlag = kFALSE;
 324               }
 325             }
 326           }
 327           if (refitFlag){ // refit
 328
 329             fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
 330
 331           } // end of refit
 332         } // end of vertex
 333       } //end of acceptance
 334     } // end of MC D*
 335     else {
 336       AliDebug(3,"Problems in filling the container");
 337       continue;
 338     }
 339   } // end of MC loop
 340
 341   //rec
 342   AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
 343
 344   //D* and D0 prongs needed to MatchToMC method
 345   Int_t pdgDgDStartoD0pi[2]={421,211};
 346   Int_t pdgDgD0toKpi[2]={321,211};
 347
 348   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
 349
 350     // D* candidates
 351     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
 352
 353     // D0 from the reco cascade
 354     AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 355     Bool_t unsetvtx=kFALSE;
 356
 357     // needed for pointing angle
 358     if(!D0particle->GetOwnPrimaryVtx()) {
 359       D0particle->SetOwnPrimaryVtx(vtx1);
 360       unsetvtx=kTRUE;
 361     }
 362
 363     // find associated MC particle for D* ->D0toKpi
 364     Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 365
 366     // find D0->Kpi ... needed in the following
 367     Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
 368
 369     if (mcLabel == -1 || mcD0Label ==-1)
 370       {
 371         AliDebug(2,"No MC particle found");
 372       }
 373     else {
 374
 375       // the D* and the D0 in MC
 376       AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
 377       AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
 378
 379       if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
 380         AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
 381         continue;
 382       }
 383
 384       // soft pion
 385       AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
 386
 387       //D0tokpi
 388       AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
 389       AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
 390
 391       // check if associated MC v0 passes the cuts
 392       if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) {
 393         AliDebug(2, "Skipping the particles due to cuts");
 394         continue;
 395       }
 396
 397       // fill the container...
 398       Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
 399       Double_t rapidity = dstarD0pi->YDstar();
 400       Double_t cosThetaStar = 9999.;
 401       Double_t pTpi  = 0.;
 402       Double_t pTK   = 0.;
 403       Double_t dca   = (D0particle->GetDCA())*1E4;
 404       Double_t d0pi  = 0.;
 405       Double_t d0K   = 0.;
 406       Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
 407       Double_t cosPointingAngle = D0particle->CosPointingAngle();
 408       Double_t phi = D0particle->Phi();
 409
 410       // Select D0 cutting variables
 411       Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
 412
 413       // D0 related variables
 414       if (pdgCode > 0){
 415
 416         cosThetaStar = D0particle->CosThetaStarD0();
 417         pTpi    = D0particle->PtProng(0);
 418         pTK     = D0particle->PtProng(1);
 419         d0pi    = (D0particle->Getd0Prong(0))*1E4;
 420         d0K     = (D0particle->Getd0Prong(1))*1E4;
 421
 422       }
 423       else {
 424
 425         cosThetaStar = D0particle->CosThetaStarD0bar();
 426         pTpi    = D0particle->PtProng(1);
 427         pTK     = D0particle->PtProng(0);
 428         d0pi    = (D0particle->Getd0Prong(1))*1E4;
 429         d0K     = (D0particle->Getd0Prong(0))*1E4;
 430
 431       }
 432
 433       // ct of the D0 from D*
 434       Double_t cT = D0particle->CtD0();
 435
 436       containerInput[0]  = pt;
 437       containerInput[1]  = rapidity;
 438       containerInput[2]  = cosThetaStar;
 439       containerInput[3]  = track2->Pt();
 440       containerInput[4]  = D0particle->Pt();
 441       containerInput[5]  = cT*1.E4;  // in micron
 442       containerInput[6]  = dca;  // in micron
 443       containerInput[7]  = d0pi;  // in micron
 444       containerInput[8]  = d0K;  // in micron
 445       containerInput[9]  = d0xd0;  // in micron^2
 446       containerInput[10] = cosPointingAngle;  // in micron
 447       containerInput[11] = phi;
 448       containerInput[12] = zPrimVertex;    // z of reconstructed of primary vertex
 449       containerInput[13] = pTpi;  // D0 pion
 450       containerInput[14] = pTK;   // D0 kaon
 451
 452       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
 453
 454       // refit in ITS and TPC for D0 daughters
 455       if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
 456         continue;
 457       }
 458
 459       // reft in ITS for soft pion
 460       if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
 461         continue;
 462       }
 463
 464       // cut in acceptance for the soft pion and for the D0 daughters
 465       Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
 466       Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
 467
 468       // soft pion acceptance ... is it fine 0.9?????
 469       Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
 470
 471       if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
 472         AliDebug(2,"D* reco daughters in acceptance");
 473         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
 474
 475         if(fAcceptanceUnf){
 476           Double_t fill[4]; //fill response matrix
 477
 478           // dimensions 0&1 : pt,eta (Rec)
 479           fill[0] = pt ;
 480           fill[1] = rapidity;
 481           // dimensions 2&3 : pt,eta (MC)
 482           fill[2] =  mcVtxHFDStar->Pt();
 483           fill[3] =  mcVtxHFDStar->Y();
 484           fCorrelation->Fill(fill);
 485         }
 486
 487         // cut on the min n. of clusters in ITS for the D0 and soft pion
 488         Int_t ncls0=0,ncls1=0,ncls2=0;
 489         for(Int_t l=0;l<6;l++) {
 490           if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
 491           if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
 492           if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
 493         }
 494         // see AddTask for soft pion ITS clusters request
 495         AliDebug(2, Form("n clusters = %d", ncls0));
 496
 497         if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
 498           fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
 499
 500           // D0 cuts optimized for D* analysis
 501           Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
 502
 503           // needed for cuts
 504           Double_t theD0pt = D0particle->Pt();
 505
 506           if (theD0pt <= 1){ // first bin not optimized
 507             cuts[0] = 400;
 508             cuts[1] = 0.8;
 509             cuts[2] = 0.21;
 510             cuts[3] = 500;
 511             cuts[4] = 500;
 512             cuts[5] = -20000;
 513             cuts[6] = 0.6;
 514           }
 515           else if (theD0pt > 1 && theD0pt <= 2){
 516             cuts[0] = 200;
 517             cuts[1] = 0.7;
 518             cuts[2] = 0.8;
 519             cuts[3] = 210;
 520             cuts[4] = 210;
 521             cuts[5] = -20000;
 522             cuts[6] = 0.9;
 523           }
 524           else if (theD0pt > 2 && theD0pt <= 3){
 525             cuts[0] = 400;
 526             cuts[1] = 0.8;
 527             cuts[2] = 0.8;
 528             cuts[3] = 420;
 529             cuts[4] = 350;
 530             cuts[5] = -8500;
 531             cuts[6] = 0.9;
 532           }
 533           else if (theD0pt > 3 && theD0pt <= 5){
 534             cuts[0] = 160;
 535             cuts[1] = 1.0;
 536             cuts[2] = 1.2;
 537             cuts[3] = 560;
 538             cuts[4] = 420;
 539             cuts[5] = -8500;
 540             cuts[6] = 0.9;
 541           }
 542           else if (theD0pt > 5){
 543             cuts[0] = 800;
 544             cuts[1] = 1.0;
 545             cuts[2] = 1.2;
 546             cuts[3] = 700;
 547             cuts[4] = 700;
 548             cuts[5] = 10000;
 549             cuts[6] = 0.9;
 550           }
 551           if (dca < cuts[0]
 552               && TMath::Abs(cosThetaStar) < cuts[1]
 553               && pTpi > cuts[2]
 554               && pTK > cuts[2]
 555               && TMath::Abs(d0pi) < cuts[3]
 556               && TMath::Abs(d0K) < cuts[4]
 557               && d0xd0 < cuts[5]
 558               && cosPointingAngle > cuts[6]
 559             ){
 560
 561             AliDebug(2,"Particle passed D* selection cuts");
 562             fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;
 563
 564             if(!fAcceptanceUnf){ // unfolding
 565
 566               Double_t fill[4]; //fill response matrix
 567
 568               // dimensions 0&1 : pt,eta (Rec)
 569               fill[0] = pt ;
 570               fill[1] = rapidity;
 571               // dimensions 2&3 : pt,eta (MC)
 572               fill[2] =  mcVtxHFDStar->Pt();
 573               fill[3] =  mcVtxHFDStar->Y();
 574
 575               fCorrelation->Fill(fill);
 576             }
 577           }
 578         }
 579       }
 580     }
 581     if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
 582   } // end loop on D*->Kpipi
 583
 584   fHistEventsProcessed->Fill(0);
 585
 586   PostData(1,fHistEventsProcessed) ;
 587   PostData(2,fCFManager->GetParticleContainer()) ;
 588   PostData(3,fCorrelation) ;
 589 }
 590
 591
 592 //___________________________________________________________________________
 593 void AliCFTaskForDStarAnalysis::Terminate(Option_t*)
 594 {
 595   // The Terminate()
 596   AliAnalysisTaskSE::Terminate();
 597
 598   // draw correlation matrix
 599   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 600   if(!cont) {
 601     printf("CONTAINER NOT FOUND\n");
 602     return;
 603   }
 604 }
 605
 606 //___________________________________________________________________________
 607 void AliCFTaskForDStarAnalysis::UserCreateOutputObjects() {
 608
 609   //
 610   // useroutput
 611   //
 612
 613   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 614
 615   //slot #1
 616   OpenFile(1);
 617   fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
 618 }
 619
 620 //________________________________________________________________________________
 621 Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
 622
 623   //
 624   // fill the D* and D0 MC container
 625   //
 626
 627   Bool_t isDStar = kFALSE;
 628
 629   if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
 630
 631   // getting the daughters
 632   Int_t daughter0 = mcPart->GetDaughter(0);
 633   Int_t daughter1 = mcPart->GetDaughter(1);
 634
 635   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
 636   if (daughter0 == 0 || daughter1 == 0) {
 637     AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
 638     return isDStar;
 639   }
 640
 641   if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
 642     AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
 643     return isDStar;
 644   }
 645
 646   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 647   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 648   if (!mcPartDaughter0 || !mcPartDaughter1) {
 649     AliWarning("D*: At least one Daughter Particle not found in tree, skipping");
 650     return isDStar;
 651   }
 652
 653   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
 654         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
 655       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
 656         TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
 657     AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
 658     return isDStar;
 659   }
 660
 661   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 662   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 663
 664   // getting vertex from daughters
 665   mcPartDaughter0->XvYvZv(vtx2daughter0);
 666   mcPartDaughter1->XvYvZv(vtx2daughter1);
 667
 668   // check if the secondary vertex is the same for both
 669   if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
 670     AliError("The D* daughters have different secondary vertex, skipping the track");
 671     return isDStar;
 672   }
 673
 674   AliAODMCParticle* neutralDaugh = mcPartDaughter0;
 675
 676   Double_t VectorD0[2] ={0.,0.};
 677
 678   if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
 679     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 680     return isDStar;
 681   }
 682   // get the pT of the daughters
 683
 684   Double_t pTpi = 0.;
 685   Double_t pTD0 = 0.;
 686
 687   if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
 688     pTpi = mcPartDaughter0->Pt();
 689     pTD0 = mcPartDaughter1->Pt();
 690   }
 691   else {
 692     pTpi = mcPartDaughter1->Pt();
 693     pTD0 = mcPartDaughter0->Pt();
 694   }
 695
 696   vectorMC[0] = mcPart->Pt();
 697   vectorMC[1] = mcPart->Y() ;
 698   vectorMC[2] = 0;
 699   vectorMC[3] = pTpi ;
 700   vectorMC[4] = pTD0 ;
 701   vectorMC[5] = 0;
 702   vectorMC[6] = mcPart->Phi() ;
 703   vectorMC[7] = VectorD0[0] ;
 704   vectorMC[8] = VectorD0[1] ;
 705
 706   isDStar = kTRUE;
 707
 708   return isDStar;
 709 }
 710 //________________________________________________________________________________________________
 711
 712 Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
 713
 714   //
 715   // chack wether D0 is decaing into kpi
 716   //
 717
 718   Bool_t isHadronic = kFALSE;
 719
 720   Int_t daughterD00 = neutralDaugh->GetDaughter(0);
 721   Int_t daughterD01 = neutralDaugh->GetDaughter(1);
 722
 723   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
 724   if (daughterD00 == 0 || daughterD01 == 0) {
 725     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 726     return isHadronic;
 727   }
 728
 729   if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
 730     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 731     return isHadronic;
 732   }
 733
 734   AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
 735   AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
 736   if (!mcPartDaughterD00 || !mcPartDaughterD01) {
 737     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
 738     return isHadronic;
 739   }
 740
 741   if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
 742         TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
 743       !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
 744         TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
 745     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
 746     return isHadronic;
 747   }
 748
 749   Double_t pTD0pi = 0;
 750   Double_t pTD0K = 0;
 751
 752
 753   if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
 754     pTD0pi = mcPartDaughterD00->Pt();
 755     pTD0K = mcPartDaughterD01->Pt();
 756   }
 757   else {
 758     pTD0pi = mcPartDaughterD01->Pt();
 759     pTD0K  = mcPartDaughterD00->Pt();
 760   }
 761
 762   isHadronic = kTRUE;
 763
 764   VectorD0[0] = pTD0pi;
 765   VectorD0[1] = pTD0K;
 766
 767   return isHadronic;
 768
 769 }