TPC/AliBarrelRec_TPCparam.C

   1 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   2   //
   3   // Macro to create AliESDs.root using parametrized TPC tracking
   4   // and AliITStrackerSA (MI + SA)
   5   //
   6   // Input files:
   7   // - galice.root
   8   // - Kinematics.root
   9   // - TrackReferences.root
  10   // - ITS.RecPoints.root (use AliRecontruction class)
  11   // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
  12   //
  13   // A. Dainese - INFN Legnaro
  14   //
  15
  16
  17   Int_t  collcode = 1; // pp collisions
  18   Bool_t useMeanVtx = kFALSE;
  19
  20   AliGeomManager::LoadGeometry("geometry.root");
  21
  22   if (gAlice) {
  23     delete AliRunLoader::Instance();
  24     delete gAlice;
  25     gAlice=0;
  26   }
  27   AliRunLoader *rl = AliRunLoader::Open("galice.root");
  28   if (rl == 0x0) {
  29     cerr<<"Can not open session"<<endl;
  30     return;
  31   }
  32   Int_t retval = rl->LoadgAlice();
  33   if (retval) {
  34     cerr<<"LoadgAlice returned error"<<endl;
  35     delete rl;
  36     return;
  37   }
  38   retval = rl->LoadHeader();
  39   if (retval) {
  40     cerr<<"LoadHeader returned error"<<endl;
  41     delete rl;
  42     return;
  43   }
  44   gAlice=rl->GetAliRun();
  45
  46
  47
  48   // Get field from galice.root
  49   AliMagF *fiel = (AliMagF*)gAlice->Field();
  50   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
  51   // Set the conversion constant between curvature and Pt
  52   AliTracker::SetFieldMap(fiel,kTRUE);
  53
  54   /**** The TPC corner ********************/
  55
  56   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
  57   tpcTrackerPar.Init();
  58
  59   //**** Switch on the PID class (mandatory!)
  60   AliPID pid;
  61
  62   /**** The ITS corner ********************/
  63
  64   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  65   if (itsl == 0x0) {
  66     cerr<<"Cannot get the ITS loader"<<endl;
  67     return;
  68   }
  69   itsl->LoadRecPoints("read");
  70
  71   AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam();
  72   AliITSReconstructor::SetRecoParam(itsRecoParam);
  73
  74   // Instance of the ITS tracker
  75   AliITStrackerSA itsTracker(0);
  76   Int_t ITSclusters[6] = {1,1,1,1,1,1};
  77   itsTracker.SetLayersNotToSkip(ITSclusters);
  78
  79   // Primary vertex reconstruction in pp
  80   AliESDVertex *initVertex = 0;
  81   TFile *invtx = new TFile("AliESDVertexMean.root");
  82   if(collcode==1 && useMeanVtx) {
  83     // open the mean vertex
  84     initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  85   } else {
  86     Double_t pos[3]={0.5,0.5,0.};
  87     Double_t err[3]={3.,3.,5.};
  88     initVertex = new AliESDVertex(pos,err);
  89   }
  90   invtx->Close();
  91   delete invtx;
  92   AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
  93   vertexer->SetVtxStart(initVertex);
  94   vertexer->SetDebug(0);
  95   delete initVertex;
  96   initVertex=0;
  97
  98   /***** The TREE for ESD is born *****/
  99   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
 100   AliESDEvent *event=0; AliESDEvent *eventTPCin=0;
 101   event = new AliESDEvent();
 102   event->CreateStdContent();
 103   event->WriteToTree(esdTree);
 104
 105   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
 106   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
 107   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
 108
 109   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
 110   AliESDVertex *vertexSPD = new AliESDVertex();
 111   Char_t zver[100];
 112   Double_t vtx[3]={0,0,0};
 113   Double_t sigmavtx[3]={0.07,0.07,0.1};
 114
 115
 116   //<---------------------------------- The Loop over events begins
 117   TStopwatch timer;
 118   Int_t trc;
 119   for(Int_t i=firstEvent; i<=lastEvent; i++) {
 120
 121     cout<<" Processing event number : "<<i<<endl;
 122     //AliESDEvent *event = new AliESDEvent();
 123     event->SetRunNumber(gAlice->GetRunNumber());
 124     event->SetEventNumberInFile(i);
 125     event->SetMagneticField(gAlice->Field()->SolenoidField());
 126     rl->GetEvent(i);
 127
 128     //***** Primary vertex from SPD from file
 129     sprintf(zver,"Event%d/Vertex",i);
 130     vertexSPD = (AliESDVertex*)ppZ->Get(zver);
 131     if(!vertexSPD) {
 132       esdTree->Fill(); event->Reset();
 133       continue;
 134     }
 135     event->SetVertex(vertexSPD);
 136     vertexSPD->GetXYZ(vtx);
 137     vertexSPD->GetSigmaXYZ(sigmavtx);
 138
 139     //***** TPC tracking
 140     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
 141       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
 142       esdTree->Fill(); event->Reset();
 143       continue;
 144     }
 145
 146     // make a copy of the ESD at this stage
 147     eventTPCin = event;
 148
 149     //***** ITS tracking
 150     itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
 151     //    itsl->LoadRecPoints("read");
 152     TTree *itsTree=itsl->TreeR();
 153     if (!itsTree) {
 154       cerr<<"Can't get the ITS cluster tree !\n";
 155       esdTree->Fill(); event->Reset();
 156       return;
 157     }
 158     itsTracker.UnloadClusters();
 159     itsTracker.LoadClusters(itsTree);
 160     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
 161       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
 162       esdTree->Fill(); event->Reset();
 163       continue;
 164     }
 165
 166     // Bring kTPCin-tracks back to the TPC inner wall
 167     BackToTPCInnerWall(event,eventTPCin);
 168
 169     // refit inward in ITS:
 170     // - refit without vertex constraint
 171     // - propagate through beam pipe to local x = 0
 172     itsTracker.RefitInward(event);
 173
 174
 175     //***** Vertex from ESD tracks
 176     if(collcode==1) { // pp
 177       AliESDVertex *vertexTrks =
 178         (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 179       event->SetPrimaryVertex(vertexTrks);
 180     }
 181
 182     esdTree->Fill();
 183     event->Reset();
 184
 185   }//<-----------------------------------The Loop over events ends here
 186   timer.Stop(); timer.Print();
 187
 188   // The AliESDs.root is born
 189   TFile *ef = TFile::Open("AliESDs.root","RECREATE");
 190   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
 191
 192   //Write the tree and close everything
 193   esdTree->Write();
 194   delete esdTree;
 195   ef->Close();
 196
 197   delete vertexer;
 198   delete rl;
 199
 200   return;
 201 }
 202 //--------------------------------------------------------------------------
 203 void BackToTPCInnerWall(AliESDEvent *event,AliESDEvent *eventTPC) {
 204
 205   Int_t ntracks = eventTPC->GetNumberOfTracks();
 206   AliESDtrack *esdTrackTPC = 0;
 207
 208   // create relation between event and eventTPC
 209   Int_t labelsTPC[100000000];
 210   for(Int_t tr = 0; tr<ntracks; tr++) {
 211     esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
 212     labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
 213   }
 214
 215   ntracks = event->GetNumberOfTracks();
 216   AliESDtrack *esdTrack = 0;
 217   esdTrackTPC = 0;
 218   Int_t indexTPC;
 219
 220   // loop on tracks
 221   for(tr = 0; tr<ntracks; tr++) {
 222     esdTrack = (AliESDtrack*)event->GetTrack(tr);
 223     // set to kITSout the tracks that don't have kTPCin
 224     // (they've been found by AliITStrackerSA)
 225     if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
 226       esdTrack->SetStatus(AliESDtrack::kITSout);
 227       continue;
 228     }
 229
 230     // skip tracks that don't have kITSin
 231     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
 232
 233     indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
 234     esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
 235
 236     AliITStrackMI *itsTrack = 0;
 237     try {
 238       itsTrack = new AliITStrackMI(*esdTrackTPC);
 239       esdTrack = 0;
 240     }
 241     catch (const Char_t *msg) {
 242         Warning("ToTPCInnerWall",msg);
 243         continue;
 244     }
 245     itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
 246     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
 247
 248     delete itsTrack;
 249   } // end loop on tracks
 250
 251   return;
 252 }
 253 //--------------------------------------------------------------------------