Analysis QA task exploiting symmetries (charge, phi, eta+/-, ...)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 5 Dec 2009 03:09:22 +0000 (03:09 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 5 Dec 2009 03:09:22 +0000 (03:09 +0000)
E. Sicking, MvL, AM ..

PWG1/AliAnalysisTaskQASym.cxx [new file with mode: 0644]
PWG1/AliAnalysisTaskQASym.h [new file with mode: 0644]

diff --git a/PWG1/AliAnalysisTaskQASym.cxx b/PWG1/AliAnalysisTaskQASym.cxx
new file mode 100644 (file)
index 0000000..f44e913
--- /dev/null
@@ -0,0 +1,932 @@
+#include "TChain.h"
+#include "TTree.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "TList.h"
+#include "TParticle.h"
+#include "TParticlePDG.h"
+#include "TProfile.h"
+#include "TNtuple.h"
+#include "TFile.h"
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+
+#include "AliESDEvent.h"
+
+#include "AliESDVertex.h"
+#include "AliESDInputHandler.h"
+#include "AliESDtrackCuts.h"
+#include "AliMultiplicity.h"
+
+
+#include "AliAnalysisTaskQASym.h"
+#include "AliExternalTrackParam.h"
+#include "AliTrackReference.h"
+
+#include "AliHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliGenDPMjetEventHeader.h"
+
+// Analysis Task for basic QA on the ESD
+
+// Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing,
+//          Andreas Morsch, Eva Sicking
+
+ClassImp(AliAnalysisTaskQASym)
+
+//________________________________________________________________________
+AliAnalysisTaskQASym::AliAnalysisTaskQASym(const char *name) 
+  : AliAnalysisTaskSE(name) 
+    ,fFieldOn(kTRUE)
+    ,fHists(0)
+    ,fHistRECpt(0)
+    ,fEta(0)
+    ,fEtaPt(0)
+    ,fQPt(0)
+    ,fDca(0)
+    ,fqRec(0)
+    ,fsigmaPt(0)
+  
+    ,fRecPtPos(0)
+    ,fRecPtNeg(0)
+    ,fRecPhiPos(0)
+    ,fRecPhiNeg(0)
+    ,fRecEtaPos(0)
+    ,fRecEtaNeg(0)
+    ,fRecEtaPtPos(0)
+    ,fRecEtaPtNeg(0)
+    ,fRecDcaPos(0)
+    ,fRecDcaNeg(0)
+    ,fRecDcaNegInv(0)
+    ,fRecDPos(0)
+    ,fRecDNeg(0)
+
+    ,fRecQPtPosEta(0)
+    ,fRecQPtNegEta(0)
+    ,fRecPtPosEta(0)
+    ,fRecPtNegEta(0)
+    ,fRecPhiPosEta(0)
+    ,fRecPhiNegEta(0)
+    ,fRecDcaPosEta(0)
+    ,fRecDcaNegEta(0)
+    ,fRecDPosEta(0)
+    ,fRecDNegEta(0)
+
+    ,fRecPtPosVz(0)
+    ,fRecPtNegVz(0)
+    ,fRecEtaPosVz(0)
+    ,fRecEtaNegVz(0)
+    ,fRecPhiPosVz(0)
+    ,fRecPhiNegVz(0)
+    ,fSignedDcaPosVz(0)
+    ,fSignedDcaNegVz(0)
+    ,fRecQPtPosEtaVz(0)
+    ,fRecQPtNegEtaVz(0)
+    ,fRecEtaPtPosVz(0)
+    ,fRecEtaPtNegVz(0)
+
+
+    ,fDeltaPhiAll(0)
+    ,fDeltaPhiLeading(0) 
+    ,fDiffDcaD(0)
+    ,fPhiRec(0)
+    ,fThetaRec(0)
+    ,fNumber(0)
+    ,fVx(0)
+    ,fVy(0)
+    ,fVz(0)
+    ,test(0)
+
+
+    ,fCuts(0)
+
+   
+  
+{
+    // Constructor
+  for(Int_t i = 0;i<18;++i){
+    fRecPtTpcSector[i] = 0;
+    fRecEtaTpcSector[i] = 0;
+    fSignedDcaTpcSector[i] = 0;
+    fRecQPtTpcSector[i] = 0;
+    fRecEtaPtTpcSector[i] = 0;
+  }
+
+  for(Int_t i = 0;i< 7;++i){
+    fRecPtPosLadder[i] = 0;
+    fRecPtNegLadder[i] = 0;
+    fRecPhiPosLadder[i] = 0;
+    fRecPhiNegLadder[i] = 0;
+    fRecEtaPosLadder[i] = 0;
+    fRecEtaNegLadder[i] = 0;
+    fSignDcaPos[i] = 0;
+    fSignDcaNeg[i] = 0;
+    fSignDcaNegInv[i] = 0;
+    fPtSigmaPos[i] =0;
+    fPtSigmaNeg[i] =0;
+    fqPtRec[i] =0;
+    fDcaSigmaPos[i] =0;
+    fDcaSigmaNeg[i] =0;
+  }
+
+  DefineOutput(1,  TList::Class()); 
+
+  
+  
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskQASym::UserCreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+
+  Double_t range = 300.;
+  Double_t pt = 20.;
+
+  fHists = new TList();
+  test   = new TNtuple("test","test",  
+                         "pt:phi:theta:x:y:z:charge");
+  fHistRECpt   = new TH1F("fHistRECpt", 
+                         " p_{T}",
+                         100, 0., pt);
+  fEta   = new TH1F("fEta", 
+                         " #eta",
+                         200, -2., 2.);
+  fThetaRec   = new TH1F("fThetaRec", 
+                         " #theta",
+                         180, 0., TMath::Pi());
+  fPhiRec   = new TH1F("fPhiRec", 
+                         " #phi",
+                         180, 0., 2*TMath::Pi());
+  fNumber   = new TH1F("fNumber", 
+                      "number of tracks per event",
+                      50, 0.5, 49.5);
+  fVx   = new TH1F("fVx", 
+                  "X of vertex",
+                  500, -500., 500.);
+  fVy   = new TH1F("fVy", 
+                  "Y of vertex",
+                  500, -500., 500.);
+  fVz   = new TH1F("fVz", 
+                  "Z of vertex",
+                  500, -500., 500.);
+
+  fEtaPt   = new TH1F("fEtaPt", 
+                         " #eta/p_{T} ",
+                         100, -1., 1.);
+
+  fQPt   = new TH1F("fQPt", 
+                         " charge/p_{T} ",
+                         100, -1., 1.);
+
+  fDca   = new TH1F("fDca", 
+                         " dca ",
+                         200, -range, range);
+
+
+  fqRec    = new TH1F("fqRec",   
+                         " charge all reconstructed particle",
+                         21, -9.5, 10.5);
+  
+  fsigmaPt    = new TH1F("fsigmaPt",   
+                         "Log_{10}(#sigma_{p_{T}})",
+                         200, -2., 8.);
+
+
+
+
+  //------------
+  for(Int_t ITSlayer_case=0;ITSlayer_case<7;ITSlayer_case++){
+
+    fSignDcaPos[ITSlayer_case]   = new TH1F(Form("fSignDcaPos%d", ITSlayer_case),  
+                                  " Signed dca", 
+                                  200, -range, range);
+    fSignDcaPos[ITSlayer_case]->GetXaxis()->SetTitle("dca");
+    fSignDcaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
+   
+    fSignDcaNeg[ITSlayer_case]   = new TH1F(Form("fSignDcaNeg%d", ITSlayer_case),  
+                                  " Signed dcas",
+                                  200, -range, range);
+    fSignDcaNeg[ITSlayer_case]->GetXaxis()->SetTitle("dca");
+    fSignDcaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
+
+    fSignDcaNegInv[ITSlayer_case]   = new TH1F(Form("fSignDcaNegInv%d", ITSlayer_case),  
+                                  " inverse Signed dca ",
+                                  200, -range, range);
+    fSignDcaNegInv[ITSlayer_case]->GetXaxis()->SetTitle("-dca");
+    fSignDcaNegInv[ITSlayer_case]->GetYaxis()->SetTitle("");
+
+
+
+
+    fPtSigmaPos[ITSlayer_case]   = new TH1F(Form("fPtSigmaPos%d", ITSlayer_case),  
+                                  " #sigma_{pT} ",
+                                  208, -2., 8.);
+    fPtSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
+    fPtSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("");
+    
+    
+    fPtSigmaNeg[ITSlayer_case]   = new TH1F(Form("fPtSigmaNeg%d",ITSlayer_case),  
+                                 " #sigma_{pT}",
+                                  208, -2., 8.);
+    fPtSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("Log_{10}(#sigma_{pT})");
+    fPtSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("");
+
+
+
+
+
+    fqPtRec[ITSlayer_case]   = new TH1F(Form("fqPtRec%d",ITSlayer_case),  
+                                 "q/ p_{T}",
+                                  200, -100., 100.);
+    fqPtRec[ITSlayer_case]->GetXaxis()->SetTitle("q_{tr}/p_{T, tr} (GeV/c)");
+    fqPtRec[ITSlayer_case]->GetYaxis()->SetTitle("");
+
+  
+   
+
+
+    fDcaSigmaPos[ITSlayer_case]   = new TH2F(Form("fDcaSigmaPos%d", ITSlayer_case),  
+                                  " p_{T} shift vs #sigma_{pT} ",
+                                  200, -range, range,200, -4., 4. );
+    fDcaSigmaPos[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA)");
+    fDcaSigmaPos[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
+    
+    
+    fDcaSigmaNeg[ITSlayer_case]   = new TH2F(Form("fDcaSigmaNeg%d", ITSlayer_case),  
+                                  " p_{T} shift vs #sigma_{pT} ",
+                                  200, -range, range,200, -4., 4. );
+    fDcaSigmaNeg[ITSlayer_case]->GetXaxis()->SetTitle("signed DCA");
+    fDcaSigmaNeg[ITSlayer_case]->GetYaxis()->SetTitle("log_{10}(#sigma_{pT})");
+
+
+    
+    
+    
+    // YIELDs---------- positive and negative particles
+    
+    fRecPtPos   = new TH1F("fRecPtPos", 
+                          " p_{T}",
+                          100, 0.,pt);
+    fRecPtPos->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    fRecPtNeg   = new TH1F("fRecPtNeg", 
+                          " p_{T} ",
+                          100, 0., pt);
+    fRecPtNeg->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+    
+    fRecPhiPos   = new TH1F("fRecPhiPos", 
+                           "#phi",
+                           361, 0., 360.);
+    fRecPhiPos->GetXaxis()->SetTitle("#phi (deg)");
+  
+    fRecPhiNeg   = new TH1F("fRecPhiNeg", 
+                           "#phi ",
+                           361, 0., 360.);
+    fRecPhiNeg->GetXaxis()->SetTitle("#phi (deg)");
+    
+    fRecEtaPos   = new TH1F("fRecEtaPos", 
+                           "#eta",
+                           200, -2., 2.);
+    fRecEtaPos->GetXaxis()->SetTitle("#eta");
+
+    fRecEtaNeg   = new TH1F("fRecEtaNeg", 
+                           "#eta",
+                           200, -2., 2.);
+    fRecEtaNeg->GetXaxis()->SetTitle("#eta");
+    
+    fRecEtaPtPos   = new TH1F("fRecEtaPtPos", 
+                             "#eta/p_{T}",
+                             200, -0.1, .1);
+    fRecEtaPtPos->GetXaxis()->SetTitle("#eta/p_{T}");
+
+    fRecEtaPtNeg   = new TH1F("fRecEtaPtNeg", 
+                             "#eta/p_{T}",
+                             200, -.1, .1);
+    fRecEtaPtNeg->GetXaxis()->SetTitle("#eta/p_{T}");
+
+    fRecDcaPos   = new TH1F("fRecDcaPos", 
+                        " dca",
+                          100, -range, range);
+    fRecDcaPos->GetXaxis()->SetTitle("dca (cm)");
+    fRecDcaNeg   = new TH1F("fRecDcaNeg", 
+                          " dca",
+                          100, -range, range);
+    fRecDcaNeg->GetXaxis()->SetTitle("dca (cm)");
+
+    fRecDcaNegInv   = new TH1F("fRecDcaNegInv", 
+                          " dca",
+                          100, -range, range);
+    fRecDcaNegInv->GetXaxis()->SetTitle("dca (cm)");
+
+
+    fRecDPos   = new TH1F("fRecDPos", 
+                        " d",
+                          100, -range, range);
+    fRecDPos->GetXaxis()->SetTitle("d (cm)");
+    fRecDNeg   = new TH1F("fRecDNeg", 
+                          "d",
+                          100, -range, range);
+    fRecDNeg->GetXaxis()->SetTitle("d (cm)");
+
+
+    //  YIELDs ---------------- positive and negative eta
+    
+    
+    fRecQPtPosEta   = new TH1F("fRecQPtPosEta", 
+                              "q/p_{T}",
+                              200, -0.5, 0.5);
+    fRecQPtPosEta->GetXaxis()->SetTitle("q/p_{T} ");
+
+    fRecQPtNegEta   = new TH1F("fRecQPtNegEta", 
+                              "q/p_{T}",
+                              200, -0.5, 0.5);
+    fRecQPtNegEta->GetXaxis()->SetTitle("q/p_{T}");
+    
+    fRecPtPosEta   = new TH1F("fRecPtPosEta", 
+                             " p_{T} ",
+                             100, 0., pt);
+    fRecPtPosEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+    fRecPtNegEta   = new TH1F("fRecPtNegEta", 
+                             " p_{T}",
+                             100, 0., pt);
+    fRecPtNegEta->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+    
+    fRecPhiPosEta   = new TH1F("fRecPhiPosEta", 
+                           "#phi",
+                           361, 0., 360);
+    fRecPhiPosEta->GetXaxis()->SetTitle("#phi (deg)");
+
+    fRecPhiNegEta   = new TH1F("fRecPhiNegEta", 
+                           "#phi ",
+                           361, 0, 360);
+    fRecPhiNegEta->GetXaxis()->SetTitle("#phi (deg)");
+
+    fRecDcaPosEta   = new TH1F("fRecDcaPosEta", 
+                        " dca ",
+                          100, -range, range);
+    fRecDcaPosEta->GetXaxis()->SetTitle("dca (cm)");
+    fRecDcaNegEta   = new TH1F("fRecDcaNegEta", 
+                          " dca",
+                          100, -range, range);
+    fRecDcaNegEta->GetXaxis()->SetTitle("dca (cm)");
+
+    fRecDPosEta   = new TH1F("fRecDPosEta", 
+                        " d",
+                          100, -range, range);
+    fRecDPosEta->GetXaxis()->SetTitle("d (cm)");
+    fRecDNegEta   = new TH1F("fRecDNegEta", 
+                          "d",
+                          100, -5., 5.);
+    fRecDNegEta->GetXaxis()->SetTitle("d (cm)");
+
+
+    
+    //  YIELDs ---------------- for TPC sectors
+    for(Int_t sector=0; sector<18;sector++){
+      
+
+      fRecPtTpcSector[sector]   = new TH1F(Form("fRecPtTpcSector%02d",sector), 
+                                          Form("p_{T} distribution: TPC sector %d",
+                                               sector),100, 0., pt);
+      fRecPtTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+      fRecEtaTpcSector[sector]   = new TH1F(Form("fRecEtaTpcSector%02d",sector), 
+                                          Form("#eta distribution: TPC sector %d",
+                                               sector),200, -2., 2.);
+      fRecEtaTpcSector[sector]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+     
+
+      fSignedDcaTpcSector[sector]   = new TH1F(Form("fSignedDcaTpcSector%02d",sector), 
+                                          Form("dca distribution: TPC sector %d",
+                                               sector),200, -range, range );
+      fSignedDcaTpcSector[sector]->GetXaxis()->SetTitle("dca");
+
+      fRecQPtTpcSector[sector]   = new TH1F(Form("fRecQPtTpcSector%02d",sector), 
+                                          Form("Q/ p_{T} distribution: TPC sector %d",
+                                               sector),100, -1., 1.);
+      fRecQPtTpcSector[sector]->GetXaxis()->SetTitle("Q/p_{T} (GeV/c)");
+
+      fRecEtaPtTpcSector[sector]   = new TH1F(Form("fRecEtaPtTpcSector%02d",sector), 
+                                          Form("#eta/ p_{T} distribution: TPC sector %d",
+                                               sector),100, -1., 1.);
+      fRecEtaPtTpcSector[sector]->GetXaxis()->SetTitle("#eta/p_{T} (GeV/c)");
+    }
+    // YIELDS ITS ladder
+    for(Int_t i=0;i<7;i++){
+      fRecPtPosLadder[i]   = new TH1F(Form("fRecPtPosLadder%d", i), 
+                            " p_{T} distribution",
+                            100, 0., pt);
+      fRecPtPosLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+      fRecPtNegLadder[i]   = new TH1F(Form("fRecPtNegLadder%d",i), 
+                            " p_{T} distribution ",
+                            100, 0., pt);
+      fRecPtNegLadder[i]->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+
+
+      fRecPhiPosLadder[i]   = new TH1F(Form("fRecPhiPosLadder%d",i), 
+                                "#phi distribution: all pos eta",
+                                361, 0., 360);
+      fRecPhiPosLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
+      
+      fRecPhiNegLadder[i]   = new TH1F(Form("fRecPhiNegLadder%d", i),
+                                "#phi distribution: all neg eta",
+                                361, 0, 360);
+      fRecPhiNegLadder[i]->GetXaxis()->SetTitle("#phi (deg)");
+
+
+
+      fRecEtaPosLadder[i]   = new TH1F(Form("fRecEtaPosLadder%d",i), 
+                                      "#eta distribution",
+                                200, -2., 2.);
+      fRecEtaPosLadder[i]->GetXaxis()->SetTitle("#eta)");
+      
+      fRecEtaNegLadder[i]   = new TH1F(Form("fRecEtaNegLadder%d", i),
+                                "#eta distribution",
+                                200, -2., 2.);
+      fRecEtaNegLadder[i]->GetXaxis()->SetTitle("#eta");
+    }
+
+    Double_t vzmax = 30;
+
+    fRecPtPosVz = new TH2F("fRecPtPosVz", 
+                          "p_{T} distribution vs Vz()",
+                          100, -1., 2., 200,-vzmax,vzmax);
+    fRecPtPosVz->GetXaxis()->SetTitle("log_{10}(p_{T})");
+    
+    fRecPtNegVz = new TH2F("fRecPtNegVz",
+                          "p_{T} distribution vs Vz()",
+                          100, -1., 2.,200,-vzmax,vzmax);
+    fRecPtNegVz->GetXaxis()->SetTitle("Log_{10}(p_{T})");
+    
+   
+    fRecEtaPosVz= new TH2F("fRecEtaPosVz", 
+                         "#eta distribution vs Vz()",
+                         100, -2., 2., 200,-vzmax,vzmax);
+    fRecEtaPosVz->GetXaxis()->SetTitle("#eta");
+    fRecEtaNegVz = new TH2F("fRecEtaNegVz",
+                          "#eta distribution vs Vz()",
+                          100, -2., 2.,200,-vzmax,vzmax);
+    fRecEtaNegVz->GetXaxis()->SetTitle("#eta");
+
+    fRecPhiPosVz= new TH2F("fRecPhiPosVz", 
+                         "#eta distribution vs Vz()",
+                         361, 0., 360., 200,-vzmax,vzmax);
+    fRecPhiPosVz->GetXaxis()->SetTitle("#phi (deg)");
+    fRecPhiNegVz = new TH2F("fRecPhiNegVz",
+                          "dca vs Vz()",
+                          361, 0., 360.,200,-vzmax,vzmax);
+    fRecPhiNegVz->GetXaxis()->SetTitle("#phi (deg)");
+
+    fSignedDcaPosVz= new TH2F("fSignedDcaPosVz", 
+                         "#eta distribution vs Vz()",
+                         200, -range, range, 200,-vzmax,vzmax);
+    fSignedDcaPosVz->GetXaxis()->SetTitle("dca (cm)");
+    fSignedDcaNegVz = new TH2F("fSignedDcaNegVz",
+                          "dca vs Vz()",
+                          200, -range, range,200,-vzmax,vzmax);
+    fSignedDcaNegVz->GetXaxis()->SetTitle("dca (cm)");
+
+    fRecQPtPosEtaVz= new TH2F("fRecQPtPosEtaVz",
+                          " Q/p_{T} distribution vs Vz()",
+                          100, -1., 1., 200,-vzmax,vzmax);
+    fRecQPtPosEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
+    fRecQPtNegEtaVz = new TH2F("fRecQPtNegEtaVz",
+                           " Q/p_{T} distribution vs Vz()",
+                           100, -1., 1.,200,-vzmax,vzmax);
+    fRecQPtNegEtaVz->GetXaxis()->SetTitle("Q/p_{T}");
+
+    fRecEtaPtPosVz= new TH2F("fRecEtaPtPosVz",
+                          " #eta/p_{T} distribution vs Vz()",
+                          100, -1., 1., 200,-vzmax,vzmax);
+    fRecEtaPtPosVz->GetXaxis()->SetTitle("#eta/p_{T");
+    fRecEtaPtNegVz = new TH2F("fRecEtaPtNegVz",
+                           " #eta/p_{T} distribution vs Vz()",
+                           100, -1., 1.,200,-vzmax,vzmax);
+    fRecEtaPtNegVz->GetXaxis()->SetTitle("#eta/p_{T}");
+
+    //new
+    fDeltaPhiAll = new TH1F("fDeltaPhiAll",
+                           " #Delta #phi",200,-360,360);
+    fDeltaPhiAll->GetXaxis()->SetTitle("#Delta #phi");
+
+
+    fDeltaPhiLeading = new TH2F("fDeltaPhiLeading",
+                           " #Delta #phi",361,-360,360, 361,0, 360);
+    fDeltaPhiLeading->GetXaxis()->SetTitle("#Delta #phi (deg.)");
+    fDeltaPhiLeading->GetYaxis()->SetTitle("#phi_{leading particle} (deg.)");
+
+    fDiffDcaD    = new TH1F("fDiffDcaD",   
+                           "dca-d",
+                           200, -5., 5.);
+    
+  }
+
+  fHists->SetOwner();
+
+  fHists->Add(test);
+  fHists->Add(fHistRECpt);
+  fHists->Add(fEta);
+  fHists->Add(fThetaRec);
+  fHists->Add(fPhiRec);
+  fHists->Add(fNumber);
+  fHists->Add(fVx);
+  fHists->Add(fVy);
+  fHists->Add(fVz);
+
+  fHists->Add(fEtaPt);
+  fHists->Add(fQPt);
+  fHists->Add(fDca);
+
+  fHists->Add(fDeltaPhiAll);
+  fHists->Add(fDeltaPhiLeading);
+  fHists->Add(fDiffDcaD);
+
+  fHists->Add(fqRec);
+  fHists->Add(fsigmaPt);
+
+  fHists->Add(fRecPtPos);
+  fHists->Add(fRecPtNeg);
+  fHists->Add(fRecPhiPos);
+  fHists->Add(fRecPhiNeg);
+  fHists->Add(fRecEtaPos);
+  fHists->Add(fRecEtaNeg);
+  fHists->Add(fRecEtaPtPos);
+  fHists->Add(fRecEtaPtNeg);
+  fHists->Add(fRecDcaPos);
+  fHists->Add(fRecDcaNeg);
+  fHists->Add(fRecDcaNegInv);
+  fHists->Add(fRecDPos);
+  fHists->Add(fRecDNeg);
+
+
+  fHists->Add(fRecQPtPosEta);
+  fHists->Add(fRecQPtNegEta);
+  fHists->Add(fRecPtPosEta);
+  fHists->Add(fRecPtNegEta);
+  fHists->Add(fRecPhiPosEta);
+  fHists->Add(fRecPhiNegEta);
+  fHists->Add(fRecDcaPosEta);
+  fHists->Add(fRecDcaNegEta);
+  fHists->Add(fRecDPosEta);
+  fHists->Add(fRecDNegEta);
+
+
+  for(Int_t i=0;i<18;i++){
+    fHists->Add(fRecPtTpcSector[i]);
+    fHists->Add(fRecEtaTpcSector[i]);
+    fHists->Add(fSignedDcaTpcSector[i]);
+    fHists->Add(fRecQPtTpcSector[i]);
+    fHists->Add(fRecEtaPtTpcSector[i]);
+  }
+
+  for(Int_t i=0;i<7;i++){
+    fHists->Add(fRecPtPosLadder[i]);
+    fHists->Add(fRecPtNegLadder[i]);
+    fHists->Add(fRecPhiPosLadder[i]);
+    fHists->Add(fRecPhiNegLadder[i]);
+    fHists->Add(fRecEtaPosLadder[i]);
+    fHists->Add(fRecEtaNegLadder[i]);
+  }
+
+  fHists->Add(fRecPtPosVz);
+  fHists->Add(fRecPtNegVz);
+  fHists->Add(fRecEtaPosVz);
+  fHists->Add(fRecEtaNegVz);
+  fHists->Add(fRecPhiPosVz);
+  fHists->Add(fRecPhiNegVz);
+  fHists->Add(fSignedDcaPosVz);
+  fHists->Add(fSignedDcaNegVz);
+  fHists->Add(fRecQPtPosEtaVz);
+  fHists->Add(fRecQPtNegEtaVz);
+  fHists->Add(fRecEtaPtPosVz);
+  fHists->Add(fRecEtaPtNegVz);
+
+
+  for(Int_t i=0;i<7;i++){
+    fHists->Add(fSignDcaPos[i]);
+    fHists->Add(fSignDcaNeg[i]);
+    fHists->Add(fSignDcaNegInv[i]);
+    fHists->Add(fPtSigmaPos[i]);
+    fHists->Add(fPtSigmaNeg[i]);
+    fHists->Add(fqPtRec[i]);
+  
+    fHists->Add(fDcaSigmaPos[i]);
+    fHists->Add(fDcaSigmaNeg[i]);
+
+  } 
+  
+    
+    
+  for (Int_t i=0; i<fHists->GetEntries(); ++i) {
+    TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));
+    if (h1){
+      // Printf("%s ",h1->GetName());
+      h1->Sumw2();
+    }
+  }
+    // BKC
+
+
+}
+
+//__________________________________________________________
+
+void AliAnalysisTaskQASym::UserExec(Option_t *) 
+{
+  AliVEvent *event = InputEvent();
+  if (!event) {
+     Printf("ERROR: Could not retrieve event");
+     return;
+  }
+
+
+  if(Entry()==0){
+    AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
+    if(esd){
+      Printf("We are reading from ESD");
+    }
+   
+  }
+
+  Printf("There are %d tracks in this event", event->GetNumberOfTracks());
+
+  
+  Int_t   LeadingTrack  =   0;
+  Float_t LeadingEnergy = -20.;
+  Float_t LeadingPhi    =   0;//TMath::Pi();
+
+
+  if(event->GetNumberOfTracks()!=0) fNumber->Fill(event->GetNumberOfTracks());
+
+  for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
+    
+
+    AliVParticle *track = event->GetTrack(iTrack);
+    AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);
+    if (!track) {
+      Printf("ERROR: Could not receive track %d", iTrack);
+      continue;
+    }
+    
+    if (!fCuts->AcceptTrack(esdtrack)) continue;
+    const AliExternalTrackParam * tpcP = esdtrack->GetTPCInnerParam();
+    if (!tpcP) continue;
+   
+  
+    
+    if(tpcP->E()>LeadingEnergy){
+      LeadingTrack=iTrack;
+      LeadingEnergy=tpcP->E();
+      LeadingPhi=tpcP->Phi();
+    }
+   
+
+
+    
+    //propagate to dca
+    esdtrack->PropagateToDCA(event->GetPrimaryVertex(),
+                            event->GetMagneticField(), 10000.);
+    
+    
+    // if(tpcP->Pt()<2.)continue;
+
+
+    fqRec->Fill(tpcP->Charge());
+  
+
+    Double_t sigmapt = tpcP->GetSigma1Pt2();
+    sigmapt= sqrt(sigmapt);
+    sigmapt= sigmapt *(tpcP->Pt()*tpcP->Pt()); 
+
+    if(sigmapt == 0.)continue;
+    fsigmaPt->Fill(TMath::Log10(sigmapt));
+   
+
+   
+  
+
+    
+
+    // hits in ITS layer
+    Int_t cas=-1;
+    if(esdtrack->HasPointOnITSLayer(0)) 
+      cas=0;
+    else if(!esdtrack->HasPointOnITSLayer(0)
+           &&  esdtrack->HasPointOnITSLayer(1)) 
+      cas=1;
+    else if(!esdtrack->HasPointOnITSLayer(0)
+           && !esdtrack->HasPointOnITSLayer(1) 
+           &&  esdtrack->HasPointOnITSLayer(2)) 
+      cas=2;
+    else if(!esdtrack->HasPointOnITSLayer(0)
+           && !esdtrack->HasPointOnITSLayer(1) 
+           && !esdtrack->HasPointOnITSLayer(2)
+           &&  esdtrack->HasPointOnITSLayer(3)) 
+      cas=3;
+    else if(!esdtrack->HasPointOnITSLayer(0)
+           && !esdtrack->HasPointOnITSLayer(1) 
+           && !esdtrack->HasPointOnITSLayer(2)
+           && !esdtrack->HasPointOnITSLayer(3)
+           &&  esdtrack->HasPointOnITSLayer(4)) 
+      cas=4;
+    else if(   !esdtrack->HasPointOnITSLayer(0)
+           && !esdtrack->HasPointOnITSLayer(1)
+           && !esdtrack->HasPointOnITSLayer(2)
+           && !esdtrack->HasPointOnITSLayer(3)
+           && !esdtrack->HasPointOnITSLayer(4) 
+           &&  esdtrack->HasPointOnITSLayer(5)) 
+      cas=5;
+    else 
+      cas=6;
+  
+   
+   
+    //------------------- 
+
+    sdca_tr = (tpcP->Py()*tpcP->Xv()
+              - tpcP->Px()*tpcP->Yv())/tpcP->Pt();
+  
+
+
+    fqPtRec[cas]->Fill(tpcP->Charge()/tpcP->Pt());
+    
+    
+    test->Fill(tpcP->Pt(),tpcP->Phi(), tpcP->Theta(), 
+              tpcP->Xv(),tpcP->Yv(), tpcP->Zv(), tpcP->Charge());
+    fHistRECpt->Fill(tpcP->Pt());
+    fEta->Fill(tpcP->Eta());
+    fThetaRec->Fill(tpcP->Theta());
+    fPhiRec->Fill(tpcP->Phi());
+    fVx->Fill(tpcP->Xv());
+    fVy->Fill(tpcP->Yv());
+    fVz->Fill(tpcP->Zv());
+  
+
+    fEtaPt->Fill(tpcP->Eta()/tpcP->Pt());
+    fQPt->Fill(tpcP->Charge()/tpcP->Pt());
+    fDca->Fill(sdca_tr);
+
+
+
+
+    esdtrack->GetImpactParameters(xy,z);
+    fDiffDcaD->Fill(sdca_tr+xy);
+   
+    
+    //for positive particles
+    if(tpcP->Charge()>0){
+      fRecPtPos->Fill(tpcP->Pt());
+      fRecPtPosLadder[cas]->Fill(tpcP->Pt());
+      fRecPtPosVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
+      fRecPhiPos->Fill(TMath::RadToDeg()*tpcP->Phi());
+    
+     
+      fRecPhiPosLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
+      fRecPhiPosVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
+      fSignedDcaPosVz->Fill(sdca_tr,tpcP->Zv());
+
+      fRecEtaPos->Fill(tpcP->Eta());
+      fRecEtaPosLadder[cas]->Fill(tpcP->Eta());
+      fRecEtaPtPos->Fill(tpcP->Eta()/tpcP->Pt());
+      fRecEtaPosVz->Fill(tpcP->Eta(),tpcP->Zv());
+      fRecEtaPtPosVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
+     
+      fRecDcaPos->Fill(sdca_tr);
+      fRecDPos->Fill(xy);
+      fSignDcaPos[cas]->Fill(sdca_tr);
+    
+     
+      fDcaSigmaPos[cas]->Fill(sdca_tr, TMath::Log10(sigmapt));
+    
+      fPtSigmaPos[cas]->Fill(TMath::Log10(sigmapt));
+    }
+    //and negative particles
+    else {
+      fRecPtNeg->Fill(tpcP->Pt());
+      fRecPtNegLadder[cas]->Fill(tpcP->Pt());
+      fRecPtNegVz->Fill(TMath::Log10(tpcP->Pt()),tpcP->Zv());
+           
+      fRecPhiNeg->Fill(TMath::RadToDeg()*tpcP->Phi());
+      fRecPhiNegLadder[cas]->Fill(TMath::RadToDeg()*tpcP->Phi());
+      fRecPhiNegVz->Fill(TMath::RadToDeg()*tpcP->Phi(),tpcP->Zv());
+      fSignedDcaNegVz->Fill(sdca_tr,tpcP->Zv());
+      fRecEtaPtNegVz->Fill(tpcP->Eta()/tpcP->Pt(),tpcP->Zv());
+
+      fRecEtaNeg->Fill(tpcP->Eta());
+      fRecEtaNegLadder[cas]->Fill(tpcP->Eta());
+      fRecEtaPtNeg->Fill(tpcP->Eta()/tpcP->Pt());
+      fRecEtaNegVz->Fill(tpcP->Eta(),tpcP->Zv());
+     
+      fRecDcaNeg->Fill(sdca_tr);
+      fRecDcaNegInv->Fill(-sdca_tr);
+      fRecDNeg->Fill(xy);
+      fSignDcaNeg[cas]->Fill(sdca_tr);
+      fSignDcaNegInv[cas]->Fill(-sdca_tr);
+     
+     
+      fDcaSigmaNeg[cas]->Fill(sdca_tr,TMath::Log10(sigmapt));
+   
+      fPtSigmaNeg[cas]->Fill(TMath::Log10(sigmapt));
+    }
+    
+
+
+    //all particles with positive eta
+    if(tpcP->Eta()>0){
+      fRecQPtPosEta->Fill(tpcP->Charge()/tpcP->Pt());
+      fRecPtPosEta->Fill(tpcP->Pt());
+      fRecPhiPosEta->Fill(TMath::RadToDeg()*tpcP->Phi());
+      fRecQPtPosEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
+      fRecDcaPosEta->Fill(sdca_tr);
+      fRecDPosEta->Fill(xy);
+    }
+    //all particles with negative eta (and eta==0)
+    else{
+      fRecQPtNegEta->Fill(tpcP->Charge()/tpcP->Pt());
+      fRecPtNegEta->Fill(tpcP->Pt());
+      fRecPhiNegEta->Fill(TMath::RadToDeg()*tpcP->Phi());
+      fRecQPtNegEtaVz->Fill(tpcP->Charge()/tpcP->Pt(),tpcP->Zv());
+      fRecDcaNegEta->Fill(sdca_tr);
+      fRecDNegEta->Fill(xy);
+
+    }
+     
+
+    //spectren detected by TPC sectors
+    //pt cut on 1 GeV/c ?!
+    // if(tpcP->Pt()<1.) continue;
+    fRecPtTpcSector[Int_t(tpcP->Phi()*
+                         TMath::RadToDeg()/20)]->Fill(tpcP->Pt());
+    fRecEtaTpcSector[Int_t(tpcP->Phi()*
+                         TMath::RadToDeg()/20)]->Fill(tpcP->Eta());
+    fSignedDcaTpcSector[Int_t(tpcP->Phi()*
+                         TMath::RadToDeg()/20)]->Fill(sdca_tr); 
+    fRecQPtTpcSector[Int_t(tpcP->Phi()*
+                         TMath::RadToDeg()/20)]->Fill(tpcP->Charge()/tpcP->Pt());
+    fRecEtaPtTpcSector[Int_t(tpcP->Phi()*
+                         TMath::RadToDeg()/20)]->Fill(tpcP->Eta()/tpcP->Pt());
+     
+
+
+
+
+
+  // another track loop
+    for (Int_t iTrack2 = 0; iTrack2 < event->GetNumberOfTracks(); iTrack2++) {
+      
+      if(LeadingTrack==iTrack2) continue;
+
+      AliVParticle *track2 = event->GetTrack(iTrack2);
+      AliESDtrack* esdtrack2 =  dynamic_cast<AliESDtrack*>(track2);
+      if (!track2) {
+       Printf("ERROR: Could not receive track %d", iTrack);
+       continue;
+      }
+      if (!fCuts->AcceptTrack(esdtrack2)) continue;
+      //propagate to dca
+      esdtrack2->PropagateToDCA(event->GetPrimaryVertex(),
+                              event->GetMagneticField(), 10000.);
+      fDeltaPhiLeading->Fill((LeadingPhi-esdtrack2->Phi())*TMath::RadToDeg(),
+                            LeadingPhi*TMath::RadToDeg() );
+
+     
+
+    }//second track loop
+  }//first track loop
+
+  
+
+  // Post output data.
+  // PostData(1, fHistPt);
+  PostData(1, fHists);
+}      
+
+
+
+
+
+//________________________________________________________________________
+void AliAnalysisTaskQASym::Terminate(Option_t *) 
+{
+
+
+}  
+
+
+
+
+
diff --git a/PWG1/AliAnalysisTaskQASym.h b/PWG1/AliAnalysisTaskQASym.h
new file mode 100644 (file)
index 0000000..b087aaa
--- /dev/null
@@ -0,0 +1,150 @@
+#ifndef AliAnalysisTaskQASym_cxx
+#define AliAnalysisTaskQASym_cxx
+
+class TH1F;
+class TH2F;
+class TList;
+class TNtuple;
+
+class AliESDEvent;
+class AliESDtrack;
+class AliESDtrackCuts;
+
+
+#include "AliAnalysisTaskSE.h"
+#include "TFile.h"
+#include "TNtuple.h"
+
+class AliAnalysisTaskQASym : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskQASym(const char *name = "AliAnalysisTaskQASym");
+  virtual ~AliAnalysisTaskQASym() {}
+  
+  virtual void   UserCreateOutputObjects();
+  virtual void   UserExec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+
+  
+  
+  virtual void   SetCuts(AliESDtrackCuts* cuts)
+     {fCuts = cuts;}
+
+  virtual void   SetFieldOn(Bool_t b = kTRUE){fFieldOn = b;} 
+  
+ private:
+
+  Bool_t      fFieldOn;
+
+  TList       *fHists;          // List of histos
+
+  //old
+  TH1F        *fHistRECpt;      // pt 
+  TH1F        *fEta;            // eta
+  TH1F        *fEtaPt;          // eta over pt 
+  TH1F        *fQPt;            // charge over pt 
+  TH1F        *fDca;            // distance of closest approach
+  TH1F        *fqPtRec[7];      // charge over pt divided for ITS layer cases
+  TH1F        *fqRec;           // reconstrcuted charge
+  TH1F        *fsigmaPt;        // sigma_pT
+  TH2F        *fDcaSigmaPos[7]; // distance of closest approach against sigma_pT
+  TH2F        *fDcaSigmaNeg[7]; // distance of closest approach against sigma_pT
+
+
+  //positive und negative particles
+  TH1F        *fRecPtPos;      // pt of pos partice
+  TH1F        *fRecPtNeg;      // pt of neg particle
+  TH1F        *fRecPhiPos;     // phi of pos. particle
+  TH1F        *fRecPhiNeg;     // phi of neg. particle
+  TH1F        *fRecEtaPos;     // eta of neg. particle
+  TH1F        *fRecEtaNeg;     // eta of neg. particle
+  TH1F        *fRecEtaPtPos;   // eta over pt of neg. particle
+  TH1F        *fRecEtaPtNeg;   // eta over pt of neg. particle
+  TH1F        *fRecDcaPos;     // distance of closest approach of neg. particle
+  TH1F        *fRecDcaNeg;     // distance of closest of neg. particle
+  TH1F        *fRecDcaNegInv;  // invers dca of neg. particle
+  TH1F        *fRecDPos;       // impact parameter of neg. particle
+  TH1F        *fRecDNeg;       // impact parameter of neg. particle
+
+
+
+
+  // two sides of TPC -> Eta/Theta
+  TH1F        *fRecQPtPosEta; 
+  TH1F        *fRecQPtNegEta;
+  TH1F        *fRecPtPosEta;  
+  TH1F        *fRecPtNegEta; 
+  TH1F        *fRecPhiPosEta;
+  TH1F        *fRecPhiNegEta;
+  TH1F        *fRecDcaPosEta;
+  TH1F        *fRecDcaNegEta;
+  TH1F        *fRecDPosEta;
+  TH1F        *fRecDNegEta;
+
+
+  // sectors of TPC (with pt>xGeV?), TODO: extent to TPC standalone tracks
+  TH1F        *fRecPtTpcSector[18];
+  TH1F        *fRecEtaTpcSector[18];
+  TH1F        *fRecQPtTpcSector[18];
+  TH1F        *fRecEtaPtTpcSector[18];
+  TH1F        *fSignedDcaTpcSector[18];
+
+
+  // 7 different case of hit in ITS ladders
+  TH1F        *fSignDcaPos[7];
+  TH1F        *fSignDcaNeg[7];
+  TH1F        *fSignDcaNegInv[7];
+  TH1F        *fPtSigmaPos[7];
+  TH1F        *fPtSigmaNeg[7];
+  TH1F        *fRecPtPosLadder[7];
+  TH1F        *fRecPtNegLadder[7];
+  TH1F        *fRecPhiPosLadder[7];
+  TH1F        *fRecPhiNegLadder[7];
+  TH1F        *fRecEtaPosLadder[7];
+  TH1F        *fRecEtaNegLadder[7];
+
+
+  // 2D: all measures as function of z of vertex
+  TH2F        *fRecPtPosVz;
+  TH2F        *fRecPtNegVz;
+  TH2F        *fRecEtaPosVz;
+  TH2F        *fRecEtaNegVz;
+  TH2F        *fRecPhiPosVz;
+  TH2F        *fRecPhiNegVz;
+  TH2F        *fSignedDcaPosVz;
+  TH2F        *fSignedDcaNegVz;
+  TH2F        *fRecQPtPosEtaVz;
+  TH2F        *fRecQPtNegEtaVz;
+  TH2F        *fRecEtaPtPosVz;
+  TH2F        *fRecEtaPtNegVz;
+
+
+  //high
+  TH1F * fDeltaPhiAll;
+  TH2F * fDeltaPhiLeading;
+  TH1F * fDiffDcaD;
+
+  //sim
+  TH1F * fPhiRec;
+  TH1F * fThetaRec;
+  TH1F * fNumber;
+  TH1F * fVx;
+  TH1F * fVy;
+  TH1F * fVz;
+  TNtuple * test;
+
+
+
+  Double_t  sdca_tr;
+  Float_t xy, z;
+
+  AliESDtrackCuts* fCuts;                      // List of cuts
+  AliAnalysisTaskQASym(const AliAnalysisTaskQASym&); // not implemented
+  AliAnalysisTaskQASym& operator=(const AliAnalysisTaskQASym&); // not implemented
+  
+  ClassDef(AliAnalysisTaskQASym, 1); // Basic QA exploiting symmetries
+};
+
+#endif