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