PWG1/AliAlignmentDataFilterITS.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 extract from ESD tracks the AliTrackPointArrays
  19 // with ITS points for selected tracks. This are the input data for alignment
  20 //
  21 // Author: A.Dainese, andrea.dainese@pd.infn.it
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TTree.h>
  25 #include <TChain.h>
  26 #include <TNtuple.h>
  27 #include <TList.h>
  28 #include <TH1F.h>
  29 #include <TH2F.h>
  30 #include <TVector3.h>
  31 #include <TGeoManager.h>
  32
  33 #include "AliLog.h"
  34 #include "AliGeomManager.h"
  35 #include "AliITSReconstructor.h"
  36 #include "AliITSgeomTGeo.h"
  37 #include "AliTrackPointArray.h"
  38 #include "AliESDInputHandler.h"
  39 #include "AliESDVertex.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDfriend.h"
  43 #include "AliAnalysisTask.h"
  44 #include "AliAnalysisManager.h"
  45 #include "AliAlignmentDataFilterITS.h"
  46
  47
  48 ClassImp(AliAlignmentDataFilterITS)
  49
  50
  51 //________________________________________________________________________
  52 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
  53 AliAnalysisTask(name,"task"),
  54 fOnlySPDFO(kFALSE),
  55 fGeometryFileName("geometry.root"),
  56 fITSRecoParam(0),
  57 fESD(0),
  58 fESDfriend(0),
  59 fListOfHistos(0),
  60 fspTree(0),
  61 fHistNevents(0),
  62 fHistNpoints(0),
  63 fHistPt(0),
  64 fHistLayer0(0),
  65 fHistLayer1(0),
  66 fHistLayer2(0),
  67 fHistLayer3(0),
  68 fHistLayer4(0),
  69 fHistLayer5(0),
  70 fntExtra(0),
  71 fntCosmicMatching(0)
  72 {
  73   // Constructor
  74
  75   // Define input and output slots here
  76   // Input slot #0 works with a TChain
  77   DefineInput(0, TChain::Class());
  78   // Output slot #0 writes into a TTree
  79   DefineOutput(0,TTree::Class());  //My private output
  80   // Output slot #1 writes into a TList
  81   DefineOutput(1,TList::Class());  //My private output
  82 }
  83
  84 //________________________________________________________________________
  85 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
  86 {
  87   // Destructor
  88   if (fListOfHistos) {
  89     delete fListOfHistos;
  90     fListOfHistos = 0;
  91   }
  92   if (fspTree) {
  93     delete fspTree;
  94     fspTree = 0;
  95   }
  96   if (fHistNevents) {
  97     delete fHistNevents;
  98     fHistNevents = 0;
  99   }
 100   if (fHistNpoints) {
 101     delete fHistNpoints;
 102     fHistNpoints = 0;
 103   }
 104   if (fHistPt) {
 105     delete fHistPt;
 106     fHistPt = 0;
 107   }
 108   if (fHistLayer0) {
 109     delete fHistLayer0;
 110     fHistLayer0 = 0;
 111   }
 112   if (fHistLayer1) {
 113     delete fHistLayer1;
 114     fHistLayer1 = 0;
 115   }
 116   if (fHistLayer2) {
 117     delete fHistLayer2;
 118     fHistLayer2 = 0;
 119   }
 120   if (fHistLayer3) {
 121     delete fHistLayer3;
 122     fHistLayer3 = 0;
 123   }
 124   if (fHistLayer4) {
 125     delete fHistLayer4;
 126     fHistLayer4 = 0;
 127   }
 128   if (fHistLayer5) {
 129     delete fHistLayer5;
 130     fHistLayer5 = 0;
 131   }
 132   if (fntExtra) {
 133     delete fntExtra;
 134     fntExtra = 0;
 135   }
 136   if (fntCosmicMatching) {
 137     delete fntCosmicMatching;
 138     fntCosmicMatching = 0;
 139   }
 140 }
 141
 142 //________________________________________________________________________
 143 void AliAlignmentDataFilterITS::ConnectInputData(Option_t *)
 144 {
 145   // Connect ESD
 146   // Called once
 147
 148   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 149   if(!tree) {
 150     printf("ERROR: Could not read chain from input slot 0\n");
 151   } else {
 152     // Disable all branches and enable only the needed ones
 153
 154     tree->SetBranchStatus("fTriggerMask", 1);
 155     tree->SetBranchStatus("fSPDVertex*", 1);
 156
 157     tree->SetBranchStatus("ESDfriend*", 1);
 158     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
 159
 160     tree->SetBranchStatus("fSPDMult*", 1);
 161
 162     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 163
 164     if(!esdH) {
 165       printf("ERROR: Could not get ESDInputHandler\n");
 166     } else {
 167       fESD = esdH->GetEvent();
 168     }
 169   }
 170
 171   return;
 172 }
 173
 174 //________________________________________________________________________
 175 void AliAlignmentDataFilterITS::Init()
 176 {
 177   // Initialization
 178
 179   return;
 180 }
 181
 182 //________________________________________________________________________
 183 void AliAlignmentDataFilterITS::CreateOutputObjects()
 184 {
 185   // Create the output container
 186   //
 187
 188   // load the geometry
 189   if(!gGeoManager) {
 190     printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data());
 191     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 192     if(!gGeoManager) {
 193       printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n");
 194       return;
 195     }
 196   }
 197
 198   // Several histograms are more conveniently managed in a TList
 199   fListOfHistos = new TList();
 200   fListOfHistos->SetOwner();
 201
 202   fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5);
 203   fHistNevents->Sumw2();
 204   fHistNevents->SetMinimum(0);
 205   fListOfHistos->Add(fHistNevents);
 206
 207   fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
 208   fHistNpoints->Sumw2();
 209   fHistNpoints->SetMinimum(0);
 210   fListOfHistos->Add(fHistNpoints);
 211
 212   fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
 213   fHistPt->Sumw2();
 214   fHistPt->SetMinimum(0);
 215   fListOfHistos->Add(fHistPt);
 216
 217
 218   Float_t zmax=14.;
 219   Int_t nbinsphi=20,nbinsz=4;
 220   fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 221   fListOfHistos->Add(fHistLayer0);
 222   zmax=14.;
 223   nbinsphi=40;nbinsz=4;
 224   fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 225   fListOfHistos->Add(fHistLayer1);
 226   zmax=22.;
 227   nbinsphi=14;nbinsz=6;
 228   fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 229   fListOfHistos->Add(fHistLayer2);
 230   zmax=29.5;
 231   nbinsphi=22;nbinsz=8;
 232   fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 233   fListOfHistos->Add(fHistLayer3);
 234   zmax=45.;
 235   nbinsphi=34;nbinsz=23;
 236   fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 237   fListOfHistos->Add(fHistLayer4);
 238   zmax=51.;
 239   nbinsphi=38;nbinsz=26;
 240   fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 241   fListOfHistos->Add(fHistLayer5);
 242
 243
 244   fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu");
 245   fListOfHistos->Add(fntExtra);
 246
 247   fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
 248   fListOfHistos->Add(fntCosmicMatching);
 249
 250   fspTree = new TTree("spTree","Tree with ITS track points");
 251   const AliTrackPointArray *array = 0;
 252   Float_t curv,curverr;
 253   fspTree->Branch("SP","AliTrackPointArray",&array);
 254   fspTree->Branch("curv",&curv);
 255   fspTree->Branch("curverr",&curverr);
 256
 257   return;
 258 }
 259
 260 //________________________________________________________________________
 261 void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
 262 {
 263   // Execute analysis for current event:
 264   // write ITS AliTrackPoints for selected tracks to fspTree
 265
 266   // check the geometry
 267   if(!gGeoManager) {
 268     printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 269     return;
 270   }
 271
 272   // check if we have AliITSRecoParam
 273   if(!GetRecoParam()) {
 274     if(!fITSRecoParam) {
 275       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
 276       return;
 277     }
 278   }
 279
 280
 281   if(!fESD) {
 282     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
 283     return;
 284   }
 285   if(!fESDfriend) {
 286     printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
 287     return;
 288   }
 289   // attach ESDfriend
 290   fESD->SetESDfriend(fESDfriend);
 291
 292   // Post the data for slot 0
 293   fHistNevents->Fill(0);
 294
 295   // Process event as Cosmic or Collision
 296   //if(esd->GetEventType()== ???? ) {
 297   printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
 298   if(GetRecoParam()->GetAlignFilterCosmics()) {
 299     FilterCosmic(fESD);
 300   } else {
 301     FilterCollision(fESD);
 302   }
 303
 304   PostData(1,fListOfHistos);
 305
 306   return;
 307 }
 308
 309 //________________________________________________________________________
 310 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
 311 {
 312   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 313   // of top and bottom track, merge the two tracks, if requested)
 314   //
 315
 316   // Set branch addresses for space points tree
 317   AliTrackPointArray *arrayForTree=0;
 318   Float_t curv,curverr;
 319   fspTree->SetBranchAddress("SP",&arrayForTree);
 320   fspTree->SetBranchAddress("curv",&curv);
 321   fspTree->SetBranchAddress("curverr",&curverr);
 322
 323   TString triggeredClass = fESD->GetFiredTriggerClasses();
 324   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
 325
 326
 327   Int_t ntracks = esd->GetNumberOfTracks();
 328   if(ntracks<2) return;
 329
 330   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
 331
 332   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 333
 334   Int_t *goodtracksArray = new Int_t[ntracks];
 335   Float_t *phiArray = new Float_t[ntracks];
 336   Float_t *thetaArray = new Float_t[ntracks];
 337   Int_t *nclsArray = new Int_t[ntracks];
 338   Int_t ngt=0;
 339   Int_t itrack=0;
 340   for (itrack=0; itrack < ntracks; itrack++) {
 341     AliESDtrack *track = esd->GetTrack(itrack);
 342     if (!track) continue;
 343
 344
 345     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 346
 347     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 348     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 349
 350     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 351     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 352
 353     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 354        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 355
 356     goodtracksArray[ngt] = itrack;
 357     phiArray[ngt]        = phi;
 358     thetaArray[ngt]      = theta;
 359     nclsArray[ngt]       = track->GetNcls(0);
 360     ngt++;
 361   }
 362
 363   if(ngt<2) {
 364     delete [] goodtracksArray; goodtracksArray=0;
 365     delete [] phiArray; phiArray=0;
 366     delete [] thetaArray; thetaArray=0;
 367     delete [] nclsArray; nclsArray=0;
 368     return;
 369   }
 370
 371   // check matching of the two tracks from the muon
 372   Float_t min = 10000000.;
 373   Int_t maxCls = 0;
 374   Int_t good1 = -1, good2 = -1;
 375   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 376     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 377       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 378       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 379       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 380       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 381       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
 382         maxCls = nclsArray[itr1]+nclsArray[itr2];
 383         min = deltatheta+deltaphi;
 384         good1 = goodtracksArray[itr1];
 385         good2 = goodtracksArray[itr2];
 386       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
 387         if(deltatheta+deltaphi < min) {
 388           min = deltatheta+deltaphi;
 389           good1 = goodtracksArray[itr1];
 390           good2 = goodtracksArray[itr2];
 391         }
 392       }
 393     }
 394   }
 395
 396   delete [] goodtracksArray; goodtracksArray=0;
 397   delete [] phiArray; phiArray=0;
 398   delete [] thetaArray; thetaArray=0;
 399   delete [] nclsArray; nclsArray=0;
 400
 401   if(good1<0) return;
 402   AliDebug(2,"ok track matching");
 403
 404   // track1 will be the inward track (top)
 405   // track2 the outward (bottom)
 406   AliESDtrack *track1=0;
 407   AliESDtrack *track2=0;
 408   AliESDtrack *track = esd->GetTrack(good1);
 409   if(track->Py()>0) {
 410     track1 = esd->GetTrack(good1);
 411     track2 = esd->GetTrack(good2);
 412   } else {
 413     track1 = esd->GetTrack(good2);
 414     track2 = esd->GetTrack(good1);
 415   }
 416
 417   AliTrackPoint point;
 418   const AliTrackPointArray *array=0;
 419   Int_t ipt,volId,modId,layerId,lay,lad,det;
 420   Int_t jpt=0;
 421   Bool_t layerOK[6][2];
 422   Int_t nclsTrk[2]={0,0};
 423
 424   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
 425
 426   for(itrack=0; itrack<2; itrack++) {
 427     if(itrack==0) {
 428       track = track1;
 429     } else {
 430       track = track2;
 431     }
 432     array = track->GetTrackPointArray();
 433     if(!array) {
 434       AliWarning("No tracks points avaialble");
 435       continue;
 436     }
 437     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 438       array->GetPoint(point,ipt);
 439       volId = point.GetVolumeID();
 440       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 441       AliDebug(2,Form("%d %d\n",ipt,layerId-1));
 442       if(point.IsExtra() &&
 443          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 444       if(layerId>6) continue;
 445       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
 446       // check minAngleWrtITSModulePlanes
 447       Double_t p[3]; track->GetDirection(p);
 448       TVector3 pvec(p[0],p[1],p[2]);
 449       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 450       TVector3 normvec(rot[1],rot[4],rot[7]);
 451       Double_t angle = pvec.Angle(normvec);
 452       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 453       angle = 0.5*TMath::Pi()-angle;
 454       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
 455       layerOK[layerId-1][itrack]=kTRUE;
 456       jpt++;
 457       nclsTrk[itrack]++;
 458     }
 459   }
 460   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
 461
 462   // read ITS cluster maps
 463   Int_t map1[6],map2[6];
 464   for(Int_t ilay=0;ilay<6;ilay++) {
 465     map1[ilay]=0; map2[ilay]=0;
 466     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
 467     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
 468   }
 469   AliDebug(2,Form("ITS map 1: %d %d %d %d %d %d pt %f\n",map1[0],map1[1],map1[2],map1[3],map1[4],map1[5],track1->Pt()));
 470   AliDebug(2,Form("ITS map 2: %d %d %d %d %d %d pt %f\n",map2[0],map2[1],map2[2],map2[3],map2[4],map2[5],track2->Pt()));
 471   Int_t idx1[12],idx2[12];
 472   track1->GetITSclusters(idx1);
 473   track2->GetITSclusters(idx2);
 474   AliDebug(2,Form("cls idx 1 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx1[0],idx1[1],idx1[2],idx1[3],idx1[4],idx1[5],idx1[6],idx1[7],idx1[8],idx1[9],idx1[10],idx1[11]));
 475   AliDebug(2,Form("cls idx 2 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx2[0],idx2[1],idx2[2],idx2[3],idx2[4],idx2[5],idx2[6],idx2[7],idx2[8],idx2[9],idx2[10],idx2[11]));
 476
 477
 478   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
 479   AliDebug(2,Form(" Total points %d, accepted\n",jpt));
 480   fHistNpoints->Fill(jpt);
 481   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
 482
 483   Float_t d0z0mu[2];
 484   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
 485   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
 486
 487   Float_t dzOverlap[2];
 488   Float_t curvArray[2],curverrArray[2];
 489   Double_t globExtra[3],locExtra[3];
 490   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks())
 491     arrayForTree = new AliTrackPointArray(jpt);
 492
 493   jpt=0;
 494   for(itrack=0; itrack<2; itrack++) {
 495     if(itrack==0) {
 496       track = track1;
 497     } else {
 498       track = track2;
 499     }
 500     curvArray[itrack] = track->GetC(esd->GetMagneticField());
 501     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 502
 503     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 504       jpt=0;
 505       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
 506     }
 507     array = track->GetTrackPointArray();
 508     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 509       array->GetPoint(point,ipt);
 510       volId = point.GetVolumeID();
 511       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 512       if(layerId>6 || !layerOK[layerId-1][itrack]) continue;
 513       arrayForTree->AddPoint(jpt,&point);
 514       jpt++;
 515       switch(layerId) {
 516       case 1:
 517         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 518         break;
 519       case 2:
 520         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 521         break;
 522       case 3:
 523         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 524         break;
 525       case 4:
 526         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 527         break;
 528       case 5:
 529         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 530         break;
 531       case 6:
 532         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 533         break;
 534       }
 535       // Post the data for slot 0
 536       if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
 537       if(!point.IsExtra() ||
 538          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 539       nclsTrk[itrack]--;
 540       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 541       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 542       globExtra[0]=point.GetX();
 543       globExtra[1]=point.GetY();
 544       globExtra[2]=point.GetZ();
 545       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 546       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
 547       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 548       AliTrackPoint pointT;
 549       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 550       for(Int_t lll=0;lll<ipt;lll++) {
 551         array->GetPoint(pointT,lll);
 552         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 553         if(layerIdT!=layerId) continue;
 554         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 555         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 556         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 557         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 558         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 559         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 560         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 561         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 562         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 563         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 564         fntExtra->Fill((Float_t)nclsTrk[itrack],(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0mu[0],d0z0mu[1]);
 565       }
 566     }
 567
 568     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 569       curv = curvArray[itrack];
 570       curverr = curverrArray[itrack];
 571       fspTree->Fill();
 572     }
 573   }
 574
 575   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 576     curv = 0.5*(curvArray[0]+curvArray[1]);
 577     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
 578     fspTree->Fill();
 579   }
 580   PostData(0,fspTree);
 581
 582   if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
 583   // fill ntuple with track-to-track matching
 584   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
 585   Float_t d0[2],z0[2];
 586   Float_t sigmad0[2],sigmaz0[2];
 587   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
 588   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
 589   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
 590   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
 591   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
 592   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
 593
 594   dphi = phimu - (phiout+TMath::Pi());
 595   dtheta = thetamu - (TMath::Pi()-thetaout);
 596   if(rotymu>0) {
 597     droty = rotymu - (rotyout+TMath::Pi());
 598   } else {
 599     droty = rotymu - (rotyout-TMath::Pi());
 600   }
 601
 602   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
 603
 604   track1->Propagate(alpha,0.,esd->GetMagneticField());
 605   track2->Propagate(alpha,0.,esd->GetMagneticField());
 606   d0[0] = track1->GetY();
 607   z0[0] = track1->GetZ();
 608   d0[1] = track2->GetY();
 609   z0[1] = track2->GetZ();
 610   Float_t dxy = -(d0[0]-d0[1]);
 611   Float_t dz  = z0[0]-z0[1];
 612   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
 613   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
 614   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
 615   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
 616   /*
 617   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
 618   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
 619   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
 620   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
 621   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
 622   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
 623   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
 624   Float_t dxaty0 = x1aty0-x2aty0;
 625   */
 626   fntCosmicMatching->Fill((Float_t)nclsTrk[0],(Float_t)nclsTrk[1],track1->Pt(),track2->Pt(),sigmad0[0],sigmad0[1],sigmaz0[0],sigmaz0[1],dxy,dz,phimu,thetamu,TMath::Abs(d0z0mu[0]),d0z0mu[1]);
 627
 628   return;
 629 }
 630
 631 //________________________________________________________________________
 632 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
 633 {
 634   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 635   // of top and bottom track, merge the two tracks, if requested)
 636   //
 637
 638   // Set branch addresses for space points tree
 639   AliTrackPointArray *arrayForTree=0;
 640   Float_t curv,curverr;
 641   fspTree->SetBranchAddress("SP",&arrayForTree);
 642   fspTree->SetBranchAddress("curv",&curv);
 643   fspTree->SetBranchAddress("curverr",&curverr);
 644
 645   Int_t ntracks = esd->GetNumberOfTracks();
 646
 647   if(ntracks==0) return;
 648
 649   if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return;
 650
 651   Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos);
 652
 653   Int_t ncls=0;
 654   Double_t pt=-10000.;
 655   Double_t d0z0[2],covd0z0[3];
 656   const AliTrackPointArray *array = 0;
 657
 658   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 659     AliESDtrack * track = esd->GetTrack(itrack);
 660     if (!track) continue;
 661
 662     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 663
 664     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 665     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 666
 667     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 668        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 669
 670     pt = track->Pt();
 671     ncls = track->GetNcls(0);
 672     Double_t maxd=10000.;
 673     track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
 674
 675     // read ITS cluster map
 676     Int_t map[6];
 677     for(Int_t ilay=0;ilay<6;ilay++) {
 678       map[ilay]=0;
 679       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
 680     }
 681     AliDebug(2,Form("ITS map : %d %d %d %d %d %d pt %f\n",map[0],map[1],map[2],map[3],map[4],map[5],track->Pt()));
 682     Int_t idx[12];
 683     track->GetITSclusters(idx);
 684     AliDebug(2,Form("cls idx %d %d %d %d %d %d %d %d %d %d %d %d\n",idx[0],idx[1],idx[2],idx[3],idx[4],idx[5],idx[6],idx[7],idx[8],idx[9],idx[10],idx[11]));
 685
 686
 687     AliTrackPoint point;
 688     Int_t ipt,volId,modId,layerId,lay,lad,det;
 689     Int_t jpt=0;
 690     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
 691
 692     array = track->GetTrackPointArray();
 693     if(!array) continue;
 694     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 695       array->GetPoint(point,ipt);
 696       volId = point.GetVolumeID();
 697       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 698       if(layerId<1 || layerId>6) continue;
 699       if(point.IsExtra() &&
 700          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 701       layerOK[layerId-1]=kTRUE;
 702       jpt++;
 703     }
 704
 705     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 706
 707     fHistNpoints->Fill(jpt);
 708     fHistPt->Fill(pt);
 709     PostData(1,fListOfHistos);
 710
 711     Float_t dzOverlap[2];
 712     Double_t globExtra[3],locExtra[3];
 713     arrayForTree = new AliTrackPointArray(jpt);
 714     jpt=0;
 715     array = track->GetTrackPointArray();
 716     if(!array) continue;
 717     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 718       array->GetPoint(point,ipt);
 719       volId = point.GetVolumeID();
 720       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 721       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
 722       if(!point.IsExtra() ||
 723          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 724       ncls--;
 725       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 726       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 727       globExtra[0]=point.GetX();
 728       globExtra[1]=point.GetY();
 729       globExtra[2]=point.GetZ();
 730       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 731       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 732       AliTrackPoint pointT;
 733       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 734       for(Int_t lll=0;lll<ipt;lll++) {
 735         array->GetPoint(pointT,lll);
 736         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 737         if(layerIdT!=layerId) continue;
 738         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 739         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 740         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 741         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 742         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 743         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 744         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 745         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 746         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 747         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 748         fntExtra->Fill((Float_t)ncls,(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0[0],d0z0[1]);
 749       }
 750       arrayForTree->AddPoint(jpt,&point);
 751       switch(layerId) {
 752       case 1:
 753         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 754         break;
 755       case 2:
 756         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 757         break;
 758       case 3:
 759         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 760         break;
 761       case 4:
 762         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 763         break;
 764       case 5:
 765         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 766         break;
 767       case 6:
 768         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 769         break;
 770       }
 771       jpt++;
 772     }
 773
 774     curv = track->GetC(esd->GetMagneticField());
 775     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 776
 777     fspTree->Fill();
 778
 779   } // end of tracks loop
 780
 781   PostData(0,fspTree);
 782
 783   return;
 784 }
 785
 786 //________________________________________________________________________
 787 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
 788 {
 789   // Terminate analysis
 790   //
 791   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
 792
 793   fspTree = dynamic_cast<TTree*> (GetOutputData(0));
 794   if (!fspTree) {
 795     printf("ERROR: fspTree not available\n");
 796     return;
 797   }
 798
 799   fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
 800   if (!fListOfHistos) {
 801     printf("ERROR: fListOfHistos not available\n");
 802     return;
 803   }
 804
 805   fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
 806   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
 807   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
 808   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
 809   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
 810   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
 811   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
 812   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
 813   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
 814   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
 815   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
 816
 817
 818
 819   return;
 820 }
 821