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   // load the geometry
 134   if(!gGeoManager) {
 135     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 136     if(!gGeoManager) {
 137       printf("AliTrackMatchingTPCITSCosmics::CreateOutputObjects(): no geometry loaded \n");
 138       return;
 139     }
 140   }
 141
 142   // Several histograms are more conveniently managed in a TList
 143   fList = new TList();
 144   fList->SetOwner();
 145
 146   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);
 147   fList->Add(fHistEvCount);
 148
 149   fntTrks = new TNtuple("fntTrks","TPC-ITS matching","ptTPC:nClsTPC:nClsITSSA:nClsITSMI:SSD1:SSD2:phi:z:dx:dy:dz:drphi:dphi:dtgl");
 150   fList->Add(fntTrks);
 151   fntPairs = new TNtuple("fntPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:nClsITSMIin:nClsITSMIout:dxyITSSA:dzITSSA:dxyITSMI:dzITSMI:ptITSSAin:ptITSSAout:ptITSMIin:ptITSMIout:sigmad0MIin:sigmad0MIout");
 152   fList->Add(fntPairs);
 153   fntITSPairs = new TNtuple("fntITSPairs","pairs at vertex","nClsITSSAin:nClsITSSAout:dxyITSSA:dzITSSA:ptITSSAin:ptITSSAout:d0mu");
 154   fList->Add(fntITSPairs);
 155
 156   return;
 157 }
 158
 159 //________________________________________________________________________
 160 void AliTrackMatchingTPCITSCosmics::Exec(Option_t */*option*/)
 161 {
 162   // Execute analysis for current event:
 163   // write ITS AliTrackPoints for selected tracks to fspTree
 164
 165   // check the geometry
 166   if(!gGeoManager) {
 167     printf("AliTrackMatchingTPCITSCosmics::Exec(): no geometry loaded \n");
 168     return;
 169   }
 170
 171   if(!fESD) {
 172     printf("AliTrackMatchingTPCITSCosmics::Exec(): no ESD \n");
 173     return;
 174   }
 175
 176
 177   fHistEvCount->Fill(0);
 178   PostData(0,fList);
 179
 180   // check if event is triggered by HLT
 181   Bool_t hltTrigg=kFALSE;
 182   if(fReadHLTESD) {
 183     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 184     AliESDEvent *hltESD = esdH->GetHLTEvent();
 185     if(!hltESD) {
 186       printf("AliTrackMatchingTPCITSCosmics::Exec(): no HLT ESD \n");
 187       return;
 188     }
 189     if(hltESD->IsHLTTriggerFired()) {fHistEvCount->Fill(5);hltTrigg=kTRUE;}
 190   }
 191
 192
 193   TString triggeredClass = fESD->GetFiredTriggerClasses();
 194   if(triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) {
 195     fHistEvCount->Fill(1);
 196   } else {
 197     if(fOnlySPDFO) {  PostData(0,fList); return; }
 198   }
 199
 200   Int_t ntracks = fESD->GetNumberOfTracks();
 201   printf("CONTR: %d\n",fESD->GetVertex()->GetNContributors());
 202   if(fESD->GetVertex()->GetNContributors()>=0) {
 203     fHistEvCount->Fill(2);
 204   }
 205
 206   Int_t nTrksITSSA=0,nTrksTPC=0,nTrksTPCITS=0;
 207   Int_t idxITSSA[10000],idxTPC[10000];
 208   Int_t idxITSSAin=-1,idxITSSAout=-1;
 209   Int_t idxTPCin=-1,idxTPCout=-1;
 210   for(Int_t itr=0; itr<ntracks; itr++) {
 211     AliESDtrack *t = fESD->GetTrack(itr);
 212     Int_t nClsTPC = t->GetNcls(1);
 213     Int_t nClsITS = t->GetNcls(0);
 214     if(nClsITS>=2 && nClsTPC==0) { // ITS SA
 215       idxITSSA[nTrksITSSA]=itr;
 216       nTrksITSSA++;
 217       //printf("%d\n",itr);
 218     }
 219     if(nClsTPC>=50) {              // TPC
 220       idxTPC[nTrksTPC]=itr;
 221       nTrksTPC++;
 222     }
 223     if(nClsTPC>=50 && nClsITS>=2) {// TPC+ITS
 224       nTrksTPCITS++;
 225       /*printf("  --> TPC+ITS: %d + %d (%d,%d,%d,%d,%d,%d) pt = %f\n",nClsTPC,nClsITS,
 226         (Int_t)(t->HasPointOnITSLayer(0)),
 227         (Int_t)(t->HasPointOnITSLayer(1)),
 228         (Int_t)(t->HasPointOnITSLayer(2)),
 229         (Int_t)(t->HasPointOnITSLayer(3)),
 230         (Int_t)(t->HasPointOnITSLayer(4)),
 231         (Int_t)(t->HasPointOnITSLayer(5)),
 232         t->Pt());*/
 233     }
 234   }
 235
 236   //printf("nTrksTPC %d  nTrksITSSA %d\n",nTrksTPC,nTrksITSSA);
 237   if(nTrksITSSA>0) fHistEvCount->Fill(4);
 238   if(hltTrigg && nTrksITSSA>0) fHistEvCount->Fill(7);
 239   if(nTrksTPC>0)   fHistEvCount->Fill(3);
 240   if(hltTrigg && nTrksTPC>0)   fHistEvCount->Fill(6);
 241   if(nTrksTPCITS>0)   fHistEvCount->Fill(8);
 242   if(hltTrigg && nTrksTPCITS>0)   fHistEvCount->Fill(9);
 243
 244   if(nTrksTPC>2 || nTrksTPC==0 || nTrksITSSA>2 || nTrksITSSA==0) { PostData(0,fList); return; }
 245
 246   for(Int_t itr=0; itr<nTrksTPC; itr++) {
 247     AliESDtrack *t = fESD->GetTrack(idxTPC[itr]);
 248     if(t->Py()>0) idxTPCin=idxTPC[itr];
 249     if(t->Py()<0) idxTPCout=idxTPC[itr];
 250   }
 251
 252
 253   for(Int_t itr=0; itr<nTrksITSSA; itr++) {
 254     AliESDtrack *t = fESD->GetTrack(idxITSSA[itr]);
 255     if(t->Py()>0) idxITSSAin=idxITSSA[itr];
 256     if(t->Py()<0) idxITSSAout=idxITSSA[itr];
 257   }
 258
 259   Double_t xyzITS[3],xyzTPC[3];
 260
 261   // analysis for inward track
 262   if(idxITSSAin>-1 && idxTPCin>-1) {
 263     //printf(" %d\n",idxITSSAin);
 264     AliESDtrack *tITS = fESD->GetTrack(idxITSSAin);
 265     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 266     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 267     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 268     AliESDtrack *tTPC = fESD->GetTrack(idxTPCin);
 269     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 270     innerTPC->Rotate(tITS->GetAlpha());
 271     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 272     fntTrks->Fill(innerTPC->Pt(),
 273                  tTPC->GetNcls(1),
 274                  tITS->GetNcls(0),
 275                  tTPC->GetNcls(0),
 276                  tITS->HasPointOnITSLayer(4),
 277                  tITS->HasPointOnITSLayer(5),
 278                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 279                  xyzITS[2],
 280                  xyzITS[0]-xyzTPC[0],
 281                  xyzITS[1]-xyzTPC[1],
 282                  xyzITS[2]-xyzTPC[2],
 283                  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])),
 284                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 285                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt());
 286     delete innerTPC; innerTPC=NULL;
 287   }
 288
 289   // analysis for outward track
 290   if(idxITSSAout>-1 && idxTPCout>-1) {
 291     //printf(" %d\n",idxITSSAout);
 292     AliESDtrack *tITS = fESD->GetTrack(idxITSSAout);
 293     const AliExternalTrackParam *outerITS = tITS->GetOuterParam();
 294     xyzITS[0]=-1000.;xyzITS[1]=-1000.;xyzITS[2]=-1000.;
 295     if(outerITS) outerITS->GetXYZAt(50.,fESD->GetMagneticField(),xyzITS);
 296     AliESDtrack *tTPC = fESD->GetTrack(idxTPCout);
 297     AliExternalTrackParam *innerTPC = new AliExternalTrackParam(*(tTPC->GetInnerParam()));
 298     innerTPC->Rotate(tITS->GetAlpha());
 299     innerTPC->GetXYZAt(50.,fESD->GetMagneticField(),xyzTPC);
 300     fntTrks->Fill(innerTPC->Pt(),
 301                  tTPC->GetNcls(1),
 302                  tITS->GetNcls(0),
 303                  tTPC->GetNcls(0),
 304                  tITS->HasPointOnITSLayer(4),
 305                  tITS->HasPointOnITSLayer(5),
 306                  TMath::ATan2(xyzITS[1],xyzITS[0]),
 307                  xyzITS[2],
 308                  xyzITS[0]-xyzTPC[0],
 309                  xyzITS[1]-xyzTPC[1],
 310                  xyzITS[2]-xyzTPC[2],
 311                  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])),
 312                  TMath::ATan2(tITS->Py(),tITS->Px())-TMath::ATan2(innerTPC->Py(),innerTPC->Px()),
 313                  tITS->Pz()/tITS->Pt()-innerTPC->Pz()/innerTPC->Pt()
 314                  );
 315     delete innerTPC; innerTPC=NULL;
 316   }
 317
 318   // inward-outward track ITS only
 319   if(idxITSSAin>-1 && idxITSSAout>-1) {
 320     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 321     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 322     if(tITSin->HasPointOnITSLayer(0) &&
 323        tITSin->HasPointOnITSLayer(1) &&
 324        tITSout->HasPointOnITSLayer(0) &&
 325        tITSout->HasPointOnITSLayer(1)) {
 326       Double_t alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 327       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 328       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 329       Float_t d0[2],z0[2];
 330       d0[0] = tITSin->GetY();
 331       z0[0] = tITSin->GetZ();
 332       d0[1] = tITSout->GetY();
 333       z0[1] = tITSout->GetZ();
 334       Float_t dxyITS = -(d0[0]-d0[1]);
 335       Float_t dzITS = z0[0]-z0[1];
 336       fntITSPairs->Fill(tITSin->GetNcls(0),tITSout->GetNcls(0),dxyITS,dzITS,tITSin->Pt(),tITSout->Pt(),TMath::Abs(d0[0]));
 337     }
 338   }
 339
 340   // inward-outward track TPC+ITS
 341   if(idxITSSAin>-1 && idxTPCin>-1 && idxITSSAout>-1 && idxTPCout>-1) {
 342     AliESDtrack *tTPCin = fESD->GetTrack(idxTPCin);
 343     AliESDtrack *tITSin = fESD->GetTrack(idxITSSAin);
 344     AliESDtrack *tTPCout = fESD->GetTrack(idxTPCout);
 345     AliESDtrack *tITSout = fESD->GetTrack(idxITSSAout);
 346     if(tTPCin->HasPointOnITSLayer(0) &&
 347        tTPCin->HasPointOnITSLayer(1) &&
 348        tITSin->HasPointOnITSLayer(0) &&
 349        tITSin->HasPointOnITSLayer(1) &&
 350        tTPCout->HasPointOnITSLayer(0) &&
 351        tTPCout->HasPointOnITSLayer(1) &&
 352        tITSout->HasPointOnITSLayer(0) &&
 353        tITSout->HasPointOnITSLayer(1)) {
 354       Double_t alpha = TMath::ATan2(tTPCin->Py(),tTPCin->Px());
 355       tTPCin->Propagate(alpha,0.,fESD->GetMagneticField());
 356       tTPCout->Propagate(alpha,0.,fESD->GetMagneticField());
 357       Float_t d0[2],z0[2];
 358       d0[0] = tTPCin->GetY();
 359       z0[0] = tTPCin->GetZ();
 360       d0[1] = tTPCout->GetY();
 361       z0[1] = tTPCout->GetZ();
 362       Float_t dxyTPC = -(d0[0]-d0[1]);
 363       Float_t dzTPC = z0[0]-z0[1];
 364       alpha = TMath::ATan2(tITSin->Py(),tITSin->Px());
 365       tITSin->Propagate(alpha,0.,fESD->GetMagneticField());
 366       tITSout->Propagate(alpha,0.,fESD->GetMagneticField());
 367       d0[0] = tITSin->GetY();
 368       z0[0] = tITSin->GetZ();
 369       d0[1] = tITSout->GetY();
 370       z0[1] = tITSout->GetZ();
 371       Float_t dxyITS = -(d0[0]-d0[1]);
 372       Float_t dzITS = z0[0]-z0[1];
 373       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());
 374     }
 375   }
 376
 377   PostData(0,fList);
 378
 379   return;
 380 }
 381
 382 //________________________________________________________________________
 383 void AliTrackMatchingTPCITSCosmics::Terminate(Option_t */*option*/)
 384 {
 385   // Terminate analysis
 386   //
 387   AliDebug(2,"AliTrackMatchingTPCITSCosmics: Terminate() \n");
 388
 389   fList = dynamic_cast<TList*> (GetOutputData(0));
 390   if (!fList) {
 391     printf("ERROR: fList not available\n");
 392     return;
 393   }
 394
 395   return;
 396 }
 397