]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/Correlations/DPhi/AliAnalysisTaskDptDptCorrelations.cxx
Completed changes needed because of previous commit
[u/mrichter/AliRoot.git] / PWGCF / Correlations / DPhi / AliAnalysisTaskDptDptCorrelations.cxx
index 2f94b4150958e2c0a51cdf5167cd65e13295fef6..c03c7bad742daa346605297e802453a6c81bc05c 100644 (file)
@@ -22,6 +22,7 @@
 #include "TH3D.h"
 #include "THnSparse.h"
 #include "TCanvas.h"
+#include "TRandom.h"
 
 #include <TROOT.h>
 #include <TChain.h>
@@ -36,7 +37,6 @@
 #include <TH1D.h>
 #include <TH2D.h>
 #include <TH3D.h>
-#include <TRandom.h>
 #include "AliAnalysisManager.h"
 
 #include "AliAODHandler.h"
@@ -111,6 +111,8 @@ _dedxMax              ( 100000),
 _nClusterMin          ( 80), 
 _trackFilterBit       (0),
 fNSigmaCut            (3.),
+_tpcnclus             ( 50),
+_chi2ndf              (5.),
 _field    ( 1.),
 _nTracks  ( 0 ),
 _mult0    ( 0 ),
@@ -121,7 +123,7 @@ _mult4    ( 0 ),
 _mult4a    ( 0 ),
 _mult5    ( 0 ),
 _mult6    ( 0 ),
-arraySize ( 2000),
+arraySize ( 2500),
 _id_1(0),       
 _charge_1(0),    
 _iEtaPhi_1(0),    
@@ -394,6 +396,8 @@ _dedxMax              ( 100000),
 _nClusterMin          ( 80), 
 _trackFilterBit       ( 0),
 fNSigmaCut            ( 3.),
+_tpcnclus             ( 50),
+_chi2ndf              (5.),
 _field    ( 1.),
 _nTracks  ( 0 ),
 _mult0    ( 0 ),
@@ -404,7 +408,7 @@ _mult4    ( 0 ),
 _mult4a    ( 0 ),
 _mult5    ( 0 ),
 _mult6    ( 0 ),
-arraySize ( 2000),
+arraySize ( 2500),
 _id_1(0),       
 _charge_1(0),    
 _iEtaPhi_1(0),    
@@ -693,7 +697,7 @@ void AliAnalysisTaskDptDptCorrelations::UserCreateOutputObjects()
   _py_1       = new float[arraySize];   
   _pz_1       = new float[arraySize];   
   _correction_1 = new float[arraySize];
-    
+  _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.);
@@ -713,7 +717,7 @@ void AliAnalysisTaskDptDptCorrelations::UserCreateOutputObjects()
     _py_2       = new float[arraySize];   
     _pz_2       = new float[arraySize];   
     _correction_2 = new float[arraySize];
-        
+    _dedx_2       = new float[arraySize];
     __n1_2_vsPt              = getDoubleArray(_nBins_pt_2,        0.);
     __n1_2_vsEtaPhi          = getDoubleArray(_nBins_etaPhi_2,    0.);
     __s1pt_2_vsEtaPhi        = getDoubleArray(_nBins_etaPhi_2,    0.);
@@ -896,7 +900,7 @@ void  AliAnalysisTaskDptDptCorrelations::createHistograms()
   name = "m4"; _m4      = createHisto1D(name,name,_nBins_M4, _min_M4, _max_M4, _title_m4, _title_counts);
   name = "m5"; _m5      = createHisto1D(name,name,_nBins_M5, _min_M5, _max_M5, _title_m5, _title_counts);
   name = "m6"; _m6      = createHisto1D(name,name,_nBins_M6, _min_M6, _max_M6, _title_m6, _title_counts);
-  name = "zV"; _vertexZ = createHisto1D(name,name,1200, -12, 12, "z-Vertex (cm)", _title_counts);
+  name = "zV"; _vertexZ = createHisto1D(name,name,100, -10, 10, "z-Vertex (cm)", _title_counts);
   
 
   name = "Eta";     _etadis   = createHisto1F(name,name, 200, -1.0, 1.0, "#eta","counts");
@@ -904,7 +908,7 @@ void  AliAnalysisTaskDptDptCorrelations::createHistograms()
   name = "DCAz";    _dcaz     = createHisto1F(name,name, 500, -5.0, 5.0, "dcaZ","counts");
   name = "DCAxy";   _dcaxy    = createHisto1F(name,name, 500, -5.0, 5.0, "dcaXY","counts");
 
-  //name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
+  name = "Nclus1";   _Ncluster1    = createHisto1F(name,name, 200, 0, 200, "Ncluster1","counts");
   //name = "Nclus2";   _Ncluster2    = createHisto1F(name,name, 200, 0, 200, "Ncluster2","counts");
   
   if (_singlesOnly)
@@ -948,7 +952,7 @@ void  AliAnalysisTaskDptDptCorrelations::createHistograms()
     name = s2PtNNwName+pair_12_Name + vsM;    _s2PtNNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgSumPtN_12);       
     name = s2NPtNwName+pair_12_Name + vsM;    _s2NPtNw_12_vsM       = createProfile(name,name, _nBins_M4, _min_M4, _max_M4, _title_m4, _title_AvgNSumPt_12);     
     
-    name = "mInv";     _invMass     = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts"); 
+    name = "mInv";     _invMass     = createHisto1F(name,name, 50, 0.41, 0.55, "M_{inv}","counts"); 
     name = "mInvElec"; _invMassElec = createHisto1F(name,name, 500, 0., 1.000, "M_{inv}","counts"); 
     }
   
@@ -1011,19 +1015,20 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
   
   int    k1,k2;
   int    iPhi, iEta, iEtaPhi, iPt, charge;
-  float  q, phi, pt, eta, corr, corrPt, px, py, pz;
+  float  q, phi, pt, eta, corr, corrPt, px, py, pz, dedx;
   int    ij;
   int    id_1, q_1, iEtaPhi_1, iPt_1;
-  float  pt_1, px_1, py_1, pz_1, corr_1;
+  float  pt_1, px_1, py_1, pz_1, corr_1, dedx_1;
   int    id_2, q_2, iEtaPhi_2, iPt_2;
-  float  pt_2, px_2, py_2, pz_2, corr_2;
+  float  pt_2, px_2, py_2, pz_2, corr_2, dedx_2;
   float  ptpt;
   int    iVertex, iVertexP1, iVertexP2;
   int    iZEtaPhiPt;
-  float  massElecSq = 2.5e-7;
+  float  massElecSq = 1.94797849000000016e-02;
   //double b[2];
   //double bCov[3];
   const  AliAODVertex* vertex;
+  int    nClus;
   bool   bitOK;
   
   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
@@ -1074,15 +1079,15 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
   float vertexX  = -999;
   float vertexY  = -999;
   float vertexZ  = -999;
-  float vertexXY = -999;
-  float dcaZ     = -999;
-  float dcaXY    = -999;
+  //float vertexXY = -999;
+  //float dcaZ     = -999;
+  //float dcaXY    = -999;
   float centrality = -999;
   
   if(fAODEvent)
     {
       //Centrality
-      AliCentrality* centralityObject =  fAODEvent->GetHeader()->GetCentralityP();
+      AliCentrality* centralityObject =  ((AliVAODHeader*)fAODEvent->GetHeader())->GetCentralityP();
       if (centralityObject)
        {
          //cout << "AliAnalysisTaskDptDptCorrelations::UserExec(Option_t *option) - 6" << endl;
@@ -1124,26 +1129,29 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
       _eventAccounting->Fill(2);// count all events with right centrality
       
       // filter on z and xy vertex
-      vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertexSPD();
-      if (!vertex || vertex->GetNContributors()<1)
+      vertex = (AliAODVertex*) fAODEvent->GetPrimaryVertex();
+      // Double_t V[2];
+      //vertex->GetXYZ(V);      
+
+      if(vertex)
        {
-         vertexZ  = -999;
-         vertexXY = -999;
-         AliInfo("AliAnalysisTaskDptDptCorrelations::UserExec(Option_t *option) - No valid vertex object or poor vertex");
-       }
-      else
-       { 
-         vertexX = vertex->GetX();
-         vertexY = vertex->GetY();
-         vertexZ = vertex->GetZ();
-         vertexXY = sqrt(vertexX*vertexX+vertexY*vertexY);
+         Double32_t fCov[6];
+         vertex->GetCovarianceMatrix(fCov);
+         if(vertex->GetNContributors() > 0)
+           {
+             if(fCov[5] != 0)
+               {
+                 vertexX = vertex->GetX();
+                 vertexY = vertex->GetY();
+                 vertexZ = vertex->GetZ();
+                 
+                 if(TMath::Abs(vertexZ) > 10)
+                   {
+                     return;
+                   } // Z-Vertex Cut  
+               }
+           }
        }
-      if (!vertex ||
-         vertexZ  < _vertexZMin  || 
-         vertexZ  > _vertexZMax  ||
-         vertexXY < _vertexXYMin || 
-         vertexXY > _vertexXYMax)
-       return;
       
       _vertexZ->Fill(vertexZ);
       
@@ -1159,8 +1167,9 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
       //====================== 
       
       //*********************************************************
-      TExMap *trackMap = new TExMap();//Mapping matrix----                                            
-      //1st loop track                                                                                
+       TExMap *trackMap = new TExMap();//Mapping matrix----                                            
+
+      //1st loop track for Global tracks                                                                                
       for(Int_t i = 0; i < _nTracks; i++)
        {
          AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAODEvent->GetTrack(i));
@@ -1170,7 +1179,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
          }
          Int_t gID = aodTrack->GetID();
          if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks                       
-       }
+         }
            
       AliAODTrack* newAodTrack;
       
@@ -1184,11 +1193,12 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
          }
          
          bitOK  = t->TestFilterBit(_trackFilterBit);
-         if (!bitOK) continue; //128bit
-                 
-         Int_t gID = t->GetID();
-         newAodTrack = gID >= 0 ?t : fAODEvent->GetTrack(trackMap->GetValue(-1-gID));
+         if (!bitOK) continue; //128bit or 272bit
          
+         Int_t gID = t->GetID();
+         newAodTrack = gID >= 0 ?t : dynamic_cast<AliAODTrack*>(fAODEvent->GetTrack(trackMap->GetValue(-1-gID)));
+         if(!newAodTrack) AliFatal("Not a standard AOD?");
          q      = t->Charge();
          charge = int(q);
          phi    = t->Phi();
@@ -1197,10 +1207,23 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
          py     = t->Py();
          pz     = t->Pz();
          eta    = t->Eta();
-         dcaXY = t->DCA(); 
-         dcaZ  = t->ZAtDCA();  
+         dedx   = t->GetTPCsignal();
+         //dcaXY = t->DCA(); 
+         //dcaZ  = t->ZAtDCA();  
+         nClus  = t->GetTPCNcls();       
+         
+          if ( nClus<_nClusterMin ) continue;
          
-         Double_t nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
+         _Ncluster1->Fill(nClus);
+         
+         /*
+         //cuts on more than 0 shared cluster (suggested by Michael)
+         if(t->GetTPCnclsS() > 0){
+         continue;
+         }*/
+         
+         //for Global tracks
+          Double_t nsigmaelectron = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kElectron));
          Double_t nsigmapion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kPion));
          Double_t nsigmakaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kKaon));
          Double_t nsigmaproton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,(AliPID::EParticleType)AliPID::kProton));
@@ -1212,55 +1235,45 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
             && nsigmakaon   > fNSigmaCut
             && nsigmaproton > fNSigmaCut ) continue;
          
-         if(charge == 0) continue;
-         
-         /*if (newAodTrack->PropagateToDCA(vertex, _field, 100., b, bCov))
-           {
-             AliAODTrack* clone = (AliAODTrack*)newAodTrack->Clone("trk_clone");
-             if (clone->PropagateToDCA(vertex, _field, 100., b, bCov))
-               {
-                 dcaXY = b[0];
-                 dcaZ = b[1];
-               }
-             else{
-               dcaXY = -999999;
-               dcaZ = -999999;
-             }
-             delete clone;
-             }*/
-
 
+         if(charge == 0) continue;
          // Kinematics cuts used                                                                                        
-         if( pt < _min_pt_1 || pt > _max_pt_1)      continue;
-         if( eta < _min_eta_1 || eta > _max_eta_1)  continue;
+         if( pt < _min_pt_1 || pt > _max_pt_1) continue;
+         if( eta < _min_eta_1 || eta > _max_eta_1) continue;
          
+         /*      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;
+               
+         Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
+         
+         if (DCAZ     <  _dcaZMin || 
+             DCAZ     >  _dcaZMax ||
+             DCAXY    >  _dcaXYMax ) continue; 
+         */
+
+           //------- Eff. test---------- //just for checking
+         //Double_t yy = (1 - 0.7)/1.8;
+         //Double_t zz = (pt - 0.2);
+         //Double_t effValue = 0.7 + yy*zz;
+         //Double_t R = gRandom->Rndm();
+          //if(R > effValue) continue;
+         //---------------------------   
          
-         /*Double_t pos[3];
-
-         // if not constrained track the DCA is stored as position (at the vertex)
-         if(!t->GetXYZ(pos)){
-           dcaXY = TMath::Sign(1.0,pos[0])*TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
-           dcaZ  = pos[2];
-           }*/
-
-                 
-         if (dcaZ     <  _dcaZMin ||
-             dcaZ     >  _dcaZMax ||
-             dcaXY    <  _dcaXYMin ||
-             dcaXY    >  _dcaXYMax
-             )continue; 
          //==== QA ===========================
-         _dcaz->Fill(dcaZ);
-         _dcaxy->Fill(dcaXY);
+         //_dcaz->Fill(DCAZ);
+         //_dcaxy->Fill(DCAXY);
          _etadis->Fill(eta);
-         _phidis->Fill(phi);
+         //_phidis->Fill(phi);
          //===================================
          //*************************************************
                  
          //Particle 1
-         if (_requestedCharge_1 == charge
+         if (_requestedCharge_1 == charge && dedx >=  _dedxMin && dedx < _dedxMax)
            {
-             
              iPhi   = int( phi/_width_phi_1);
              
              if (iPhi<0 || iPhi>=_nBins_phi_1 ) 
@@ -1329,7 +1342,9 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
                }
            }
          
-         if (!_sameFilter && _requestedCharge_2 == charge)  
+         if (!_sameFilter && _requestedCharge_2 == charge && 
+             dedx >=  _dedxMin && dedx < _dedxMax)
+              
            {
              
              iPhi   = int( phi/_width_phi_2);
@@ -1444,7 +1459,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
              iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
              corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
              pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
-             //dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
+             dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
              //1 and 2
              for (int i2=i1+1; i2<k1; i2++)
                {        
@@ -1457,7 +1472,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
                      iPt_2     = _iPt_1[i2];        ////cout << "      iPt_1:" << iPt_1 << endl;
                      corr_2    = _correction_1[i2]; ////cout << "     corr_1:" << corr_1 << endl;
                      pt_2      = _pt_1[i2];         ////cout << "       pt_1:" << pt_1 << endl;
-                     //dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
+                     dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
                      corr      = corr_1*corr_2;
                      if (q_2>q_1 || (q_1>0 && q_2>0 && pt_2<=pt_1) || (q_1<0 && q_2<0 && pt_2>=pt_1))
                        {
@@ -1499,7 +1514,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
              iPt_1     = _iPt_1[i1];          ////cout << "      iPt_1:" << iPt_1 << endl;
              corr_1    = _correction_1[i1];   ////cout << "     corr_1:" << corr_1 << endl;
              pt_1      = _pt_1[i1];           ////cout << "       pt_1:" << pt_1 << endl;
-             //dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
+             dedx_1    = _dedx_1[i1];         ////cout << "     dedx_1:" << dedx_1 << endl;
              //1 and 2
              for (int i2=i1+1; i2<k1; i2++)
                {        
@@ -1512,7 +1527,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
                      iPt_2     = _iPt_1[i2];        ////cout << "      iPt_2:" << iPt_2 << endl;
                      corr_2    = _correction_1[i2]; ////cout << "     corr_2:" << corr_2 << endl;
                      pt_2      = _pt_1[i2];         ////cout << "       pt_2:" << pt_2 << endl;
-                     //dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
+                     dedx_2    = _dedx_1[i2];       ////cout << "     dedx_2:" << dedx_2 << endl;
                      corr      = corr_1*corr_2;
                      if ( q_2<q_1 || (q_1>0 && q_2>0 && pt_2>=pt_1) || (q_1<0 && q_2<0 && pt_2<=pt_1))
                        {
@@ -1569,7 +1584,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
              px_1      = _px_1[i1];          ////cout << "      px_1:" << px_1 << endl;
              py_1      = _py_1[i1];          ////cout << "      py_1:" << py_1 << endl;
              pz_1      = _pz_1[i1];          ////cout << "      pz_1:" << pz_1 << endl;
-             //dedx_1    = _dedx_1[i1];        ////cout << "     dedx_1:" << dedx_1 << endl;
+             dedx_1    = _dedx_1[i1];        ////cout << "     dedx_1:" << dedx_1 << endl;
              
              //1 and 2
              for (int i2=0; i2<k2; i2++)
@@ -1586,7 +1601,7 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
                      px_2      = _px_2[i2];          ////cout << "      px_2:" << px_2 << endl;
                      py_2      = _py_2[i2];          ////cout << "      py_2:" << py_2 << endl;
                      pz_2      = _pz_2[i2];          ////cout << "      pz_2:" << pz_2 << endl;
-                     //dedx_2    = _dedx_2[i2];        ////cout << "     dedx_2:" << dedx_2 << endl;
+                     dedx_2    = _dedx_2[i2];        ////cout << "     dedx_2:" << dedx_2 << endl;
                      
                      
                      if (_rejectPairConversion)
@@ -1596,6 +1611,14 @@ void  AliAnalysisTaskDptDptCorrelations::UserExec(Option_t */*option*/)
                          float mInvSq = 2*(massElecSq + sqrt(e1Sq*e2Sq) - px_1*px_2 - py_1*py_2 - pz_1*pz_2 );
                          float mInv = sqrt(mInvSq);
                          _invMass->Fill(mInv);
+                         if (mInv<0.51)
+                           {
+                             if (dedx_1>75. && dedx_2>75.)
+                               {
+                                 //_invMassElec->Fill(mInv);
+                                 if (mInv<0.05) continue;
+                               }
+                           }
                        }
                      
                      corr      = corr_1*corr_2;