PWGPP/AliAnalysisTaskITSTPCalignment.cxx

   1 ////////////////////////////////////////////////////////////////////////////////
   2 //  AliAnalysisTaskITSTPCalignment
   3 //  Runs the relative ITS TPC alignment procedure and TPC vdrift calib
   4 //  Origin: Mikolaj Krzewicki, mikolaj.krzewicki@cern.ch
   5 ////////////////////////////////////////////////////////////////////////////////
   6
   7 #include <TChain.h>
   8 #include <TTree.h>
   9 #include <TH2.h>
  10 #include <TProfile.h>
  11 #include <TCanvas.h>
  12 #include <TArrayI.h>
  13 #include <TObjArray.h>
  14 #include <TGraphErrors.h>
  15 #include <AliMagF.h>
  16 #include <AliAnalysisTaskSE.h>
  17 #include <AliMCEventHandler.h>
  18 #include <AliAnalysisManager.h>
  19 #include <AliESDEvent.h>
  20 #include <AliESDfriend.h>
  21 #include <AliMCEvent.h>
  22 #include <AliStack.h>
  23 #include <AliExternalTrackParam.h>
  24 #include <AliESDtrack.h>
  25 #include <AliESDfriendTrack.h>
  26 #include <AliESDInputHandler.h>
  27 #include "AliRelAlignerKalman.h"
  28 #include "AliRelAlignerKalmanArray.h"
  29 #include "AliAnalysisTaskITSTPCalignment.h"
  30 #include "AliLog.h"
  31
  32 ClassImp(AliAnalysisTaskITSTPCalignment)
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment():
  36     AliAnalysisTaskSE(),
  37     fArrayITSglobal(0),
  38     fArrayITSsa(0),
  39     fDebugTree(0),
  40     fAligner(0),
  41     fList(0),
  42     fFillDebugTree(kFALSE),
  43     fDoQA(kTRUE),
  44     fT0(0),
  45     fTend(0),
  46     fSlotWidth(0),
  47     fMinPt(0.4),
  48     fMinPointsVol1(3),
  49     fMinPointsVol2(80),
  50     fRejectOutliers(kTRUE),
  51     fOutRejSigma(1.),
  52     fRejectOutliersSigma2Median(kFALSE),
  53     fOutRejSigma2Median(5.),
  54     fOutRejSigmaOnMerge(10.),
  55     fUseITSoutGlobalTrack(kFALSE),
  56     fUseITSoutITSSAtrack(kTRUE)
  57 {
  58   //dummy ctor
  59   // Define input and output slots here
  60   // Input slot #0 works with a TChain
  61   //DefineInput(0, TChain::Class());
  62   //DefineOutput(0, TTree::Class());
  63   //DefineOutput(1, TList::Class());
  64   //DefineOutput(2, AliRelAlignerKalmanArray::Class());
  65 }
  66
  67 //________________________________________________________________________
  68 AliAnalysisTaskITSTPCalignment::AliAnalysisTaskITSTPCalignment(const char *name):
  69     AliAnalysisTaskSE(name),
  70     fArrayITSglobal(0),
  71     fArrayITSsa(0),
  72     fDebugTree(0),
  73     fAligner(0),
  74     fList(0),
  75     fFillDebugTree(kFALSE),
  76     fDoQA(kTRUE),
  77     fT0(0),
  78     fTend(0),
  79     fSlotWidth(0),
  80     fMinPt(0.4),
  81     fMinPointsVol1(3),
  82     fMinPointsVol2(80),
  83     fRejectOutliers(kTRUE),
  84     fOutRejSigma(1.),
  85     fRejectOutliersSigma2Median(kFALSE),
  86     fOutRejSigma2Median(5.),
  87     fOutRejSigmaOnMerge(10.),
  88     fUseITSoutGlobalTrack(kFALSE),
  89     fUseITSoutITSSAtrack(kTRUE)
  90 {
  91   // Constructor
  92
  93   // Define input and output slots here
  94   // Input slot #0 is a TChain
  95   // Output slot #0 is a TTree
  96   DefineOutput(1, TList::Class());
  97   DefineOutput(2, AliRelAlignerKalmanArray::Class());
  98   DefineOutput(3, AliRelAlignerKalmanArray::Class());
  99 }
 100
 101 //________________________________________________________________________
 102 Bool_t AliAnalysisTaskITSTPCalignment::UserNotify()
 103 {
 104   //ON INPUT FILE/TREE CHANGE
 105
 106   //fArray->Print();
 107   if (fFillDebugTree)
 108   {
 109     if (fAligner->GetNUpdates()>0) fDebugTree->Fill();
 110   }
 111
 112   return kTRUE;
 113 }
 114
 115 //________________________________________________________________________
 116 void AliAnalysisTaskITSTPCalignment::UserCreateOutputObjects()
 117 {
 118   // Create output objects
 119
 120   fArrayITSglobal = new AliRelAlignerKalmanArray();
 121   fArrayITSglobal->SetName("array ITS global");
 122   fArrayITSglobal->SetupArray(fT0,fTend,fSlotWidth);
 123   fArrayITSglobal->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge);
 124   //set up the template
 125   AliRelAlignerKalman* templ = fArrayITSglobal->GetAlignerTemplate();
 126   templ->SetRejectOutliers(fRejectOutliers);
 127   templ->SetOutRejSigma(fOutRejSigma);
 128   templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
 129   templ->SetOutRejSigma2Median(fOutRejSigma2Median);
 130   fArrayITSsa = new AliRelAlignerKalmanArray();
 131   fArrayITSsa->SetName("array ITS SA");
 132   fArrayITSsa->SetupArray(fT0,fTend,fSlotWidth);
 133   fArrayITSsa->SetOutRejSigmaOnMerge(fOutRejSigmaOnMerge);
 134   //set up the template
 135   templ = fArrayITSsa->GetAlignerTemplate();
 136   templ->SetRejectOutliers(fRejectOutliers);
 137   templ->SetOutRejSigma(fOutRejSigma);
 138   templ->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
 139   templ->SetOutRejSigma2Median(fOutRejSigma2Median);
 140
 141   fList = new TList();
 142   fList->SetName("QA");
 143   fList->SetOwner(kTRUE);
 144
 145   TH2F* pZYAResidualsHistBpos = new TH2F("fZYAResidualsHistBpos","z-r\\phi residuals side A (z>0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
 146   pZYAResidualsHistBpos->SetXTitle("\\deltaz [cm]");
 147   pZYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 148   TH2F* pZYCResidualsHistBpos = new TH2F("fZYCResidualsHistBpos","z-r\\phi residuals side C (z<0), Bpos", 100, -4., 4., 100, -1.5, 1.5 );
 149   pZYCResidualsHistBpos->SetXTitle("\\deltaz [cm]");
 150   pZYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 151   TH2F* pLPAResidualsHistBpos = new TH2F("fLPAResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
 152   pLPAResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
 153   pLPAResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
 154   TH2F* pLPCResidualsHistBpos = new TH2F("fLPCResidualsHistBpos", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bpos", 100, -.05, 0.05, 100, -0.05, 0.05 );
 155   pLPCResidualsHistBpos->SetXTitle("\\deltasin(\\phi)");
 156   pLPCResidualsHistBpos->SetYTitle("\\deltatan(\\lambda)");
 157   TH2F* pPhiYAResidualsHistBpos = new TH2F("fPhiYAResidualsHistBpos","\\phi-\\deltar\\phi side A (z>0), Bpos",36,0.0,6.3,100,-1.5,1.5);
 158   pPhiYAResidualsHistBpos->SetXTitle("\\phi [rad]");
 159   pPhiYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 160   TH2F* pPhiYCResidualsHistBpos = new TH2F("fPhiYCResidualsHistBpos","\\phi-\\deltar\\phi side C (z<0), Bpos",36,0.0,6.3,100,-1.5,1.5);
 161   pPhiYCResidualsHistBpos->SetXTitle("\\phi [rad]");
 162   pPhiYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 163   TH2F* pPhiZAResidualsHistBpos = new TH2F("fPhiZAResidualsHistBpos","\\phi-\\deltaz side A (z>0), Bpos",36,0.0,6.3,100,-3.0,3.0);
 164   pPhiZAResidualsHistBpos->SetXTitle("\\phi [rad]");
 165   pPhiZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 166   TH2F* pPhiZCResidualsHistBpos = new TH2F("fPhiZCResidualsHistBpos","\\phi-\\deltaz side C (z<0), Bpos",36,0.0,6.3,100,-3.0,3.0);
 167   pPhiZCResidualsHistBpos->SetXTitle("\\phi [rad]");
 168   pPhiZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 169   TH2F* pPtYAResidualsHistBpos = new TH2F("fPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",20,.3,8.,80,-1.5,1.5);
 170   pPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
 171   pPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 172   TH2F* pPtYCResidualsHistBpos = new TH2F("fPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",20,.3,8.,80,-1.5,1.5);
 173   pPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
 174   pPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 175   TH2F* pPtZAResidualsHistBpos = new TH2F("fPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",20,.3,8.,80,-3.0,3.0);
 176   pPtZAResidualsHistBpos->SetXTitle("Pt");
 177   pPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 178   TH2F* pPtZCResidualsHistBpos = new TH2F("fPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",20,.3,8.,80,-3.0,3.0);
 179   pPtZCResidualsHistBpos->SetXTitle("Pt");
 180   pPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 181   TH2F* pLowPtYAResidualsHistBpos = new TH2F("fLowPtYAResidualsHistBpos","Pt-\\deltar\\phi side A (z>0), Bpos",100,0.0,1.5,80,-1.5,1.5);
 182   pLowPtYAResidualsHistBpos->SetXTitle("Pt [rad]");
 183   pLowPtYAResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 184   TH2F* pLowPtYCResidualsHistBpos = new TH2F("fLowPtYCResidualsHistBpos","Pt-\\deltar\\phi side C (z<0), Bpos",100,0.0,1.5,80,-1.5,1.5);
 185   pLowPtYCResidualsHistBpos->SetXTitle("Pt [rad]");
 186   pLowPtYCResidualsHistBpos->SetYTitle("\\deltar\\phi [cm]");
 187   TH2F* pLowPtZAResidualsHistBpos = new TH2F("fLowPtZAResidualsHistBpos","Pt-\\deltaz side A (z>0), Bpos",100,0.0,1.5,80,-3.0,3.0);
 188   pLowPtZAResidualsHistBpos->SetXTitle("Pt");
 189   pLowPtZAResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 190   TH2F* pLowPtZCResidualsHistBpos = new TH2F("fLowPtZCResidualsHistBpos","Pt-\\deltaz side C (z<0), Bpos",100,0.0,1.5,80,-3.0,3.0);
 191   pLowPtZCResidualsHistBpos->SetXTitle("Pt");
 192   pLowPtZCResidualsHistBpos->SetYTitle("\\deltaz [cm]");
 193   TList* listBpos = new TList();
 194   listBpos->SetName("B+");
 195   listBpos->SetOwner(kTRUE);
 196   fList->Add(listBpos); //0
 197   listBpos->Add(pZYAResidualsHistBpos);
 198   listBpos->Add(pZYCResidualsHistBpos);
 199   listBpos->Add(pLPAResidualsHistBpos);
 200   listBpos->Add(pLPCResidualsHistBpos);
 201   listBpos->Add(pPhiYAResidualsHistBpos);
 202   listBpos->Add(pPhiYCResidualsHistBpos);
 203   listBpos->Add(pPhiZAResidualsHistBpos);
 204   listBpos->Add(pPhiZCResidualsHistBpos);
 205   listBpos->Add(pPtYAResidualsHistBpos);
 206   listBpos->Add(pPtYCResidualsHistBpos);
 207   listBpos->Add(pPtZAResidualsHistBpos);
 208   listBpos->Add(pPtZCResidualsHistBpos);
 209   listBpos->Add(pLowPtYAResidualsHistBpos);
 210   listBpos->Add(pLowPtYCResidualsHistBpos);
 211   listBpos->Add(pLowPtZAResidualsHistBpos);
 212   listBpos->Add(pLowPtZCResidualsHistBpos);
 213
 214   TH2F* pZYAResidualsHistBneg = new TH2F("fZYAResidualsHistBneg","z-r\\phi residuals side A (z>0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
 215   pZYAResidualsHistBneg->SetXTitle("\\deltaz [cm]");
 216   pZYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 217   TH2F* pZYCResidualsHistBneg = new TH2F("fZYCResidualsHistBneg","z-r\\phi residuals side C (z<0), Bneg", 100, -4., 4., 100, -1.5, 1.5 );
 218   pZYCResidualsHistBneg->SetXTitle("\\deltaz [cm]");
 219   pZYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 220   TH2F* pLPAResidualsHistBneg = new TH2F("fLPAResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
 221   pLPAResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
 222   pLPAResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
 223   TH2F* pLPCResidualsHistBneg = new TH2F("fLPCResidualsHistBneg", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bneg", 100, -.05, 0.05, 100, -0.05, 0.05 );
 224   pLPCResidualsHistBneg->SetXTitle("\\deltasin(\\phi)");
 225   pLPCResidualsHistBneg->SetYTitle("\\deltatan(\\lambda)");
 226   TH2F* pPhiYAResidualsHistBneg = new TH2F("fPhiYAResidualsHistBneg","\\phi-\\deltar\\phi side A (z>0), Bneg",36,0.0,6.3,100,-1.5,1.5);
 227   pPhiYAResidualsHistBneg->SetXTitle("\\phi [rad]");
 228   pPhiYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 229   TH2F* pPhiYCResidualsHistBneg = new TH2F("fPhiYCResidualsHistBneg","\\phi-\\deltar\\phi side C (z<0), Bneg",36,0.0,6.3,100,-1.5,1.5);
 230   pPhiYCResidualsHistBneg->SetXTitle("\\phi [rad]");
 231   pPhiYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 232   TH2F* pPhiZAResidualsHistBneg = new TH2F("fPhiZAResidualsHistBneg","\\phi-\\deltaz side A (z>0), Bneg",36,0.0,6.3,100,-3.0,3.0);
 233   pPhiZAResidualsHistBneg->SetXTitle("\\phi [rad]");
 234   pPhiZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 235   TH2F* pPhiZCResidualsHistBneg = new TH2F("fPhiZCResidualsHistBneg","\\Pt-\\deltaz side C (z<0), Bneg",36,0.0,6.3,100,-3.0,3.0);
 236   pPhiZCResidualsHistBneg->SetXTitle("\\phi [rad]");
 237   pPhiZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 238   TH2F* pPtYAResidualsHistBneg = new TH2F("fPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",20,.3,8.,80,-1.5,1.5);
 239   pPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
 240   pPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 241   TH2F* pPtYCResidualsHistBneg = new TH2F("fPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",20,.3,8.,80,-1.5,1.5);
 242   pPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
 243   pPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 244   TH2F* pPtZAResidualsHistBneg = new TH2F("fPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",20,.3,8.,80,-3.0,3.0);
 245   pPtZAResidualsHistBneg->SetXTitle("Pt");
 246   pPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 247   TH2F* pPtZCResidualsHistBneg = new TH2F("fPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",20,.3,8.,80,-3.0,3.0);
 248   pPtZCResidualsHistBneg->SetXTitle("Pt");
 249   pPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 250   TH2F* pLowPtYAResidualsHistBneg = new TH2F("fLowPtYAResidualsHistBneg","Pt-\\deltar\\phi side A (z>0), Bneg",100,0.0,1.5,80,-1.5,1.5);
 251   pLowPtYAResidualsHistBneg->SetXTitle("Pt [rad]");
 252   pLowPtYAResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 253   TH2F* pLowPtYCResidualsHistBneg = new TH2F("fLowPtYCResidualsHistBneg","Pt-\\deltar\\phi side C (z<0), Bneg",100,0.0,1.5,80,-1.5,1.5);
 254   pLowPtYCResidualsHistBneg->SetXTitle("Pt [rad]");
 255   pLowPtYCResidualsHistBneg->SetYTitle("\\deltar\\phi [cm]");
 256   TH2F* pLowPtZAResidualsHistBneg = new TH2F("fLowPtZAResidualsHistBneg","Pt-\\deltaz side A (z>0), Bneg",100,0.0,1.5,80,-3.0,3.0);
 257   pLowPtZAResidualsHistBneg->SetXTitle("Pt");
 258   pLowPtZAResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 259   TH2F* pLowPtZCResidualsHistBneg = new TH2F("fLowPtZCResidualsHistBneg","Pt-\\deltaz side C (z<0), Bneg",100,0.0,1.5,80,-3.0,3.0);
 260   pLowPtZCResidualsHistBneg->SetXTitle("Pt");
 261   pLowPtZCResidualsHistBneg->SetYTitle("\\deltaz [cm]");
 262   TList* listBneg = new TList();
 263   listBneg->SetName("B-");
 264   listBneg->SetOwner(kTRUE);
 265   fList->Add(listBneg); //1
 266   listBneg->Add(pZYAResidualsHistBneg);
 267   listBneg->Add(pZYCResidualsHistBneg);
 268   listBneg->Add(pLPAResidualsHistBneg);
 269   listBneg->Add(pLPCResidualsHistBneg);
 270   listBneg->Add(pPhiYAResidualsHistBneg);
 271   listBneg->Add(pPhiYCResidualsHistBneg);
 272   listBneg->Add(pPhiZAResidualsHistBneg);
 273   listBneg->Add(pPhiZCResidualsHistBneg);
 274   listBneg->Add(pPtYAResidualsHistBneg);
 275   listBneg->Add(pPtYCResidualsHistBneg);
 276   listBneg->Add(pPtZAResidualsHistBneg);
 277   listBneg->Add(pPtZCResidualsHistBneg);
 278   listBneg->Add(pLowPtYAResidualsHistBneg);
 279   listBneg->Add(pLowPtYCResidualsHistBneg);
 280   listBneg->Add(pLowPtZAResidualsHistBneg);
 281   listBneg->Add(pLowPtZCResidualsHistBneg);
 282
 283   TH2F* pZYAResidualsHistBnil = new TH2F("fZYAResidualsHistBnil","z-r\\phi residuals side A (z>0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
 284   pZYAResidualsHistBnil->SetXTitle("\\deltaz [cm]");
 285   pZYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 286   TH2F* pZYCResidualsHistBnil = new TH2F("fZYCResidualsHistBnil","z-r\\phi residuals side C (z<0), Bnil", 100, -3., 3., 100, -1.5, 1.5 );
 287   pZYCResidualsHistBnil->SetXTitle("\\deltaz [cm]");
 288   pZYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 289   TH2F* pLPAResidualsHistBnil = new TH2F("fLPAResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side A (z>0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
 290   pLPAResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
 291   pLPAResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
 292   TH2F* pLPCResidualsHistBnil = new TH2F("fLPCResidualsHistBnil", "sin(\\phi) tan(\\lambda) residuals side C (z<0), Bnil", 100, -.05, 0.05, 100, -0.05, 0.05 );
 293   pLPCResidualsHistBnil->SetXTitle("\\deltasin(\\phi)");
 294   pLPCResidualsHistBnil->SetYTitle("\\deltatan(\\lambda)");
 295   TH2F* pPhiYAResidualsHistBnil = new TH2F("fPhiYAResidualsHistBnil","\\phi-\\deltar\\phi side A (z>0), Bnil",36,0.0,6.3,100,-1.5,1.5);
 296   pPhiYAResidualsHistBnil->SetXTitle("\\phi [rad]");
 297   pPhiYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 298   TH2F* pPhiYCResidualsHistBnil = new TH2F("fPhiYCResidualsHistBnil","\\phi-\\deltar\\phi side C (z<0), Bnil",36,0.0,6.3,100,-1.5,1.5);
 299   pPhiYCResidualsHistBnil->SetXTitle("\\phi [rad]");
 300   pPhiYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 301   TH2F* pPhiZAResidualsHistBnil = new TH2F("fPhiZAResidualsHistBnil","\\phi-\\deltaz side A (z>0), Bnil",36,0.0,6.3,100,-3.0,3.0);
 302   pPhiZAResidualsHistBnil->SetXTitle("\\phi [rad]");
 303   pPhiZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 304   TH2F* pPhiZCResidualsHistBnil = new TH2F("fPhiZCResidualsHistBnil","\\Pt-\\deltaz side C (z<0), Bnil",36,0.0,6.3,100,-3.0,3.0);
 305   pPhiZCResidualsHistBnil->SetXTitle("\\phi [rad]");
 306   pPhiZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 307   TH2F* pPtYAResidualsHistBnil = new TH2F("fPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",20,.3,8.,80,-1.5,1.5);
 308   pPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
 309   pPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 310   TH2F* pPtYCResidualsHistBnil = new TH2F("fPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",20,.3,8.,80,-1.5,1.5);
 311   pPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
 312   pPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 313   TH2F* pPtZAResidualsHistBnil = new TH2F("fPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",20,.3,8.,80,-3.0,3.0);
 314   pPtZAResidualsHistBnil->SetXTitle("Pt");
 315   pPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 316   TH2F* pPtZCResidualsHistBnil = new TH2F("fPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",20,.3,8.,80,-3.0,3.0);
 317   pPtZCResidualsHistBnil->SetXTitle("Pt");
 318   pPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 319   TH2F* pLowPtYAResidualsHistBnil = new TH2F("fLowPtYAResidualsHistBnil","Pt-\\deltar\\phi side A (z>0), Bnil",100,0.0,1.5,80,-1.5,1.5);
 320   pLowPtYAResidualsHistBnil->SetXTitle("Pt [rad]");
 321   pLowPtYAResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 322   TH2F* pLowPtYCResidualsHistBnil = new TH2F("fLowPtYCResidualsHistBnil","Pt-\\deltar\\phi side C (z<0), Bnil",100,0.0,1.5,80,-1.5,1.5);
 323   pLowPtYCResidualsHistBnil->SetXTitle("Pt [rad]");
 324   pLowPtYCResidualsHistBnil->SetYTitle("\\deltar\\phi [cm]");
 325   TH2F* pLowPtZAResidualsHistBnil = new TH2F("fLowPtZAResidualsHistBnil","Pt-\\deltaz side A (z>0), Bnil",100,0.0,1.5,80,-3.0,3.0);
 326   pLowPtZAResidualsHistBnil->SetXTitle("Pt");
 327   pLowPtZAResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 328   TH2F* pLowPtZCResidualsHistBnil = new TH2F("fLowPtZCResidualsHistBnil","Pt-\\deltaz side C (z<0), Bnil",100,0.0,1.5,80,-3.0,3.0);
 329   pLowPtZCResidualsHistBnil->SetXTitle("Pt");
 330   pLowPtZCResidualsHistBnil->SetYTitle("\\deltaz [cm]");
 331   TList* listBnil = new TList();
 332   listBnil->SetName("B0");
 333   listBnil->SetOwner(kTRUE);
 334   fList->Add(listBnil); //2
 335   listBnil->Add(pZYAResidualsHistBnil);
 336   listBnil->Add(pZYCResidualsHistBnil);
 337   listBnil->Add(pLPAResidualsHistBnil);
 338   listBnil->Add(pLPCResidualsHistBnil);
 339   listBnil->Add(pPhiYAResidualsHistBnil);
 340   listBnil->Add(pPhiYCResidualsHistBnil);
 341   listBnil->Add(pPhiZAResidualsHistBnil);
 342   listBnil->Add(pPhiZCResidualsHistBnil);
 343   listBnil->Add(pPtYAResidualsHistBnil);
 344   listBnil->Add(pPtYCResidualsHistBnil);
 345   listBnil->Add(pPtZAResidualsHistBnil);
 346   listBnil->Add(pPtZCResidualsHistBnil);
 347   listBnil->Add(pLowPtYAResidualsHistBnil);
 348   listBnil->Add(pLowPtYCResidualsHistBnil);
 349   listBnil->Add(pLowPtZAResidualsHistBnil);
 350   listBnil->Add(pLowPtZCResidualsHistBnil);
 351
 352   TH1F* pNmatchingEff=new TH1F("pNmatchingEff","matching efficiency",50,0.,1.);
 353   fList->Add(pNmatchingEff); //3
 354
 355   TH1I* pChecks = new TH1I("some counters","some counters",10,0,10);
 356   pChecks->GetXaxis()->SetBinLabel(1+kNoESD,"no ESD");
 357   pChecks->GetXaxis()->SetBinLabel(1+kNoESDfriend,"no ESDfriend");
 358   pChecks->GetXaxis()->SetBinLabel(1+kNoFriendTrack,"no friendtrack");
 359   pChecks->GetXaxis()->SetBinLabel(1+kNoITSoutParams,"no ITS out params");
 360   pChecks->GetXaxis()->SetBinLabel(1+kESDfriend,"ESD friend");
 361   pChecks->GetXaxis()->SetBinLabel(1+kFriendsSkipBit,"friends+skipbit");
 362   pChecks->GetXaxis()->SetBinLabel(1+kFriendTrack,"friendtrack");
 363   pChecks->GetXaxis()->SetBinLabel(1+kITSoutParams,"ITS out params");
 364   fList->Add(pChecks); //4
 365
 366   fAligner = new AliRelAlignerKalman();
 367   fAligner->SetRejectOutliers(fRejectOutliers);
 368   fAligner->SetOutRejSigma(fOutRejSigma);
 369   fAligner->SetRejectOutliersSigma2Median(fRejectOutliersSigma2Median);
 370   fAligner->SetOutRejSigma2Median(fOutRejSigma2Median);
 371   fList->Add(fAligner);
 372
 373   fDebugTree = new TTree("debugTree","tree with debug info");
 374   fDebugTree->Branch("aligner","AliRelAlignerKalman",&fAligner);
 375
 376   // Post output data.
 377   PostData(0, fDebugTree);
 378   PostData(1, fList);
 379   PostData(2, fArrayITSglobal);
 380   PostData(3, fArrayITSsa);
 381 }
 382
 383 //________________________________________________________________________
 384 void AliAnalysisTaskITSTPCalignment::UserExec(Option_t *)
 385 {
 386   // Main loop
 387   // Called for each event
 388
 389   TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
 390   //check for ESD and friends
 391   AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
 392   AliESDfriend* eventFriend = dynamic_cast<AliESDfriend*>(ESDfriend());
 393   if (!event)
 394   {
 395     AliError("no ESD");
 396     pChecks->Fill(kNoESD);
 397     return;
 398   }
 399
 400   if (!eventFriend)
 401   {
 402     AliError("no ESD friend");
 403     pChecks->Fill(kNoESDfriend);
 404     return;
 405   }
 406   else
 407   {
 408     pChecks->Fill(kESDfriend);
 409   }
 410
 411   //event->SetESDfriend(eventFriend);
 412
 413   if (eventFriend->TestSkipBit())
 414   {
 415     pChecks->Fill(kFriendsSkipBit);
 416     return;
 417   }
 418
 419   //Update the parmeters
 420   AnalyzeESDevent(event);
 421 }
 422
 423 //________________________________________________________________________
 424 void AliAnalysisTaskITSTPCalignment::AnalyzeESDevent(AliESDEvent* event)
 425 {
 426   //analyze an ESD event with track matching
 427   Int_t ntracks = event->GetNumberOfTracks();
 428   if (ntracks==0) return;
 429
 430   if (fUseITSoutGlobalTrack)
 431   {
 432     AliRelAlignerKalman* alignerGlobal = fArrayITSglobal->GetAligner(event);
 433     if (alignerGlobal)
 434       alignerGlobal->AddESDevent(event);
 435   }
 436
 437   if (fUseITSoutITSSAtrack)
 438   {
 439     //get the aligner for the event
 440     //do it with matching
 441     TObjArray arrayParamITS(ntracks);
 442     TObjArray arrayParamTPC(ntracks);
 443
 444     Int_t n = FindMatchingTracks(arrayParamITS, arrayParamTPC, event);
 445     AliInfo(Form("matched tracklet pairs: %i\n",n));
 446
 447     TH1F* nMatchingEff=static_cast<TH1F*>(fList->At(3));
 448     Float_t ratio = (float)n/(float)ntracks;
 449     nMatchingEff->Fill(ratio);
 450
 451     AliRelAlignerKalman* alignerSA = NULL;
 452     if (n>0) alignerSA = fArrayITSsa->GetAligner(event);
 453
 454     for (Int_t i=0;i<n;i++)
 455     {
 456       AliExternalTrackParam* paramsITS=static_cast<AliExternalTrackParam*>(arrayParamITS[i]);
 457       AliExternalTrackParam* paramsTPC=static_cast<AliExternalTrackParam*>(arrayParamTPC[i]);
 458
 459       if (!(paramsITS&&paramsTPC)) continue;
 460
 461       //QA
 462       if (fDoQA) DoQA(paramsITS,paramsTPC);
 463
 464       //debug
 465       if (fAligner->AddTrackParams(paramsITS, paramsTPC))
 466       {
 467         fAligner->SetRunNumber(event->GetRunNumber());
 468         fAligner->SetMagField(event->GetMagneticField());
 469         fAligner->SetTimeStamp(event->GetTimeStamp());
 470       }
 471
 472       //alignment
 473       if (alignerSA) alignerSA->AddTrackParams(paramsITS, paramsTPC);
 474     }
 475     arrayParamITS.Delete();
 476     arrayParamTPC.Delete();
 477   }
 478 }
 479
 480 //________________________________________________________________________
 481 void AliAnalysisTaskITSTPCalignment::Terminate(Option_t *)
 482 {
 483   // Called once at the end of the query
 484   fArrayITSsa->Print("a");
 485   fAligner->Print();
 486 }
 487
 488 //________________________________________________________________________
 489 Int_t AliAnalysisTaskITSTPCalignment::FindMatchingTracks(TObjArray& arrITS, TObjArray& arrTPC, AliESDEvent* pESD)
 490 {
 491   //find matching tracks and return tobjarrays with the params
 492   //fUniqueID of param object contains tracknumber in event.
 493
 494   TH1I* pChecks = static_cast<TH1I*>(fList->At(4));
 495   Int_t ntracks = pESD->GetNumberOfTracks();
 496   Double_t magfield = pESD->GetMagneticField();
 497
 498   Int_t* matchedArr = new Int_t[ntracks]; //storage for index of ITS track for which a match was found
 499   for (Int_t i=0;i<ntracks;i++)
 500   {
 501     matchedArr[i]=-1;
 502   }
 503
 504   Int_t iMatched=-1;
 505   for (Int_t i=0; i<ntracks; i++)
 506   {
 507     //get track1 and friend
 508     AliESDtrack* track1 = pESD->GetTrack(i);
 509     if (!track1) continue;
 510
 511     if (track1->GetNcls(0) < fMinPointsVol1) continue;
 512
 513     if (!( ( track1->IsOn(AliESDtrack::kITSrefit)) &&
 514            (!track1->IsOn(AliESDtrack::kTPCin)) )) continue;
 515
 516     const AliESDfriendTrack* constfriendtrack1 = track1->GetFriendTrack();
 517     if (!constfriendtrack1) {AliInfo(Form("no friendTrack1\n"));pChecks->Fill(kNoFriendTrack);continue;}
 518     pChecks->Fill(kFriendTrack);
 519     AliESDfriendTrack friendtrack1(*constfriendtrack1);
 520
 521     if (!friendtrack1.GetITSOut()) {pChecks->Fill(kNoITSoutParams);continue;}
 522     pChecks->Fill(kITSoutParams);
 523     AliExternalTrackParam params1(*(friendtrack1.GetITSOut()));
 524
 525     Double_t bestd = 1000.; //best distance
 526     Bool_t newi = kTRUE; //whether we start with a new i
 527     for (Int_t j=0; j<ntracks; j++)
 528     {
 529       if (matchedArr[j]>0 && matchedArr[j]!=i) continue; //already matched, everything tried
 530       //get track2 and friend
 531       AliESDtrack* track2 = pESD->GetTrack(j);
 532       if (!track2) continue;
 533       if (track1==track2) continue;
 534       if (!(track2->IsOn(AliESDtrack::kTPCin))) continue; //must be TPC
 535       if (!(track2->IsOn(AliESDtrack::kITSout))) continue; //must have ITS
 536
 537       //if (track2->GetNcls(0) != track1->GetNcls(0)) continue;
 538       //if (track2->GetITSClusterMap() != track1->GetITSClusterMap()) continue;
 539       if (track2->GetNcls(1) < fMinPointsVol2) continue; //min 80 clusters in TPC
 540       if (track2->GetTgl() > 1.) continue; //acceptance
 541       //cut crossing tracks
 542       if (!track2->GetInnerParam()) continue;
 543       if (!track2->GetOuterParam()) continue;
 544       if (track2->GetOuterParam()->GetZ()*track2->GetInnerParam()->GetZ()<0) continue;
 545       if (track2->GetInnerParam()->GetX()>90) continue;
 546       if (TMath::Abs(track2->GetInnerParam()->GetZ())<5.) continue; //too close to membrane?
 547
 548       AliExternalTrackParam params2(*(track2->GetInnerParam()));
 549
 550       //bring to same reference plane
 551       if (!params2.Rotate(params1.GetAlpha())) continue;
 552       if (!params2.PropagateTo(params1.GetX(), magfield)) continue;
 553
 554       //pt cut, only for data with magfield
 555       if ((params2.Pt()<fMinPt)&&(TMath::Abs(magfield)>1.)) continue;
 556       //else
 557       //if (fMC)
 558       // {
 559       //   AliStack* stack = fMC->Stack();
 560       //   Int_t label = track2->GetLabel();
 561       //   if (label<0) continue;
 562       //   TParticle* particle = stack->Particle(label);
 563       //   if (particle->Pt()<fMinPt) continue;
 564       // }
 565
 566       const Double32_t* p1 = params1.GetParameter();
 567       const Double32_t* p2 = params2.GetParameter();
 568
 569       //hard cuts
 570       Double_t dy = TMath::Abs(p2[0]-p1[0]);
 571       Double_t dz = TMath::Abs(p2[1]-p1[1]);
 572       Double_t dphi = TMath::Abs(p2[2]-p1[2]);
 573       Double_t dlam = TMath::Abs(p2[3]-p1[3]);
 574       if (dy > 2.0) continue;
 575       if (dz > 10.0) continue;
 576       if (dphi > 0.1 ) continue;
 577       if (dlam > 0.1 ) continue;
 578
 579       //best match only
 580       Double_t d = TMath::Sqrt(dy*dy+dz*dz+dphi*dphi+dlam*dlam);
 581       if ( d >= bestd) continue;
 582       bestd = d;
 583       matchedArr[j]=i; //j-th track matches i-th (ITS) track
 584       if (newi) iMatched++; newi=kFALSE; //increment at most once per i
 585       params1.SetUniqueID(i); //store tracknummer
 586       params2.SetUniqueID(j);
 587       if (arrITS[iMatched] && arrTPC[iMatched])
 588       {
 589         *(arrITS[iMatched]) = params1;
 590         *(arrTPC[iMatched]) = params2;
 591       }
 592       else
 593       {
 594         arrITS[iMatched] = new AliExternalTrackParam(params1);
 595         arrTPC[iMatched] = new AliExternalTrackParam(params2);
 596       }//else
 597     }//for j
 598   }//for i
 599   delete [] matchedArr;
 600   return iMatched+1;
 601 }
 602
 603 //________________________________________________________________________
 604 void AliAnalysisTaskITSTPCalignment::DoQA(AliExternalTrackParam* paramsITS,
 605     AliExternalTrackParam* paramsTPC)
 606 {
 607   //fill qa histograms in a given list, per field direction (5,-5,0)
 608   Float_t resy = paramsTPC->GetY() - paramsITS->GetY();
 609   Float_t resz = paramsTPC->GetZ() - paramsITS->GetZ();
 610   Float_t ressnp = paramsTPC->GetSnp() - paramsITS->GetSnp();
 611   Float_t restgl = paramsTPC->GetTgl() - paramsITS->GetTgl();
 612
 613   TList* pList=NULL;
 614   Double_t field = fInputEvent->GetMagneticField();
 615   if (field >= 1.)  pList = static_cast<TList*>(fList->At(0));
 616   if (field <= -1.) pList = static_cast<TList*>(fList->At(1));
 617   if (field < 1. && field > -1.) pList = static_cast<TList*>(fList->At(2));
 618   if (!pList) return;
 619
 620   TH2F* pZYAResidualsHist = static_cast<TH2F*>(pList->At(0));
 621   TH2F* pZYCResidualsHist = static_cast<TH2F*>(pList->At(1));
 622   TH2F* pLPAResidualsHist = static_cast<TH2F*>(pList->At(2));
 623   TH2F* pLPCResidualsHist = static_cast<TH2F*>(pList->At(3));
 624   TH2F* pPhiYAResidualsHist = static_cast<TH2F*>(pList->At(4));
 625   TH2F* pPhiYCResidualsHist = static_cast<TH2F*>(pList->At(5));
 626   TH2F* pPhiZAResidualsHist = static_cast<TH2F*>(pList->At(6));
 627   TH2F* pPhiZCResidualsHist = static_cast<TH2F*>(pList->At(7));
 628   TH2F* pPtYAResidualsHist = static_cast<TH2F*>(pList->At(8));
 629   TH2F* pPtYCResidualsHist = static_cast<TH2F*>(pList->At(9));
 630   TH2F* pPtZAResidualsHist = static_cast<TH2F*>(pList->At(10));
 631   TH2F* pPtZCResidualsHist = static_cast<TH2F*>(pList->At(11));
 632   TH2F* pLowPtYAResidualsHist = static_cast<TH2F*>(pList->At(12));
 633   TH2F* pLowPtYCResidualsHist = static_cast<TH2F*>(pList->At(13));
 634   TH2F* pLowPtZAResidualsHist = static_cast<TH2F*>(pList->At(14));
 635   TH2F* pLowPtZCResidualsHist = static_cast<TH2F*>(pList->At(15));
 636
 637   if (paramsITS->GetZ() > 0.0)
 638   {
 639     pZYAResidualsHist->Fill(resz,resy);
 640     pLPAResidualsHist->Fill(restgl,ressnp);
 641     pPhiYAResidualsHist->Fill(paramsTPC->Phi(),resy);
 642     pPhiZAResidualsHist->Fill(paramsTPC->Phi(),resz);
 643     pPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
 644     pPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
 645     pLowPtYAResidualsHist->Fill(paramsTPC->Pt(),resy);
 646     pLowPtZAResidualsHist->Fill(paramsTPC->Pt(),resz);
 647   }
 648   else
 649   {
 650     pZYCResidualsHist->Fill(resz,resy);
 651     pLPCResidualsHist->Fill(restgl,ressnp);
 652     pPhiYCResidualsHist->Fill(paramsTPC->Phi(),resy);
 653     pPhiZCResidualsHist->Fill(paramsTPC->Phi(),resz);
 654     pPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
 655     pPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);
 656     pLowPtYCResidualsHist->Fill(paramsTPC->Pt(),resy);
 657     pLowPtZCResidualsHist->Fill(paramsTPC->Pt(),resz);
 658   }
 659 }