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