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