Update from Prabhat - Dpt Dpt corr
authorsjena <sjena@cern.ch>
Fri, 27 Jun 2014 14:48:08 +0000 (16:48 +0200)
committersjena <sjena@cern.ch>
Fri, 27 Jun 2014 14:48:08 +0000 (16:48 +0200)
PWGCF/Correlations/DPhi/AliDptDptInMC.cxx
PWGCF/Correlations/DPhi/AliDptDptInMC.h
PWGCF/Correlations/macros/dptdptcorrelations/AddTaskCorrMC.C

index 8e54467..a7d419a 100644 (file)
@@ -52,7 +52,7 @@
 #include "AliDptDptInMC.h"
 
 #include "AliPID.h"
-//#include "AliPIDResponse.h"
+#include "AliPIDResponse.h"
 
 #include "AliESDVertex.h"
 #include "AliESDEvent.h"
@@ -75,7 +75,7 @@
 #include "AliGenEventHeader.h"
 #include "AliPID.h"
 #include "AliAODPid.h"
-//#include "AliPIDResponse.h"
+#include "AliPIDResponse.h"
 #include "AliAODpidUtil.h"
 #include "AliPIDCombined.h"
 
@@ -89,7 +89,7 @@ AliDptDptInMC::AliDptDptInMC()
 fAODEvent(0), 
 fESDEvent(0),             //! ESD Event
 fInputHandler(0),
-//fPIDResponse(0x0),
+fPIDResponse(0x0),
 _outputHistoList(0),
 fArrayMC (0),
 _twoPi         ( 2.0 * 3.1415927),
@@ -384,7 +384,7 @@ AliDptDptInMC::AliDptDptInMC(const TString & name)
 fAODEvent(0), 
 fESDEvent(0),  
 fInputHandler(0),
-  //fPIDResponse(0),
+fPIDResponse(0),
 _outputHistoList(0),
 fArrayMC (0),
 _twoPi         ( 2.0 * 3.1415927),
@@ -781,6 +781,7 @@ void AliDptDptInMC::UserCreateOutputObjects()
   //_dedx_1     = new float[arraySize];
   
   __n1_1_vsPt              = getDoubleArray(_nBins_pt_1,        0.);
+
   __n1_1_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_1,    0.);
   __s1pt_1_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_1,    0.);
   __n1_1_vsZEtaPhiPt       = getFloatArray(_nBins_zEtaPhiPt_1,  0.);
@@ -1001,10 +1002,14 @@ void  AliDptDptInMC::createHistograms()
   name = "zV"; _vertexZ = createHisto1D(name,name,_nBins_vertexZ, _min_vertexZ, _max_vertexZ, "z-Vertex (cm)", _title_counts);
   
 
-  name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
-  name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
-  name = "DCAz";    _dcaz     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
-  name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
+  //name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
+  //name = "Phi";     _phidis   = createHisto1F(name,name, 360, 0.0, 6.4, "#phi","counts");
+  //name = "DCAz";    _dcaz     = createHisto1F(name,name, 340, -3.3, 3.3, "dcaZ","counts");
+  //name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 100, -0.1, 2.5, "dcaXY","counts");
+
+  name = "Eta";     _etadis   = createHisto1F(name,name, 250, 0.0, 2.5, "#eta","counts"); //temporaryly
+  name = "Phi";     _phidis   = createHisto1F(name,name, 250, 0.0, 2.5, "#phi","counts");
+  name = "DCAz";    _dcaz     = createHisto1F(name,name, 250, 0.0, 2.5, "dcaZ","counts");
 
   //name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
   //name = "Nclus2";   _Ncluster2    = createHisto1F(name,name, 200, 0, 200, "Ncluster2","counts");
@@ -1012,6 +1017,7 @@ void  AliDptDptInMC::createHistograms()
   if (_singlesOnly)
     {
     name = n1Name+part_1_Name+vsPt;              _n1_1_vsPt              = createHisto1F(name,name, _nBins_pt_1,  _min_pt_1,  _max_pt_1,   _title_pt_1,  _title_AvgN_1);
+
     name = n1Name+part_1_Name+vsZ+vsEtaPhi+vsPt; _n1_1_vsZVsEtaVsPhiVsPt = createHisto3F(name,name, _nBins_vertexZ,_min_vertexZ,_max_vertexZ, _nBins_etaPhi_1, 0., double(_nBins_etaPhi_1), _nBins_pt_1, _min_pt_1, _max_pt_1, "zVertex", _title_etaPhi_1,  _title_pt_1);
     //name = "dedxVsP_1";                          _dedxVsP_1              = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
     //name = "corrDedxVsP_1";                      _corrDedxVsP_1          = createHisto2F(name,name,400,-2.,2.,120,0.,120.,"p (GeV/c)", "dedx", "counts");
@@ -1150,29 +1156,11 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
     {
       return;
     }
-  //  fPIDResponse =inputHandler->GetPIDResponse();
-  //if (!fPIDResponse){
-  //AliFatal("This Task needs the PID response attached to the inputHandler");
-  //return;
-  //}
-
-  //--------------------------- //New Chage
-  fArrayMC = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
-  if (!fArrayMC) {
-    AliFatal("Error: MC particles branch not found!\n");
-    return;
-  }
-
-  AliAODMCHeader *mcHdr=NULL;
-  mcHdr=(AliAODMCHeader*)fAODEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
-  if(!mcHdr) {
-    printf("MC header branch not found!\n");
+  fPIDResponse =inputHandler->GetPIDResponse();
+  if (!fPIDResponse){
+    AliFatal("This Task needs the PID response attached to the inputHandler");
     return;
   }
-
-  //------------------------
-
-
   // count all events looked at here                                                                                                        
   _eventCount++;
 
@@ -1263,13 +1251,13 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
       //MC AOD Truth 
       if(fAnalysisType == "MCAOD")
         { 
-        _nTracks = fArrayMC->GetEntriesFast();
+         AliMCEvent* mcEvent = MCEvent();
+         _nTracks = mcEvent->GetNumberOfTracks();
           _mult3    = _nTracks;
-
           //loop over tracks starts here                                                                                                       
           for (int iTrack=0; iTrack< _nTracks; iTrack++)
             {
-             AliAODMCParticle *aodTrack = (AliAODMCParticle*) fArrayMC->At(iTrack);
+              AliAODMCParticle *aodTrack = (AliAODMCParticle*) mcEvent->GetTrack(iTrack);
 
               if (!aodTrack)
                 {
@@ -1277,7 +1265,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                   continue;
                 }
              
-             //if(!aodTrack->IsPhysicalPrimary()) continue;
+             if(!aodTrack->IsPhysicalPrimary()) continue;
 
              q      = aodTrack->Charge();
               charge = int(q);
@@ -1290,23 +1278,47 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              // Kinematics cuts from ESD track cuts                                                           
              if( pt < 0.2 || pt > 2.0)      continue;
              if( eta < -0.8 || eta > 0.8)  continue;
-             
-             // Remove neutral tracks                                                                                         
-              if(q == 0) continue;  
+
+             if(q == 0) continue;  
+             _dcaz->Fill(pt); //AllCh 
+
+             if(TMath::Abs(aodTrack->GetPdgCode()) == 211)
+               {
+                 _etadis->Fill(pt); //pion
+               } 
+             if(TMath::Abs(aodTrack->GetPdgCode()) == 321)
+               {
+                 _phidis->Fill(pt); //kaon
+               } 
+
+             //_etadis->Fill(eta);
+              //_phidis->Fill(phi); 
 
              if(fExcludeResonancesInMC)
                 {
-                 if (aodTrack->IsSecondaryFromWeakDecay()) continue;
-               }
+                  //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;                               
+                 Int_t gMotherIndex = aodTrack->GetMother();
+                  if(gMotherIndex != -1) {
+                    AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+                    if(motherTrack) {
+                      Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+                     
+                      if(pdgCodeOfMother == 311  ||
+                         pdgCodeOfMother == -311 ||
+                         pdgCodeOfMother == 310  ||
+                         pdgCodeOfMother == 3122 ||
+                         pdgCodeOfMother == -3122 ||
+                        pdgCodeOfMother == 111 ||
+                        pdgCodeOfMother == 22 ) continue;
+                    }
+                  }
+                }
              //Exclude electrons with PDG                                                                                   
                                                                                                                                
               if(fExcludeElectronsInMC) {
                 if(TMath::Abs(aodTrack->GetPdgCode()) == 11) continue;
               }
 
-             _etadis->Fill(eta);
-              _phidis->Fill(phi);
              //Particle loop 1st particle
               if (q > 0)
                 {
@@ -1349,6 +1361,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                     {
 
                       __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;                                            
+                      
                       __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;                                               
 
                     }
@@ -1469,24 +1482,22 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
            AliError("ERROR: Could not retrieve MC event");
          }
 
-         TExMap *trackMap = new TExMap();//Mapping matrix----
-          //1st loop track for Global tracks                                                              
-                                                                                                           
-          for(Int_t i = 0; i < fAODEvent->GetNumberOfTracks(); i++)
-            {
-              AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
-              if(!aodTrack) {
-                AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
-                continue;
-              }
-              Int_t gID = aodTrack->GetID();
-              if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                       
-                                                                                                           
-            }
-
-          AliAODTrack* newAodTrack;
-
-         for (int iTrack=0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++)
+         TExMap *trackMap = new TExMap();//Mapping matrix---- 
+         //1st loop track for Global tracks                                                                                                    
+         for(Int_t i = 0; i < fAODEvent->GetNumberOfTracks(); i++)
+           {
+             AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
+             if(!aodTrack) {
+               AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
+               continue;
+             }
+             Int_t gID = aodTrack->GetID();
+             if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                                                             
+           }
+         
+         AliAODTrack* newAodTrack;
+         
+         for (int iTrack=0; iTrack< fAODEvent->GetNumberOfTracks(); iTrack++)
            {
              
              AliAODTrack *t = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(iTrack));
@@ -1498,6 +1509,8 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              
              bitOK  = t->TestFilterBit(_trackFilterBit);
              if (!bitOK) continue;
+             Int_t gID = t->GetID();
+             newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
              
              q      = t->Charge();
              charge = int(q);
@@ -1510,43 +1523,92 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
              //Float_t dcaXY = t->DCA();  
              //Float_t dcaZ  = t->ZAtDCA(); 
              
-             Int_t gID = t->GetID();
-              newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
-
+             // get the electron nsigma                                                                                                
+             //Double_t nSigma = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
+             Double_t nSigmaPions   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kPion));
+              Double_t nSigmaKaons   = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kKaon));
+              //Double_t nSigmaProtons = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kProton));
+             
+             //Make the decision based on the n-sigma of electrons exclusively 
+             /* if(nSigma < fNSigmaCut
+                 && nSigmaPions   > fNSigmaCut
+                 && nSigmaKaons   > fNSigmaCut
+                 && nSigmaProtons > fNSigmaCut ) continue;
+             */
+
+             if(q == 0) continue;            
              // Kinematics cuts                                                                                 
              if( pt < 0.2 || pt > 2.0)      continue;
              if( eta < _min_eta_1 || eta > _max_eta_1)  continue;
              
-             //W/Wo Secondaries
-             //if (!AODmcTrack->IsPhysicalPrimary()) continue;
+             _dcaz->Fill(pt);
+          
+             if (nSigmaPions < fNSigmaCut)
+               {
+                 _etadis->Fill(pt);
+               }
+             if (nSigmaKaons < fNSigmaCut)
+               {
+                 _phidis->Fill(pt);
+               }
 
-             
+             /*
              Double_t pos[3];
-              newAodTrack->GetXYZ(pos);
-
-              Double_t DCAX = pos[0] - vertexX;
-              Double_t DCAY = pos[1] - vertexY;
-              Double_t DCAZ = pos[2] - vertexZ;
+             newAodTrack->GetXYZ(pos);
+             
+             Double_t DCAX = pos[0] - vertexX;
+             Double_t DCAY = pos[1] - vertexY;
+             Double_t DCAZ = pos[2] - vertexZ;
+             
+             Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+             
+             if (DCAZ     <  _dcaZMin ||
+                 DCAZ     >  _dcaZMax ||
+                 DCAXY    >  _dcaXYMax ) continue;
+             */
 
-              Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
 
-             Int_t label = TMath::Abs(t->GetLabel());
-             AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
-             if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;       
+      
+             //==== QA ===========================                                          
+             //_dcaz->Fill(DCAZ);                                                           
+             //_dcaxy->Fill(DCAXY);                                                         
+             //_etadis->Fill(eta);                                                          
+             //_phidis->Fill(phi); 
+             //===================================   
+             
+             //W/Wo Secondaries
+             //if (!AODmcTrack->IsPhysicalPrimary()) continue;
              
+             //cout<<"***************Prabhat on Weak Decay Particles ************"<<endl;
              if(fExcludeResonancesInMC)
                {
-                 if (AODmcTrack->IsSecondaryFromWeakDecay()) continue;
-                 //if (AODmcTrack->IsSecondaryFromMaterial()) continue;
+                 Int_t label = TMath::Abs(t->GetLabel());
+                 AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+                 
+                 Int_t gMotherIndex = AODmcTrack->GetMother();
+                 if(gMotherIndex != -1) {
+                   AliAODMCParticle* motherTrack = dynamic_cast<AliAODMCParticle *>(mcEvent->GetTrack(gMotherIndex));
+                   if(motherTrack) {
+                     Int_t pdgCodeOfMother = motherTrack->GetPdgCode();
+                     
+                     if(pdgCodeOfMother == 311  ||
+                        pdgCodeOfMother == -311 ||
+                        pdgCodeOfMother == 310  ||
+                        pdgCodeOfMother == 3122 ||
+                        pdgCodeOfMother == -3122 ||
+                        pdgCodeOfMother == 111 ||
+                        pdgCodeOfMother == 22 ) continue;
+                   }
+                 }
                }
-                     
-             //==== QA ===========================                                          
-             _dcaz->Fill(DCAZ);                                                           
-             _dcaxy->Fill(DCAXY);
-             _etadis->Fill(eta);                                                          
-             _phidis->Fill(phi); 
-             //===================================   
-
+             
+             Int_t label = TMath::Abs(t->GetLabel());
+             AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label);
+             if (AODmcTrack)
+               {
+                 if(TMath::Abs(AODmcTrack->GetPdgCode()) == 11) continue;
+               }
+             
              //Particle 1                                                                                                                  
              if (t->Charge() > 0)
                {
@@ -1589,6 +1651,7 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
                    {
 
                      __n1_1_vsPt[iPt]               += corr;          //cout << "step 15" << endl;                                           
+
                      __n1_1_vsZEtaPhiPt[iZEtaPhiPt] += corr;       //cout << "step 12" << endl;                                              
 
                    }
@@ -1702,13 +1765,13 @@ void  AliDptDptInMC::UserExec(Option_t */*option*/)
     } //AOD events             
 
 
-  _m0->Fill(_mult0);
-  _m1->Fill(_mult1);
-  _m2->Fill(_mult2);
-  _m3->Fill(_mult3);
-  _m4->Fill(_mult4);
-  _m5->Fill(_mult5);
-  _m6->Fill(_mult6);
+  //  _m0->Fill(_mult0);
+  //_m1->Fill(_mult1);
+  //_m2->Fill(_mult2);
+  //_m3->Fill(_mult3);
+  //_m4->Fill(_mult4);
+  //_m5->Fill(_mult5);
+  //_m6->Fill(_mult6);
   _vertexZ->Fill(vertexZ);
 
   if (_singlesOnly)
index 9154a56..39c3816 100644 (file)
@@ -6,8 +6,7 @@
 #include "AliLog.h"
 
 #include "AliPID.h"
-//#include "AliPIDResponse.h"
-
+#include "AliPIDResponse.h"
 
 class AliAODEvent;
 class AliESDEvent;
@@ -139,7 +138,7 @@ protected:
   AliESDEvent*             fESDEvent;             //! ESD Event 
   AliInputEventHandler*    fInputHandler;    //! Generic InputEventHandler 
   
-  //AliPIDResponse*          fPIDResponse;
+  AliPIDResponse*          fPIDResponse;
 
   // Histogram settings
   //TList*              _inputHistoList;
@@ -280,6 +279,7 @@ protected:
   double __s2PtNNw_12;
   
   double * __n1_1_vsPt;   //!
+
   double * __n1_1_vsEtaPhi;     //! 
   double * __s1pt_1_vsEtaPhi;    //!
   float  * __n1_1_vsZEtaPhiPt;    //!
index fb97415..676a0aa 100644 (file)
@@ -11,8 +11,8 @@
 /////////////////////////////////////////////////////////////////////////////////
 
 AliDptDptInMC *AddTaskCorrMC
-(int    singlesOnly            = 0,
- int    useWeights             = 1,
+(int    singlesOnly            = 1,
+ int    useWeights             = 0,
  int    centralityMethod       = 4,
  int    chargeSet              = 1,
  int    trackFilterBit         = 128,
@@ -22,8 +22,8 @@ AliDptDptInMC *AddTaskCorrMC
  double dcaZMax                =  3.0,
  double dcaXYMin               = -2.4,
  double dcaXYMax               =  2.4,
- int nCentrality               =  3, 
- TString anadata               = "MCAOD",
+ int nCentrality               =  5, 
+ TString anadata               = "MCAODreco",
  Bool_t NoResonances           =  kTRUE,
  Bool_t NoElectron             =  kTRUE,
  const char* taskname          = "MCTruth",
@@ -48,8 +48,10 @@ AliDptDptInMC *AddTaskCorrMC
        {
          
          minCentrality[0] = 0.0; maxCentrality[0] = 10.0; 
-         minCentrality[1] = 30.0; maxCentrality[1] = 40.0;
-         minCentrality[2] = 60.0; maxCentrality[2] = 70.0;
+         minCentrality[1] = 10.0; maxCentrality[1] = 20.0; 
+         minCentrality[2] = 20.0; maxCentrality[2] = 30.0; 
+         minCentrality[3] = 30.0; maxCentrality[3] = 40.0;
+         minCentrality[4] = 60.0; maxCentrality[4] = 70.0;
          
        }
       else