PWG1/AliTrackMatchingTPCITSCosmics.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTask to study the TPC-ITS track matching
  19 //
  20 // Author: A.Dainese, andrea.dainese@pd.infn.it
  21 /////////////////////////////////////////////////////////////
  22
  23 #include <TTree.h>
  24 #include <TChain.h>
  25 #include <TNtuple.h>
  26 #include <TList.h>
  27 #include <TH1F.h>
  28 #include <TH2F.h>
  29 #include <TVector3.h>
  30 #include <TGeoManager.h>
  31
  32 #include "AliLog.h"
  33 #include "AliGeomManager.h"
  34 #include "AliTrackPointArray.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliESDVertex.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDfriend.h"
  40 #include "AliAnalysisTask.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliExternalTrackParam.h"
  43 #include "AliTrackMatchingTPCITSCosmics.h"
  44
  45
  46 ClassImp(AliTrackMatchingTPCITSCosmics)
  47
  48
  49 //________________________________________________________________________
  50 AliTrackMatchingTPCITSCosmics::AliTrackMatchingTPCITSCosmics(const char *name):
  51 AliAnalysisTask(name,"task"),
  52 fOnlySPDFO(kFALSE),
  53 fReadHLTESD(kFALSE),
  54 fGeometryFileName("geometry.root"),
  55 fESD(0),
  56 fList(0),
  57 fHistEvCount(0),
  58 fntTrks(0),
  59 fntPairs(0),
  60 fntITSPairs(0)
  61 {
  62   // Constructor
  63
  64   // Define input and output slots here
  65   // Input slot #0 works with a TChain
  66   DefineInput(0, TChain::Class());
  67   // Output slot #0 writes into a TList
  68   DefineOutput(0,TList::Class());  //My private output
  69 }
  70
  71 //________________________________________________________________________
  72 AliTrackMatchingTPCITSCosmics::~AliTrackMatchingTPCITSCosmics()
  73 {
  74   // Destructor
  75   if (fList) {
  76     delete fList;
  77     fList = 0;
  78   }
  79   if (fHistEvCount) {
  80     delete fHistEvCount;
  81     fHistEvCount = 0;
  82   }
  83   if (fntTrks) {
  84     delete fntTrks;
  85     fntTrks = 0;
  86   }
  87   if (fntPairs) {
  88     delete fntPairs;
  89     fntPairs = 0;
  90   }
  91   if (fntITSPairs) {
  92     delete fntITSPairs;
  93     fntITSPairs = 0;
  94   }
  95 }
  96
  97 //________________________________________________________________________
  98 void AliTrackMatchingTPCITSCosmics::ConnectInputData(Option_t *)
  99 {
 100   // Connect ESD
 101   // Called once
 102
 103   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 104   if(!tree) {
 105     printf("ERROR: Could not read chain from input slot 0\n");
 106   } else {
 107     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 108
 109     if(!esdH) {
 110       printf("ERROR: Could not get ESDInputHandler\n");
 111     } else {
 112       fESD = esdH->GetEvent();
 113     }
 114   }
 115
 116   return;
 117 }
 118
 119 //________________________________________________________________________
 120 void AliTrackMatchingTPCITSCosmics::Init()
 121 {
 122   // Initialization
 123
 124   return;
 125 }
 126
 127 //________________________________________________________________________
 128 void AliTrackMatchingTPCITSCosmics::CreateOutputObjects()
 129 {
 130   // Create the output container
 131   //
 132
 133   // Several histograms are more conveniently managed in a TList
 134   fList = new TList();
 135   fList->SetOwner();
 136
 137   fHistEvCount = new TH1F("fHistEvCount","0: all, 1: SPDFO, 2: vtx, 3: TPCtrks, 4: ITStrks, 5: HLT, 6: HLT&TPCtrks, 7: HLT&ITStrks",11,-0.5,10.5);
 138   fList->Add(fHistEvCount);
 139
 140   fntTrks = new TNtuple("fntTrks","TPC-ITS matching","ptTPC:nClsTPC:nClsITSSA:nClsITSMI:SSD1:SSD2:phi:z:dx:dy:dz:drphi:dphi:dtgl");
 141   fList->Add(fntTrks);
 142   fntPairs = new TNtuple("fntPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:nClsITSMIin:nClsITSMIout:dxyITSSA:dzITSSA:dxyITSMI:dzITSMI:ptITSSAin:ptITSSAout:ptITSMIin:ptITSMIout:sigmad0MIin:sigmad0MIout");
 143   fList->Add(fntPairs);
 144   fntITSPairs = new TNtuple("fntITSPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:dxyITSSA:dzITSSA:ptITSSAin:ptITSSAout:d0mu");
 145   fList->Add(fntITSPairs);
 146
 147   return;
 148 }
 149
 150 //________________________________________________________________________
 151 void AliTrackMatchingTPCITSCosmics::Exec(Option_t */*option*/)
 152 {
 153   // Execute analysis for current event:
 154   // write ITS AliTrackPoints for selected tracks to fspTree
 155
 156   // load the geometry
 157   if(!gGeoManager) {
 158     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 159     if(!gGeoManager) {
 160       printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 161       return;
 162     }
 163   }
 164
 165   if(!fESD) {
 166     printf("AliTrackMatchingTPCITSCosmics::Exec(): no ESD \n");
 167     return;
 168   }
 169
 170
 171   fHistEvCount->Fill(0);
 172   PostData(0,fList);
 173
 174   // check if event is triggered by HLT
 175   Bool_t hltTrigg=kFALSE;
 176   if(fReadHLTESD) {
 177     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 178     AliESDEvent *hltESD = esdH->GetHLTEvent();
 179     if(!hltESD) {
 180       printf("AliTrackMatchingTPCITSCosmics::Exec(): no HLT ESD \n");
 181       return;
 182     }
 183     if(hltESD->IsHLTTriggerFired()) {fHistEvCount->Fill(5);hltTrigg=kTRUE;}
 184   }
 185
 186
 187   TString triggeredClass = fESD->GetFiredTriggerClasses();
 188   if(triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) {
 189     fHistEvCount->Fill(1);
 190   } else {
 191     if(fOnlySPDFO) {  PostData(0,fList); return; }
 192   }
 193
 194   Int_t ntracks = fESD->GetNumberOfTracks();
 195   printf("CONTR: %d\n",fESD->GetVertex()->GetNContributors());
 196   if(fESD->GetVertex()->GetNContributors()>=0) {
 197     fHistEvCount->Fill(2);
 198   }
 199
 200   Int_t nTrksITSSA=0,nTrksTPC=0,nTrksTPCITS=0;
 201   Int_t idxITSSA[10000],idxTPC[10000];
 202   Int_t idxITSSAin=-1,idxITSSAout=-1;
 203   Int_t idxTPCin=-1,idxTPCout=-1;
 204   for(Int_t itr=0; itr<ntracks; itr++) {
 205     AliESDtrack *t = fESD->GetTrack(itr);
 206     Int_t nClsTPC = t->GetNcls(1);
 207     Int_t nClsITS = t->GetNcls(0);
 208     if(nClsITS>=2 && nClsTPC==0) { // ITS SA
 209       idxITSSA[nTrksITSSA]=itr;
 210       nTrksITSSA++;
 211       //printf("%d\n",itr);
 212     }
 213     if(nClsTPC>=50) {              // TPC
 214       idxTPC[nTrksTPC]=itr;
 215       nTrksTPC++;
 216     }
 217     if(nClsTPC>=50 && nClsITS>=2) {// TPC+ITS
 218       nTrksTPCITS++;
 219       /*printf("  --> TPC+ITS: %d + %d (%d,%d,%d,%d,%d,%d) pt = %f\n",nClsTPC,nClsITS,
 220         (Int_t)(t->HasPointOnITSLayer(0)),
 221         (Int_t)(t->HasPointOnITSLayer(1)),
 222         (Int_t)(t->HasPointOnITSLayer(2)),
 223         (Int_t)(t->HasPointOnITSLayer(3)),
 224         (Int_t)(t->HasPointOnITSLayer(4)),
 225         (Int_t)(t->HasPointOnITSLayer(5)),
 226         t->Pt());*/
 227     }
 228   }
 229
 230   //printf("nTrksTPC %d  nTrksITSSA %d\n",nTrksTPC,nTrksITSSA);
 231   if(nTrksITSSA>0) fHistEvCount->Fill(4);
 232   if(hltTrigg && nTrksITSSA>0) fHistEvCount->Fill(7);
 233   if(nTrksTPC>0)   fHistEvCount->Fill(3);
 234   if(hltTrigg && nTrksTPC>0)   fHistEvCount->Fill(6);
 235   if(nTrksTPCITS>0)   fHistEvCount->Fill(8);
 236   if(hltTrigg && nTrksTPCITS>0)   fHistEvCount->Fill(9);
 237
 238   if(nTrksTPC>2 || nTrksTPC==0 || nTrksITSSA>2 || nTrksITSSA==0) { PostData(0,fList); return; }
 239
 240   for(Int_t itr=0; itr<nTrksTPC; itr++) {
 241     AliESDtrack *t = fESD->GetTrack(idxTPC[itr]);
 242     if(t->Py()>0) idxTPCin=idxTPC[itr];
 243     if(t->Py()<0) idxTPCout=idxTPC[itr];
 244   }
 245
 246
 247   for(Int_t itr=0; itr<nTrksITSSA; itr++) {
 248     AliESDtrack *t = fESD->GetTrack(idxITSSA[itr]);
 249     if(t->Py()>0) idxITSSAin=idxITSSA[itr];
 250     if(t->Py()<0) idxITSSAout=idxITSSA[itr];
 251   }
 252
 253   Double_t xyzITS[3],xyzTPC[3];
 254
 255   // analysis for inward track
 256   if(idxITSSAin>-1 && idxTPCin>-1) {
 257     //printf(" %d\n",idxITSSAin);
 258     AliESDtrack *tITS = fESD->GetTrack(idxITSSAin);
 259     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 260     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 261     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 262     AliESDtrack *tTPC = fESD->GetTrack(idxTPCin);
 263     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 264     innerTPC->Rotate(tITS->GetAlpha());
 265     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 266     fntTrks->Fill(innerTPC->Pt(),
 267                  tTPC->GetNcls(1),
 268                  tITS->GetNcls(0),
 269                  tTPC->GetNcls(0),
 270                  tITS->HasPointOnITSLayer(4),
 271                  tITS->HasPointOnITSLayer(5),
 272                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 273                  xyzITS[2],
 274                  xyzITS[0]-xyzTPC[0],
 275                  xyzITS[1]-xyzTPC[1],
 276                  xyzITS[2]-xyzTPC[2],
 277                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
 278                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 279                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt());
 280     delete innerTPC; innerTPC=NULL;
 281   }
 282
 283   // analysis for outward track
 284   if(idxITSSAout>-1 && idxTPCout>-1) {
 285     //printf(" %d\n",idxITSSAout);
 286     AliESDtrack *tITS = fESD->GetTrack(idxITSSAout);
 287     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 288     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 289     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 290     AliESDtrack *tTPC = fESD->GetTrack(idxTPCout);
 291     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 292     innerTPC->Rotate(tITS->GetAlpha());
 293     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 294     fntTrks->Fill(innerTPC->Pt(),
 295                  tTPC->GetNcls(1),
 296                  tITS->GetNcls(0),
 297                  tTPC->GetNcls(0),
 298                  tITS->HasPointOnITSLayer(4),
 299                  tITS->HasPointOnITSLayer(5),
 300                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 301                  xyzITS[2],
 302                  xyzITS[0]-xyzTPC[0],
 303                  xyzITS[1]-xyzTPC[1],
 304                  xyzITS[2]-xyzTPC[2],
 305                  0.5*(TMath::Sqrt(xyzITS[0]*xyzITS[0]+xyzITS[1]*xyzITS[1])+TMath::Sqrt(xyzTPC[0]*xyzTPC[0]+xyzTPC[1]*xyzTPC[1]))*(TMath::ATan2(xyzITS[1],xyzITS[0])-TMath::ATan2(xyzTPC[1],xyzTPC[0])),
 306                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 307                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt()
 308                  );
 309     delete innerTPC; innerTPC=NULL;
 310   }
 311
 312   // inward-outward track ITS only
 313   if(idxITSSAin>-1 && idxITSSAout>-1) {
 314     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 315     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 316     if(tITSin->HasPointOnITSLayer(0) &&
 317        tITSin->HasPointOnITSLayer(1) &&
 318        tITSout->HasPointOnITSLayer(0) &&
 319        tITSout->HasPointOnITSLayer(1)) {
 320       Double_t alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 321       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 322       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 323       Float_t d0[2],z0[2];
 324       d0[0] = tITSin->GetY();
 325       z0[0] = tITSin->GetZ();
 326       d0[1] = tITSout->GetY();
 327       z0[1] = tITSout->GetZ();
 328       Float_t dxyITS = -(d0[0]-d0[1]);
 329       Float_t dzITS = z0[0]-z0[1];
 330       fntITSPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),dxyITS,dzITS,tITSin->Pt(),tITSout->Pt(),TMath::Abs(d0[0]));
 331     }
 332   }
 333
 334   // inward-outward track TPC+ITS
 335   if(idxITSSAin>-1 && idxTPCin>-1 && idxITSSAout>-1 && idxTPCout>-1) {
 336     AliESDtrack *tTPCin = fESD->GetTrack(idxTPCin);
 337     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 338     AliESDtrack *tTPCout = fESD->GetTrack(idxTPCout);
 339     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 340     if(tTPCin->HasPointOnITSLayer(0) &&
 341        tTPCin->HasPointOnITSLayer(1) &&
 342        tITSin->HasPointOnITSLayer(0) &&
 343        tITSin->HasPointOnITSLayer(1) &&
 344        tTPCout->HasPointOnITSLayer(0) &&
 345        tTPCout->HasPointOnITSLayer(1) &&
 346        tITSout->HasPointOnITSLayer(0) &&
 347        tITSout->HasPointOnITSLayer(1)) {
 348       Double_t alpha = TMath::ATan2(tTPCin->Py(),tTPCin->Px());
 349       tTPCin->Propagate(alpha,0.,fESD->GetMagneticField());
 350       tTPCout->Propagate(alpha,0.,fESD->GetMagneticField());
 351       Float_t d0[2],z0[2];
 352       d0[0] = tTPCin->GetY();
 353       z0[0] = tTPCin->GetZ();
 354       d0[1] = tTPCout->GetY();
 355       z0[1] = tTPCout->GetZ();
 356       Float_t dxyTPC = -(d0[0]-d0[1]);
 357       Float_t dzTPC = z0[0]-z0[1];
 358       alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 359       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 360       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 361       d0[0] = tITSin->GetY();
 362       z0[0] = tITSin->GetZ();
 363       d0[1] = tITSout->GetY();
 364       z0[1] = tITSout->GetZ();
 365       Float_t dxyITS = -(d0[0]-d0[1]);
 366       Float_t dzITS = z0[0]-z0[1];
 367       fntPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),tTPCin->GetNcls(0),tTPCout->GetNcls(0),dxyITS,dzITS,dxyTPC,dzTPC,tITSin->Pt(),tITSout->Pt(),tTPCin->Pt(),tTPCout->Pt(),tTPCin->GetSigmaY2(),tTPCout->GetSigmaY2());
 368     }
 369   }
 370
 371   PostData(0,fList);
 372
 373   return;
 374 }
 375
 376 //________________________________________________________________________
 377 void AliTrackMatchingTPCITSCosmics::Terminate(Option_t */*option*/)
 378 {
 379   // Terminate analysis
 380   //
 381   AliDebug(2,"AliTrackMatchingTPCITSCosmics: Terminate() \n");
 382
 383   fList = dynamic_cast<TList*> (GetOutputData(0));
 384   if (!fList) {
 385     printf("ERROR: fList not available\n");
 386     return;
 387   }
 388
 389   return;
 390 }
 391