PWGPP/ITS/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 <TFile.h>
  26 #include <TSystem.h>
  27 #include <TChain.h>
  28 #include <TNtuple.h>
  29 #include <TList.h>
  30 #include <TH1F.h>
  31 #include <TH2F.h>
  32 #include <TMap.h>
  33 #include <TVector3.h>
  34 #include <TGeoManager.h>
  35 #include <TRandom.h>
  36
  37 #include "AliLog.h"
  38 #include "AliCDBEntry.h"
  39 #include "AliGeomManager.h"
  40 #include "AliITSReconstructor.h"
  41 #include "AliITSAlignMille2Module.h"
  42 #include "AliITSgeomTGeo.h"
  43 #include "AliTrackPointArray.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDtrack.h"
  47 #include "AliESDEvent.h"
  48 #include "AliESDfriend.h"
  49 #include "AliAnalysisTaskSE.h"
  50 #include "AliAnalysisManager.h"
  51 #include "AliAlignmentDataFilterITS.h"
  52
  53
  54 ClassImp(AliAlignmentDataFilterITS)
  55
  56
  57 //________________________________________________________________________
  58 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS():
  59 AliAnalysisTaskSE(),
  60 fOnlySPDFO(kFALSE),
  61 fDownsamplelowpt(kFALSE),
  62 fGeometryFileName("geometry.root"),
  63 fITSRecoParam(0),
  64 fESD(0),
  65 fListOfHistos(0),
  66 fspTree(0),
  67 fHistNevents(0),
  68 fHistNpoints(0),
  69 fHistPt(0),
  70 fHistLayer0(0),
  71 fHistLayer1(0),
  72 fHistLayer2(0),
  73 fHistLayer3(0),
  74 fHistLayer4(0),
  75 fHistLayer5(0),
  76 fntExtra(0),
  77 fntCosmicMatching(0)
  78 {
  79   // Constructor
  80 }
  81
  82 //________________________________________________________________________
  83 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
  84 AliAnalysisTaskSE(name),
  85 fOnlySPDFO(kFALSE),
  86 fDownsamplelowpt(kFALSE),
  87 fGeometryFileName("geometry.root"),
  88 fITSRecoParam(0),
  89 fESD(0),
  90 fListOfHistos(0),
  91 fspTree(0),
  92 fHistNevents(0),
  93 fHistNpoints(0),
  94 fHistPt(0),
  95 fHistLayer0(0),
  96 fHistLayer1(0),
  97 fHistLayer2(0),
  98 fHistLayer3(0),
  99 fHistLayer4(0),
 100 fHistLayer5(0),
 101 fntExtra(0),
 102 fntCosmicMatching(0)
 103 {
 104   // Constructor
 105
 106   // Define output slots here
 107
 108   // Output slot #1 writes into a TTree
 109   DefineOutput(1,TTree::Class());  //My private output
 110   // Output slot #2 writes into a TList
 111   DefineOutput(2,TList::Class());  //My private output
 112 }
 113
 114 //________________________________________________________________________
 115 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
 116 {
 117   // Destructor
 118   if (fListOfHistos) {
 119     delete fListOfHistos;
 120     fListOfHistos = 0;
 121   }
 122   if (fspTree) {
 123     delete fspTree;
 124     fspTree = 0;
 125   }
 126   if (fHistNevents) {
 127     delete fHistNevents;
 128     fHistNevents = 0;
 129   }
 130   if (fHistNpoints) {
 131     delete fHistNpoints;
 132     fHistNpoints = 0;
 133   }
 134   if (fHistPt) {
 135     delete fHistPt;
 136     fHistPt = 0;
 137   }
 138   if (fHistLayer0) {
 139     delete fHistLayer0;
 140     fHistLayer0 = 0;
 141   }
 142   if (fHistLayer1) {
 143     delete fHistLayer1;
 144     fHistLayer1 = 0;
 145   }
 146   if (fHistLayer2) {
 147     delete fHistLayer2;
 148     fHistLayer2 = 0;
 149   }
 150   if (fHistLayer3) {
 151     delete fHistLayer3;
 152     fHistLayer3 = 0;
 153   }
 154   if (fHistLayer4) {
 155     delete fHistLayer4;
 156     fHistLayer4 = 0;
 157   }
 158   if (fHistLayer5) {
 159     delete fHistLayer5;
 160     fHistLayer5 = 0;
 161   }
 162   if (fntExtra) {
 163     delete fntExtra;
 164     fntExtra = 0;
 165   }
 166   if (fntCosmicMatching) {
 167     delete fntCosmicMatching;
 168     fntCosmicMatching = 0;
 169   }
 170 }
 171
 172
 173 //________________________________________________________________________
 174 void AliAlignmentDataFilterITS::UserCreateOutputObjects()
 175 {
 176   // Create the output container
 177   //
 178
 179   // load the geometry
 180   if(!gGeoManager) {
 181     printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data());
 182     AliGeomManager::LoadGeometry(fGeometryFileName.Data());
 183     if(!gGeoManager) {
 184       printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n");
 185       return;
 186     }
 187   }
 188
 189   // Several histograms are more conveniently managed in a TList
 190   fListOfHistos = new TList();
 191   fListOfHistos->SetOwner();
 192
 193   fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5);
 194   fHistNevents->Sumw2();
 195   fHistNevents->SetMinimum(0);
 196   fListOfHistos->Add(fHistNevents);
 197
 198   fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
 199   fHistNpoints->Sumw2();
 200   fHistNpoints->SetMinimum(0);
 201   fListOfHistos->Add(fHistNpoints);
 202
 203   fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
 204   fHistPt->Sumw2();
 205   fHistPt->SetMinimum(0);
 206   fListOfHistos->Add(fHistPt);
 207
 208
 209   Float_t zmax=14.;
 210   Int_t nbinsphi=20,nbinsz=4;
 211   fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 212   fListOfHistos->Add(fHistLayer0);
 213   zmax=14.;
 214   nbinsphi=40;nbinsz=4;
 215   fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 216   fListOfHistos->Add(fHistLayer1);
 217   zmax=22.;
 218   nbinsphi=14;nbinsz=6;
 219   fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 220   fListOfHistos->Add(fHistLayer2);
 221   zmax=29.5;
 222   nbinsphi=22;nbinsz=8;
 223   fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 224   fListOfHistos->Add(fHistLayer3);
 225   zmax=45.;
 226   nbinsphi=34;nbinsz=23;
 227   fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 228   fListOfHistos->Add(fHistLayer4);
 229   zmax=51.;
 230   nbinsphi=38;nbinsz=26;
 231   fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global   #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
 232   fListOfHistos->Add(fHistLayer5);
 233
 234
 235   fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu:pt");
 236   fListOfHistos->Add(fntExtra);
 237
 238   fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
 239   fListOfHistos->Add(fntCosmicMatching);
 240
 241   fspTree = new TTree("spTree","Tree with ITS track points");
 242   AliTrackPointArray *array = 0;
 243   AliESDVertex *vertex = 0;
 244   Float_t curv=0,curverr=0,runNumber=0;
 245   TObjString *itsaligndata = 0;
 246   TObjString *itscalibrespsdd = 0;
 247   fspTree->Branch("SP","AliTrackPointArray",&array);
 248   fspTree->Branch("vertex","AliESDVertex",&vertex);
 249   fspTree->Branch("curv",&curv);
 250   fspTree->Branch("curverr",&curverr);
 251   fspTree->Branch("run",&runNumber);
 252   fspTree->Branch("ITSAlignData",&itsaligndata);
 253   fspTree->Branch("ITSCalibRespSDD",&itscalibrespsdd);
 254
 255   return;
 256 }
 257
 258 //________________________________________________________________________
 259 void AliAlignmentDataFilterITS::UserExec(Option_t */*option*/)
 260 {
 261   // Execute analysis for current event:
 262   // write ITS AliTrackPoints for selected tracks to fspTree
 263
 264   // check the geometry
 265   if(!gGeoManager) {
 266     printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
 267     return;
 268   }
 269
 270   // check if we have AliITSRecoParam
 271   if(!GetRecoParam()) {
 272     if(!fITSRecoParam) {
 273       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
 274       return;
 275     }
 276   }
 277
 278   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 279   if(!fESD) {
 280     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
 281     return;
 282   }
 283
 284   //AliESDfriend *esdfriend = (AliESDfriend*)(fESD->FindListObject("AliESDfriend"));
 285
 286   //if(!esdfriend) printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
 287   //fESD->SetESDfriend(esdfriend);
 288
 289   // Post the data for slot 0
 290   fHistNevents->Fill(0);
 291
 292
 293   // write field value to spTree UserInfo
 294   if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) {
 295     Double_t bz=fESD->GetMagneticField();
 296     TString bzString; bzString+=bz;
 297     TObjString *bzObjString = new TObjString(bzString);
 298     TList *bzList = new TList();
 299     bzList->SetOwner(1);
 300     bzList->SetName("BzkGauss");
 301     bzList->Add(bzObjString);
 302     fspTree->GetUserInfo()->Add(bzList);
 303   }
 304
 305   // write OCDB info to spTree UserInfo
 306   if(!((fspTree->GetUserInfo())->FindObject("cdbList"))) {
 307     TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 308     if(!tree) {
 309       printf("ERROR: Could not read chain from input slot 0\n");
 310     } else {
 311       // Get the OCDB path and the list of OCDB objects used for reco
 312       TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap");
 313       TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList");
 314
 315       //cdbList->Print();
 316       // write the list to the user info of the output tree
 317       if(!fspTree) {
 318         printf("ERROR: fspTree does not exist\n");
 319       } else {
 320         TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
 321         cdbMapCopy->SetOwner(1);
 322         cdbMapCopy->SetName("cdbMap");
 323         TIter iter1(cdbMap->GetTable());
 324
 325         TPair* pair = 0;
 326         while((pair = dynamic_cast<TPair*> (iter1.Next()))){
 327           TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
 328           TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
 329           if(keyStr && valStr) cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
 330         }
 331
 332         TList *cdbListCopy = new TList();
 333         cdbListCopy->SetOwner(1);
 334         cdbListCopy->SetName("cdbList");
 335
 336         TIter iter2(cdbList);
 337
 338         TObjString* cdbEntry=0;
 339         while((cdbEntry =(TObjString*)(iter2.Next()))) {
 340           cdbListCopy->Add(new TObjString(*cdbEntry));
 341         }
 342         cdbListCopy->Print();
 343
 344
 345         fspTree->GetUserInfo()->Add(cdbMapCopy);
 346         fspTree->GetUserInfo()->Add(cdbListCopy);
 347       }
 348     }
 349   }
 350
 351
 352
 353   // Process event as Cosmic or Collision
 354   if(fESD->GetEventSpecie()<=1) {
 355     printf("AliAlignmentDataFilterITS::Exec(): event specie not set !\n");
 356     if(GetRecoParam()->GetAlignFilterCosmics()) {
 357       FilterCosmic(fESD);
 358     } else {
 359       FilterCollision(fESD);
 360     }
 361   } else if(fESD->GetEventSpecie()==8) {
 362     FilterCosmic(fESD);
 363   } else {
 364     FilterCollision(fESD);
 365   }
 366
 367   PostData(2,fListOfHistos);
 368
 369   return;
 370 }
 371
 372 //________________________________________________________________________
 373 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
 374 {
 375   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 376   // of top and bottom track, merge the two tracks, if requested)
 377   //
 378
 379   // Set branch addresses for space points tree
 380   AliTrackPointArray *arrayForTree=0;
 381   AliESDVertex *vertexForTree=0;
 382   Float_t curv,curverr,runNumber;
 383   TObjString *itsaligndata=0;
 384   TObjString *itscalibrespsdd = 0;
 385   fspTree->SetBranchAddress("SP",&arrayForTree);
 386   fspTree->SetBranchAddress("vertex",&vertexForTree);
 387   fspTree->SetBranchAddress("curv",&curv);
 388   fspTree->SetBranchAddress("curverr",&curverr);
 389   fspTree->SetBranchAddress("run",&runNumber);
 390   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 391   fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
 392
 393
 394   runNumber = (Float_t)esd->GetRunNumber();
 395   Int_t uid=10000+esd->GetEventNumberInFile();
 396
 397   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 398   if(!esdTree) return;
 399   // Get the list of OCDB objects used for reco
 400   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 401   TIter iter2(cdbList);
 402   TObjString* cdbEntry=0;
 403   TString cdbEntryString;
 404   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 405     cdbEntryString = cdbEntry->GetString();
 406     if(cdbEntryString.Contains("ITS/Align/Data")) {
 407       itsaligndata = new TObjString(*cdbEntry);
 408       itsaligndata->SetString(itsaligndata->GetString());
 409     }
 410     if(cdbEntryString.Contains("ITS/Calib/RespSDD")) {
 411       itscalibrespsdd = new TObjString(*cdbEntry);
 412       itscalibrespsdd->SetString(itscalibrespsdd->GetString());
 413     }
 414   }
 415
 416
 417   TString triggeredClass = esd->GetFiredTriggerClasses();
 418   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
 419
 420
 421   Int_t ntracks = esd->GetNumberOfTracks();
 422   if(ntracks<2) return;
 423
 424   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
 425
 426   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
 427
 428   Int_t *goodtracksArray = new Int_t[ntracks];
 429   Float_t *phiArray = new Float_t[ntracks];
 430   Float_t *thetaArray = new Float_t[ntracks];
 431   Int_t *nclsArray = new Int_t[ntracks];
 432   Int_t ngt=0;
 433   Int_t itrack=0;
 434   for (itrack=0; itrack < ntracks; itrack++) {
 435     AliESDtrack *track = esd->GetTrack(itrack);
 436     if (!track) continue;
 437
 438
 439     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 440
 441     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 442     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 443
 444     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
 445     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
 446
 447     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 448        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 449
 450     goodtracksArray[ngt] = itrack;
 451     phiArray[ngt]        = phi;
 452     thetaArray[ngt]      = theta;
 453     nclsArray[ngt]       = track->GetNcls(0);
 454     ngt++;
 455   }
 456
 457   if(ngt<2) {
 458     delete [] goodtracksArray; goodtracksArray=0;
 459     delete [] phiArray; phiArray=0;
 460     delete [] thetaArray; thetaArray=0;
 461     delete [] nclsArray; nclsArray=0;
 462     return;
 463   }
 464
 465   // check matching of the two tracks from the muon
 466   Float_t min = 10000000.;
 467   Int_t maxCls = 0;
 468   Int_t good1 = -1, good2 = -1;
 469   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
 470     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
 471       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
 472       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 473       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
 474       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
 475       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
 476         maxCls = nclsArray[itr1]+nclsArray[itr2];
 477         min = deltatheta+deltaphi;
 478         good1 = goodtracksArray[itr1];
 479         good2 = goodtracksArray[itr2];
 480       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
 481         if(deltatheta+deltaphi < min) {
 482           min = deltatheta+deltaphi;
 483           good1 = goodtracksArray[itr1];
 484           good2 = goodtracksArray[itr2];
 485         }
 486       }
 487     }
 488   }
 489
 490   delete [] goodtracksArray; goodtracksArray=0;
 491   delete [] phiArray; phiArray=0;
 492   delete [] thetaArray; thetaArray=0;
 493   delete [] nclsArray; nclsArray=0;
 494
 495   if(good1<0) return;
 496   AliDebug(2,"ok track matching");
 497
 498   // track1 will be the inward track (top)
 499   // track2 the outward (bottom)
 500   AliESDtrack *track1=0;
 501   AliESDtrack *track2=0;
 502   AliESDtrack *track = esd->GetTrack(good1);
 503   if(track->Py()>0) {
 504     track1 = esd->GetTrack(good1);
 505     track2 = esd->GetTrack(good2);
 506   } else {
 507     track1 = esd->GetTrack(good2);
 508     track2 = esd->GetTrack(good1);
 509   }
 510
 511   AliTrackPoint point;
 512   const AliTrackPointArray *array=0;
 513   Int_t ipt,volId,modId,layerId,lay,lad,det;
 514   Int_t jpt=0;
 515   Bool_t layerOK[6][2];
 516   Int_t nclsTrk[2]={0,0};
 517
 518   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
 519
 520   for(itrack=0; itrack<2; itrack++) {
 521     if(itrack==0) {
 522       track = track1;
 523     } else {
 524       track = track2;
 525     }
 526     array = track->GetTrackPointArray();
 527     if(!array) {
 528       AliWarning("No tracks points avaialble");
 529       continue;
 530     }
 531     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 532       array->GetPoint(point,ipt);
 533       volId = point.GetVolumeID();
 534       if(volId<=0) continue;
 535       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 536       AliDebug(2,Form("%d %d %d  %f\n",ipt,layerId-1,volId,TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY())));
 537       if(point.IsExtra() &&
 538          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 539       if(layerId<1 || layerId>6) continue;
 540       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
 541       // check minAngleWrtITSModulePlanes
 542       Double_t p[3]; track->GetDirection(p);
 543       TVector3 pvec(p[0],p[1],p[2]);
 544       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
 545       TVector3 normvec(rot[1],rot[4],rot[7]);
 546       Double_t angle = pvec.Angle(normvec);
 547       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
 548       angle = 0.5*TMath::Pi()-angle;
 549       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
 550       layerOK[layerId-1][itrack]=kTRUE;
 551       jpt++;
 552       nclsTrk[itrack]++;
 553     }
 554   }
 555   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
 556
 557   // read ITS cluster maps
 558   Int_t map1[6],map2[6];
 559   for(Int_t ilay=0;ilay<6;ilay++) {
 560     map1[ilay]=0; map2[ilay]=0;
 561     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
 562     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
 563   }
 564   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()));
 565   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()));
 566   Int_t idx1[12],idx2[12];
 567   track1->GetITSclusters(idx1);
 568   track2->GetITSclusters(idx2);
 569   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]));
 570   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]));
 571
 572
 573   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
 574   AliDebug(2,Form(" Total points %d, accepted\n",jpt));
 575   fHistNpoints->Fill(jpt);
 576   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
 577
 578   Float_t d0z0mu[2];
 579   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
 580   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
 581
 582   vertexForTree = new AliESDVertex(*(esd->GetPrimaryVertexSPD()));
 583   vertexForTree->SetID(0);
 584
 585   Float_t dzOverlap[2];
 586   Float_t curvArray[2],curverrArray[2];
 587   Double_t globExtra[3],locExtra[3];
 588   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 589     arrayForTree = new AliTrackPointArray(jpt);
 590     arrayForTree->SetUniqueID(uid);
 591   }
 592   jpt=0;
 593   for(itrack=0; itrack<2; itrack++) {
 594     if(itrack==0) {
 595       track = track1;
 596     } else {
 597       track = track2;
 598     }
 599     curvArray[itrack] = track->GetC(esd->GetMagneticField());
 600     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 601
 602     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 603       jpt=0;
 604       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
 605       arrayForTree->SetUniqueID(uid);
 606     }
 607     array = track->GetTrackPointArray();
 608     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 609       array->GetPoint(point,ipt);
 610       volId = point.GetVolumeID();
 611       if(volId<=0) continue;
 612       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 613       if(layerId<1 || layerId>6 || !layerOK[layerId-1][itrack]) continue;
 614       arrayForTree->AddPoint(jpt,&point);
 615       jpt++;
 616       switch(layerId) {
 617       case 1:
 618         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 619         break;
 620       case 2:
 621         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 622         break;
 623       case 3:
 624         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 625         break;
 626       case 4:
 627         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 628         break;
 629       case 5:
 630         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 631         break;
 632       case 6:
 633         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 634         break;
 635       }
 636       // Post the data for slot 2
 637       if(jpt==1) PostData(2,fListOfHistos); // only if this is the first points
 638       if(!point.IsExtra() ||
 639          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 640       nclsTrk[itrack]--;
 641       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 642       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 643       globExtra[0]=point.GetX();
 644       globExtra[1]=point.GetY();
 645       globExtra[2]=point.GetZ();
 646       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 647       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
 648       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 649       AliTrackPoint pointT;
 650       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 651       for(Int_t lll=0;lll<ipt;lll++) {
 652         array->GetPoint(pointT,lll);
 653         if(pointT.GetVolumeID()<=0) continue;
 654         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 655         if(layerIdT!=layerId) continue;
 656         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 657         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 658         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 659         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 660         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 661         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 662         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 663         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 664         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 665         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 666         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],track->Pt());
 667       }
 668     }
 669
 670     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
 671       curv = curvArray[itrack];
 672       curverr = curverrArray[itrack];
 673       fspTree->Fill();
 674     }
 675   }
 676
 677   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
 678     curv = 0.5*(curvArray[0]-curvArray[1]);  // the "-" is because the two tracks have opposite curvature!
 679     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
 680     fspTree->Fill();
 681   }
 682   PostData(1,fspTree);
 683
 684   if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
 685   // fill ntuple with track-to-track matching
 686   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
 687   Float_t d0[2],z0[2];
 688   Float_t sigmad0[2],sigmaz0[2];
 689   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
 690   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
 691   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
 692   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
 693   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
 694   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
 695
 696   dphi = phimu - (phiout+TMath::Pi());
 697   dtheta = thetamu - (TMath::Pi()-thetaout);
 698   if(rotymu>0) {
 699     droty = rotymu - (rotyout+TMath::Pi());
 700   } else {
 701     droty = rotymu - (rotyout-TMath::Pi());
 702   }
 703
 704   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
 705
 706   track1->Propagate(alpha,0.,esd->GetMagneticField());
 707   track2->Propagate(alpha,0.,esd->GetMagneticField());
 708   d0[0] = track1->GetY();
 709   z0[0] = track1->GetZ();
 710   d0[1] = track2->GetY();
 711   z0[1] = track2->GetZ();
 712   Float_t dxy = -(d0[0]-d0[1]);
 713   Float_t dz  = z0[0]-z0[1];
 714   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
 715   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
 716   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
 717   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
 718   /*
 719   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
 720   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
 721   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
 722   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
 723   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
 724   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
 725   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
 726   Float_t dxaty0 = x1aty0-x2aty0;
 727   */
 728   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]);
 729
 730   return;
 731 }
 732
 733 //________________________________________________________________________
 734 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
 735 {
 736   // Extract ITS AliTrackPoints for Cosmics (check angular matching
 737   // of top and bottom track, merge the two tracks, if requested)
 738   //
 739
 740   // Set branch addresses for space points tree
 741   AliTrackPointArray *arrayForTree=0;
 742   AliESDVertex *vertexForTree=0;
 743   Float_t curv,curverr,runNumber;
 744   TObjString *itsaligndata=0;
 745   TObjString *itscalibrespsdd = 0;
 746   fspTree->SetBranchAddress("SP",&arrayForTree);
 747   fspTree->SetBranchAddress("vertex",&vertexForTree);
 748   fspTree->SetBranchAddress("curv",&curv);
 749   fspTree->SetBranchAddress("curverr",&curverr);
 750   fspTree->SetBranchAddress("run",&runNumber);
 751   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
 752   fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
 753
 754
 755   runNumber = (Float_t)esd->GetRunNumber();
 756   Int_t uid=20000+esd->GetEventNumberInFile();
 757
 758   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
 759   if(!esdTree) return;
 760   // Get the list of OCDB objects used for reco
 761   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
 762   TIter iter2(cdbList);
 763   TObjString* cdbEntry=0;
 764   TString cdbEntryString;
 765   while((cdbEntry =(TObjString*)(iter2.Next()))) {
 766     cdbEntryString = cdbEntry->GetString();
 767     if(cdbEntryString.Contains("ITS/Align/Data")) {
 768       itsaligndata = new TObjString(*cdbEntry);
 769       itsaligndata->SetString(itsaligndata->GetString());
 770     }
 771     if(cdbEntryString.Contains("ITS/Calib/RespSDD")) {
 772       itscalibrespsdd = new TObjString(*cdbEntry);
 773       itscalibrespsdd->SetString(itscalibrespsdd->GetString());
 774     }
 775   }
 776
 777   Int_t ntracks = esd->GetNumberOfTracks();
 778
 779   if(ntracks==0) return;
 780
 781   const AliESDVertex *vertexTracks = esd->GetPrimaryVertexTracks();
 782   if(!vertexTracks) return;
 783   if(vertexTracks->GetNContributors()<=0) return;
 784
 785   Double_t vtxpos[3]; vertexTracks->GetXYZ(vtxpos);
 786
 787   Int_t ncls=0;
 788   Double_t pt=-10000.;
 789   Double_t d0z0[2],covd0z0[3];
 790   const AliTrackPointArray *array = 0;
 791
 792   for (Int_t itrack=0; itrack < ntracks; itrack++) {
 793     AliESDtrack * track = esd->GetTrack(itrack);
 794     if (!track) continue;
 795
 796     if(fDownsamplelowpt && TMath::Abs(esd->GetMagneticField())>0.01 &&
 797        track->Pt()<gRandom->Rndm()) continue;
 798
 799     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 800
 801     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
 802     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
 803
 804     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
 805        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
 806
 807     pt = track->Pt();
 808     ncls = track->GetNcls(0);
 809     Double_t maxd=10000.;
 810     track->PropagateToDCA(vertexTracks,esd->GetMagneticField(),maxd,d0z0,covd0z0);
 811
 812     // read ITS cluster map
 813     Int_t map[6];
 814     for(Int_t ilay=0;ilay<6;ilay++) {
 815       map[ilay]=0;
 816       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
 817     }
 818     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()));
 819     Int_t idx[12];
 820     track->GetITSclusters(idx);
 821     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]));
 822
 823
 824     AliTrackPoint point;
 825     Int_t ipt,volId,modId,layerId,lay,lad,det;
 826     Int_t jpt=0;
 827     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
 828
 829     array = track->GetTrackPointArray();
 830     if(!array) {printf("no track points\n"); continue;}
 831     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 832       array->GetPoint(point,ipt);
 833       volId = point.GetVolumeID();
 834       if(volId<=0) continue;
 835       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 836       if(layerId<1 || layerId>6) continue;
 837       if(point.IsExtra() &&
 838          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
 839       layerOK[layerId-1]=kTRUE;
 840       jpt++;
 841     }
 842
 843     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
 844
 845     fHistNpoints->Fill(jpt);
 846     fHistPt->Fill(pt);
 847     PostData(2,fListOfHistos);
 848
 849     Float_t dzOverlap[2];
 850     Double_t globExtra[3],locExtra[3];
 851     arrayForTree = new AliTrackPointArray(jpt);
 852     arrayForTree->SetUniqueID(uid);
 853     jpt=0;
 854     array = track->GetTrackPointArray();
 855     if(!array) continue;
 856     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
 857       array->GetPoint(point,ipt);
 858       volId = point.GetVolumeID();
 859       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
 860       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
 861       arrayForTree->AddPoint(jpt,&point);
 862       switch(layerId) {
 863       case 1:
 864         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 865         break;
 866       case 2:
 867         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 868         break;
 869       case 3:
 870         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 871         break;
 872       case 4:
 873         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 874         break;
 875       case 5:
 876         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 877         break;
 878       case 6:
 879         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
 880         break;
 881       }
 882       jpt++;
 883       if(!point.IsExtra() ||
 884          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
 885       ncls--;
 886       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
 887       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
 888       globExtra[0]=point.GetX();
 889       globExtra[1]=point.GetY();
 890       globExtra[2]=point.GetZ();
 891       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
 892       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
 893       AliTrackPoint pointT;
 894       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
 895       for(Int_t lll=0;lll<ipt;lll++) {
 896         array->GetPoint(pointT,lll);
 897         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
 898         if(layerIdT!=layerId) continue;
 899         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
 900         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
 901         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
 902         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
 903         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
 904         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
 905         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
 906         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
 907         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
 908         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
 909         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],track->Pt());
 910       }
 911     }
 912
 913     curv = track->GetC(esd->GetMagneticField());
 914     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
 915
 916     vertexForTree = new AliESDVertex(*vertexTracks);
 917     if(vertexTracks->UsesTrack(track->GetID())) {
 918       vertexForTree->SetID(1);
 919     } else {
 920       vertexForTree->SetID(0);
 921     }
 922
 923     fspTree->Fill();
 924
 925   } // end of tracks loop
 926
 927   PostData(1,fspTree);
 928
 929   return;
 930 }
 931
 932 //________________________________________________________________________
 933 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
 934 {
 935   // Terminate analysis
 936   //
 937   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
 938
 939   fspTree = dynamic_cast<TTree*> (GetOutputData(1));
 940   if (!fspTree) {
 941     printf("ERROR: fspTree not available\n");
 942     return;
 943   }
 944
 945   fListOfHistos = dynamic_cast<TList*> (GetOutputData(2));
 946   if (!fListOfHistos) {
 947     printf("ERROR: fListOfHistos not available\n");
 948     return;
 949   }
 950
 951   fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
 952   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
 953   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
 954   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
 955   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
 956   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
 957   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
 958   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
 959   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
 960   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
 961   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
 962
 963
 964
 965   return;
 966 }
 967 //-------------------------------------------------------------------------------
 968 const AliITSRecoParam *AliAlignmentDataFilterITS::GetRecoParam() const
 969 {
 970   //
 971   // Return the ITSRecoParam object
 972   //
 973   if(AliITSReconstructor::GetRecoParam()) {
 974     return AliITSReconstructor::GetRecoParam();
 975   } else if(fITSRecoParam) {
 976     return fITSRecoParam;
 977   } else return NULL;
 978 }
 979 //--------------------------------------------------------------------------------
 980 Int_t AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom(const Char_t *fin,
 981                                                              const Char_t *fout,
 982                                                              const Char_t *fmis,
 983                                                              const Char_t *fgeo,
 984                                                              Bool_t prn)
 985 {
 986   //
 987   // Convert AliTrackPoints in fin, reconstructed with fmis, back
 988   // to ideal geometry
 989   //
 990   // M. Lunardon
 991   //
 992
 993
 994   TGeoHMatrix deltahm;
 995
 996   // Load geometry
 997   if (gSystem->AccessPathName(fgeo)) {
 998     printf("couldn't find geometry file %s - skipping...\n",fmis);
 999     return -1;
1000   }
1001
1002   TFile *geofile=TFile::Open(fgeo);
1003   TGeoManager *fgGeometry=NULL;
1004
1005   fgGeometry=(TGeoManager*)geofile->Get("ALICE");
1006
1007   if (!fgGeometry)
1008     fgGeometry=(TGeoManager*)geofile->Get("Geometry");
1009
1010   if (!fgGeometry) {
1011     AliCDBEntry *entry = (AliCDBEntry*)geofile->Get("AliCDBEntry");
1012     if (entry)
1013       fgGeometry = (TGeoManager*)entry->GetObject();
1014   }
1015
1016   if (!fgGeometry) return -1;
1017   AliGeomManager::SetGeometry(fgGeometry);
1018   if(!AliGeomManager::GetGeometry()) return -1;
1019
1020
1021   // open alignment file
1022   if (gSystem->AccessPathName(fmis)) {
1023     printf("couldn't open alignment file %s - skipping...\n",fmis);
1024     return -2;
1025   }
1026   TFile *pref = TFile::Open(fmis);
1027   if (!pref->IsOpen()) return -2;
1028
1029
1030   /// apply alignment to ideal geometry
1031   TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
1032   if (!prea) {
1033     if (pref->Get("AliCDBEntry"))
1034       prea = (TClonesArray*) ((AliCDBEntry*)pref->Get("AliCDBEntry"))->GetObject();
1035   }
1036   if (!prea) return -3;
1037   Int_t nprea=prea->GetEntriesFast();
1038   printf("Array of input misalignments with %d entries\n",nprea);
1039   AliGeomManager::ApplyAlignObjsToGeom(*prea); // apply all levels of objs
1040
1041   AliTrackPointArray *tpain=NULL;
1042   TFile *tpainfile=NULL;
1043   TTree *treein=NULL;
1044   AliTrackPoint point;
1045   AliITSAlignMille2Module *m2[2200];
1046   for (Int_t i=0; i<2198; i++)
1047     m2[i]=new AliITSAlignMille2Module(AliITSAlignMille2Module::GetVolumeIDFromIndex(i));
1048
1049   // open input file
1050   if (gSystem->AccessPathName(fin)) {
1051     printf("couldn't open file %s - skipping...\n",fin);
1052     return -4;
1053   }
1054   tpainfile = TFile::Open(fin);
1055   if (!tpainfile->IsOpen()) return -4;
1056
1057   treein=(TTree*)tpainfile->Get("spTree");
1058   if (!treein) return -5;
1059   Float_t curv,curverr,runNumber;
1060   TObjString *itsaligndata=0;
1061   TObjString *itscalibrespsdd = 0;
1062   treein->SetBranchAddress("SP", &tpain);
1063   treein->SetBranchAddress("curv", &curv);
1064   treein->SetBranchAddress("curverr", &curverr);
1065   treein->SetBranchAddress("run",&runNumber);
1066   treein->SetBranchAddress("ITSAlignData",&itsaligndata);
1067   treein->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1068
1069   int ntrks=treein->GetEntries();
1070   printf("Reading %d tracks from %s\n",ntrks,fin);
1071
1072
1073   // open output file
1074   TFile *pointsFile = TFile::Open(fout,"RECREATE");
1075   if (!pointsFile || !pointsFile->IsOpen()) {
1076     printf("Can't open output file %s !",fout);
1077     return -6;
1078   }
1079   AliTrackPointArray *array = new AliTrackPointArray();
1080
1081   // new!
1082   TTree *treeout=(TTree*)treein->Clone("spTree");
1083   treeout->Reset();
1084   treeout->SetBranchAddress("SP", &array);
1085   treeout->SetBranchAddress("curv", &curv);
1086   treeout->SetBranchAddress("curverr", &curverr);
1087   treeout->SetBranchAddress("run",&runNumber);
1088   treeout->SetBranchAddress("ITSAlignData",&itsaligndata);
1089   treeout->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1090
1091   // tracks main loop
1092   for (Int_t it=0; it<ntrks; it++) {
1093     if (!(it%5000) ) printf("...processing track n. %d\n",it);
1094
1095     treein->GetEvent(it);
1096
1097     //////////////////////////////
1098
1099     AliTrackPointArray *atp=tpain;
1100     AliTrackPointArray *atps=NULL;
1101     Int_t npts=atp->GetNPoints();
1102
1103     AliTrackPoint p;
1104     // check points in specific places
1105
1106     // build a new track
1107     atps=new AliTrackPointArray(npts);
1108
1109     Int_t npto=0;
1110     for (int i=0; i<npts; i++) {
1111       atp->GetPoint(p,i);
1112
1113       UShort_t volid=atp->GetVolumeID()[i];
1114       Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(volid);
1115
1116       if(index<0 || index>=2200) continue;
1117       // dealign point
1118       // get MODIFIED matrix
1119       TGeoHMatrix *svMatrix = m2[index]->GetSensitiveVolumeMatrix(p.GetVolumeID());
1120       //TGeoHMatrix *svOrigMatrix = mm->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1121
1122       Double_t pg[3],pl[3];
1123       pg[0]=p.GetX();
1124       pg[1]=p.GetY();
1125       pg[2]=p.GetZ();
1126       if (prn) printf("Global coordinates of measured point : X=%f  Y=%f  Z=%f \n",pg[0],pg[1],pg[2]);
1127       svMatrix->MasterToLocal(pg,pl);
1128
1129       // check that things went OK: local y should be 0.
1130       if(TMath::Abs(pl[1])>1.e-6) {
1131         printf("AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom: ERROR, local y = %f (should be zero)\n",pl[1]);
1132         return -7;
1133       }
1134
1135       if (prn) printf("Local coordinates of measured point : X=%f  Y=%f  Z=%f \n",pl[0],pl[1],pl[2]);
1136
1137       // update covariance matrix
1138       TGeoHMatrix hcov;
1139       Double_t hcovel[9];
1140       hcovel[0]=(Double_t)(p.GetCov()[0]);
1141       hcovel[1]=(Double_t)(p.GetCov()[1]);
1142       hcovel[2]=(Double_t)(p.GetCov()[2]);
1143       hcovel[3]=(Double_t)(p.GetCov()[1]);
1144       hcovel[4]=(Double_t)(p.GetCov()[3]);
1145       hcovel[5]=(Double_t)(p.GetCov()[4]);
1146       hcovel[6]=(Double_t)(p.GetCov()[2]);
1147       hcovel[7]=(Double_t)(p.GetCov()[4]);
1148       hcovel[8]=(Double_t)(p.GetCov()[5]);
1149       hcov.SetRotation(hcovel);
1150       // now rotate in local system
1151       hcov.Multiply(svMatrix);
1152       hcov.MultiplyLeft(&svMatrix->Inverse());
1153       // now hcov is LOCAL COVARIANCE MATRIX
1154
1155       /// get original matrix of sens. vol.
1156       TGeoHMatrix *svOrigMatrix = m2[index]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1157       // modify global coordinates according with pre-aligment
1158       svOrigMatrix->LocalToMaster(pl,pg);
1159       // now rotate in local system
1160       hcov.Multiply(&svOrigMatrix->Inverse());
1161       hcov.MultiplyLeft(svOrigMatrix);
1162       // hcov is back in GLOBAL RF
1163       Float_t pcov[6];
1164       pcov[0]=hcov.GetRotationMatrix()[0];
1165       pcov[1]=hcov.GetRotationMatrix()[1];
1166       pcov[2]=hcov.GetRotationMatrix()[2];
1167       pcov[3]=hcov.GetRotationMatrix()[4];
1168       pcov[4]=hcov.GetRotationMatrix()[5];
1169       pcov[5]=hcov.GetRotationMatrix()[8];
1170
1171       p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1172       if (prn) printf("New global coordinates of measured point : X=%f  Y=%f  Z=%f \n",pg[0],pg[1],pg[2]);
1173       atps->AddPoint(npto,&p);
1174       if (prn) printf("Adding point[%d] = ( %f , %f , %f )     volid = %d\n",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] );
1175       if (prn) p.Print("");
1176
1177       npto++;
1178     }
1179
1180
1181     ////////////////////////////////////////////////////////////
1182     array = atps;
1183     treeout->Fill();
1184
1185     delete atps;
1186     atps=NULL;
1187
1188   } // end loop on tracks
1189
1190   pointsFile->Write();
1191   pointsFile->Close();
1192   tpainfile->Close();
1193
1194   return 0;
1195 }