]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
dNdPt analysis
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.cxx
index d76d08caae54ea89c35a2031e6623a493ee128b2..22f9a0a85335f4a970df603a762fb9598f39d38b 100644 (file)
@@ -41,6 +41,8 @@
 #include "AliCFManager.h"
 #include "AliCFContainer.h"
 #include "AliLog.h"
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODRecoDecay.h"
 #include "AliAODRecoDecayHF.h"
@@ -174,9 +176,37 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
        }
        
        fEvents++;
-       if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
+       if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
        AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
-       fCFManager->SetEventInfo(aodEvent);
+
+       TClonesArray *arrayD0toKpi=0;
+
+       if(!aodEvent && AODEvent() && IsStandardAOD()) {
+         // In case there is an AOD handler writing a standard AOD, use the AOD 
+         // event in memory rather than the input (ESD) event.    
+         aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
+         // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+         // have to taken from the AOD event hold by the AliAODExtension
+         AliAODHandler* aodHandler = (AliAODHandler*) 
+           ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+         if(aodHandler->GetExtensions()) {
+           AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+           AliAODEvent *aodFromExt = ext->GetAOD();
+           arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+         }
+       } else {
+         arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
+       }
+
+
+       if (!arrayD0toKpi) {
+         AliError("Could not find array of HF vertices");
+         return;
+       }
+
+
+       fCFManager->SetRecEventInfo(aodEvent);
+       fCFManager->SetMCEventInfo(aodEvent);
        
        // MC-event selection
        Double_t containerInput[12] ;
@@ -273,10 +303,10 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
                                AliDebug(2, "Daughter particles in acceptance");
                                if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
-                                       AliInfo("Inconsistency with CF cut for daughter 0!");
+                                       AliDebug(2,"Inconsistency with CF cut for daughter 0!");
                                } 
                                if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
-                                       AliInfo("Inconsistency with CF cut for daughter 1!");
+                                       AliDebug(2,"Inconsistency with CF cut for daughter 1!");
                                } 
                                fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
                                icountAcc++;
@@ -289,11 +319,13 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                        printf("Vertex cut passed 2\n");
                                        icountVertex++;
                                        // check on the kTPCrefit and kITSrefit conditions of the daughters
-                                       Bool_t refitFlag = kFALSE;
+                                       Bool_t refitFlag = kTRUE;
                                        for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
                                                AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
-                                               if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1) && (((track->GetStatus()&AliESDtrack::kTPCrefit)==AliESDtrack::kTPCrefit) && ((track->GetStatus()&AliESDtrack::kITSrefit)==AliESDtrack::kITSrefit))){
-                                                       refitFlag = kTRUE;
+                                               if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
+                                                 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
+                                                   refitFlag = kFALSE;
+                                                 }
                                                }
                                        }
                                        if (refitFlag){
@@ -321,18 +353,15 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
 
        if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
        
-       AliInfo(Form("Found %i MC particles that are D0!!",icountMC));
-       AliInfo(Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
-       AliInfo(Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
-       AliInfo(Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
+       AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
+       AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
 
        // Now go to rec level
        fCountMC += icountMC;
        fCountAcc += icountAcc;
 
-       // load heavy flavour vertices
-       TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));     
-       if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
        AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
 
        Int_t pdgDgD0toKpi[2]={321,211};
@@ -394,7 +423,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                        // fill the container...
                        Double_t pt = d0tokpi->Pt();
                        Double_t rapidity = d0tokpi->YD0();
-                       
+                       Double_t invMass=0.;
                        Double_t cosThetaStar = 9999.;
                        Double_t pTpi = 0.;
                        Double_t pTK = 0.;
@@ -411,6 +440,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                pTK = d0tokpi->PtProng(1);
                                d0pi = d0tokpi->Getd0Prong(0);
                                d0K = d0tokpi->Getd0Prong(1);
+                               invMass=d0tokpi->InvMassD0();
                        }
                        else {
                                cosThetaStar = d0tokpi->CosThetaStarD0bar();
@@ -418,6 +448,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                pTK = d0tokpi->PtProng(0);
                                d0pi = d0tokpi->Getd0Prong(1);
                                d0K = d0tokpi->Getd0Prong(0);
+                               invMass=d0tokpi->InvMassD0bar();
                        }
 
                        Double_t cT = d0tokpi->CtD0();
@@ -490,57 +521,122 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                }  
                                
                                // cut on the min n. of clusters in ITS
-                               Int_t ncls0=0;
-                               for(Int_t l=0;l<6;l++) if(TESTBIT(d0tokpi->GetITSClusterMap(),l)) ncls0++;
+                               Int_t ncls0=0,ncls1=0;
+                               for(Int_t l=0;l<6;l++) {
+                                 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
+                                 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
+                               }
                                AliDebug(2, Form("n clusters = %d", ncls0));
-                               if (ncls0 >= fMinITSClusters){
+                               if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters) {
                                        fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
                                        icountRecoITSClusters++;   
                                        AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
                                        
                                        // PPR cuts 
-                                       Double_t cuts[6] = {9999999., 1.1, 0., 9999999., 9999999., 0.};
+                                       //added mass cut for D0 analysis
+                                       Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027}; 
+                                       //cuts of the D0 analysis (looser) see AliAnalysisTaskSED0Mass.cxx
                                        if (pt <= 1){
-                                               cuts[0] = 400;
-                                               cuts[1] = 0.8;
-                                               cuts[2] = 0.5;
-                                               cuts[3] = 500;
-                                               cuts[4] = -20000;
-                                               cuts[5] = 0.5;
+                                         
+                                         //coutputmassD02
+                                         cuts[0] = 400; //dca (um)
+                                         cuts[1] = 0.8; //cosTstar
+                                         cuts[2] = 0.5; //pt  (GeV/c)
+                                         cuts[3] = 500; //d0  (um)
+                                         cuts[4] = -25000; //d0xd0 (um^2)
+                                         cuts[5] = 0.7; //cosTpointing
+                                         /*
+                                         //PPR
+                                         cuts[0] = 400; //dca (um)
+                                         cuts[1] = 0.8; //cosTstar
+                                         cuts[2] = 0.5; //pt  (GeV/c)
+                                         cuts[3] = 500; //d0  (um)
+                                         cuts[4] = -20000; //d0xd0 (um^2)
+                                         cuts[5] = 0.5; //cosTpointing
+                                         */
                                        }
+                                       /* //not same cuts for pt = 1 to 2 and 1 to 3 GeV in case must match them because of poor stat
                                        else if (pt > 1 && pt <= 2){
-                                               cuts[0] = 300;
-                                               cuts[1] = 0.8;
-                                               cuts[2] = 0.6;
-                                               cuts[3] = 500;
-                                               cuts[4] = -20000;
-                                               cuts[5] = 0.6;
-                                       }
-                                       else if (pt > 2 && pt <= 3){
-                                               cuts[0] = 200;
-                                               cuts[1] = 0.8;
-                                               cuts[2] = 0.7;
-                                               cuts[3] = 500;
-                                               cuts[4] = -20000;
-                                               cuts[5] = 0.8;
+                                        
+                                         //coutputmassD02
+                                         cuts[0] = 300;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.6;
+                                         cuts[3] = 1000;
+                                         cuts[4] = -25000;
+                                         cuts[5] = 0.7;
+                                         
+                                         //PPR
+                                         cuts[0] = 300;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.6;
+                                         cuts[3] = 500;
+                                         cuts[4] = -20000;
+                                         cuts[5] = 0.6;
+                                         
+                                         }
+                                       */
+                                       else if (pt > 1 && pt <= 3){
+                                         
+                                         //coutputmassD02
+                                         cuts[0] = 200;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.7;
+                                         cuts[3] = 1000;
+                                         cuts[4] = -25000;
+                                         cuts[5] = 0.8;
+                                         /*
+                                         //PPR
+                                         cuts[0] = 300;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.6;
+                                         cuts[3] = 500;
+                                         cuts[4] = -20000;
+                                         cuts[5] = 0.6;
+                                         */
                                        }
                                        else if (pt > 3 && pt <= 5){
-                                               cuts[0] = 200;
-                                               cuts[1] = 0.8;
-                                               cuts[2] = 0.7;
-                                               cuts[3] = 500;
-                                               cuts[4] = -10000;
-                                               cuts[5] = 0.8;
+                                         
+                                         //coutputmassD02
+                                         cuts[0] = 200;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.7;
+                                         cuts[3] = 500;
+                                         cuts[4] = -15000;
+                                         cuts[5] = 0.8;
+                                         /*
+                                         //PPR
+                                         cuts[0] = 200;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.7;
+                                         cuts[3] = 500;
+                                         cuts[4] = -10000;
+                                         cuts[5] = 0.8;
+                                         */
                                        }
                                        else if (pt > 5){
-                                               cuts[0] = 200;
-                                               cuts[1] = 0.8;
-                                               cuts[2] = 0.7;
-                                               cuts[3] = 500;
-                                               cuts[4] = -5000;
-                                               cuts[5] = 0.8;
+                                         
+                                         //coutputmassD02
+                                         cuts[0] = 200;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.7;
+                                         cuts[3] = 500;
+                                         cuts[4] = -15000;
+                                         cuts[5] = 0.9;
+                                         /*
+                                         //PPR
+                                         cuts[0] = 200;
+                                         cuts[1] = 0.8;
+                                         cuts[2] = 0.7;
+                                         cuts[3] = 500;
+                                         cuts[4] = -5000;
+                                         cuts[5] = 0.8;
+                                         */
                                        }
-                                       if (dca*1E4 < cuts[0] 
+                               
+                                       Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+                                       if (TMath::Abs(invMass-mD0PDG) < cuts[6]
+                                           && dca*1E4 < cuts[0] 
                                            && TMath::Abs(cosThetaStar) < cuts[1]  
                                            && pTpi > cuts[2] 
                                            && pTK > cuts[2]  
@@ -550,7 +646,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
                                            && cosPointingAngle > cuts[5]
                                            ){
                                                
-                                               AliDebug(2,"Particle passed PPR cuts");
+                                               AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
                                                fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;   
                                                icountRecoPPR++;
                                                
@@ -1031,7 +1127,7 @@ void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
         corr2->Draw("text");
       
 
-       TFile* file_projection = new TFile("file_projection.root","RECREATE");
+       TFile* file_projection = new TFile("CFtaskHFprojection.root","RECREATE");
 
         corr1->Write();
         corr2->Write();