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 AliITStrackerV2
   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 - LNL
  14   //
  15
  16    Int_t  collcode = 1; // pp collisions
  17    Bool_t useMeanVtx = kTRUE;
  18
  19    if (gAlice) {
  20       delete gAlice->GetRunLoader();
  21       delete gAlice;
  22       gAlice=0;
  23    }
  24
  25    AliRunLoader *rl = AliRunLoader::Open("galice.root");
  26    if (rl == 0x0) {
  27       cerr<<"Can not open session"<<endl;
  28       return;
  29    }
  30    Int_t retval = rl->LoadgAlice();
  31    if (retval) {
  32       cerr<<"LoadgAlice returned error"<<endl;
  33       delete rl;
  34       return;
  35    }
  36    retval = rl->LoadHeader();
  37    if (retval) {
  38       cerr<<"LoadHeader returned error"<<endl;
  39       delete rl;
  40       return;
  41    }
  42    gAlice=rl->GetAliRun();
  43
  44
  45    TDatabasePDG *DataBase = TDatabasePDG::Instance();
  46
  47    // Get field from galice.root
  48    AliMagF *fiel = (AliMagF*)gAlice->Field();
  49    Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
  50    // Set the conversion constant between curvature and Pt
  51    AliTracker::SetFieldMap(fiel,kTRUE);
  52
  53    /**** The TPC corner ********************/
  54
  55    AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
  56    tpcTrackerPar.Init();
  57
  58    //**** Switch on the PID class (mandatory!)
  59    AliPID pid;
  60
  61    /**** The ITS corner ********************/
  62
  63    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
  64    if (itsl == 0x0) {
  65       cerr<<"Cannot get the ITS loader"<<endl;
  66       return;
  67    }
  68    itsl->LoadRecPoints("read");
  69
  70    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  71    if (!dITS) {
  72       cerr<<"Cannot find the ITS detector !"<<endl;
  73       return;
  74    }
  75    AliITSgeom *geom = dITS->GetITSgeom();
  76
  77    //An instance of the ITS tracker
  78    AliITStrackerV2 itsTracker(geom);
  79    Int_t ITSclusters[6] = {1,1,1,1,1,1};
  80    itsTracker.SetLayersNotToSkip(ITSclusters);
  81
  82    // Primary vertex reconstruction in pp
  83    AliESDVertex *initVertex = 0;
  84    if(collcode==1 && useMeanVtx) {
  85      // open the mean vertex
  86      TFile *invtx = new TFile("AliESDVertexMean.root");
  87      initVertex = (AliESDVertex*)invtx->Get("vtxmean");
  88      invtx->Close();
  89      delete invtx;
  90    } else {
  91      Double_t pos[3]={0.5,0.5,0.};
  92      Double_t err[3]={3.,3.,5.};
  93      initVertex = new AliESDVertex(pos,err);
  94    }
  95    AliVertexerTracks *vertexer = new AliVertexerTracks;
  96    vertexer->SetVtxStart(initVertex);
  97
  98
  99    /***** The TREE for ESD is born *****/
 100    TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
 101    AliESD *event=0;
 102    esdTree->Branch("ESD","AliESD",&event);
 103
 104    if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
 105    if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
 106    cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
 107
 108    TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
 109    AliESDVertex *vertexSPD = new AliESDVertex();
 110    Char_t zver[100];
 111    Double_t vtx[3];
 112    Double_t cvtx[6];
 113
 114
 115    //<---------------------------------- The Loop over events begins
 116    TStopwatch timer;
 117    Int_t trc;
 118    for(Int_t i=firstEvent; i<=lastEvent; i++) {
 119
 120      cerr<<" Processing event number : "<<i<<endl;
 121      AliESD *event = new AliESD();
 122      event->SetRunNumber(gAlice->GetRunNumber());
 123      event->SetEventNumber(i);
 124      event->SetMagneticField(gAlice->Field()->SolenoidField());
 125      rl->GetEvent(i);
 126
 127      //***** Primary vertex from SPD from file
 128      sprintf(zver,"Event%d/Vertex",i);
 129      vertexSPD = (AliESDVertex*)ppZ->Get(zver);
 130      if(!vertexSPD) {
 131        esdTree->Fill(); delete event;
 132        continue;
 133      }
 134      event->SetVertex(vertexSPD);
 135      vertexSPD->GetXYZ(vtx);
 136      vertexSPD->GetCovMatrix(cvtx);
 137
 138      //***** TPC tracking
 139      if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
 140        printf("exiting TPC tracker with code %d in event %d\n",trc,i);
 141        esdTree->Fill(); delete event;
 142        continue;
 143      }
 144
 145      //***** ITS tracking
 146      itsTracker.SetVertex(vtx,cvtx);
 147      TTree *itsTree=itsl->TreeR();
 148      if (!itsTree) {
 149         cerr<<"Can't get the ITS cluster tree !\n";
 150         esdTree->Fill(); delete event;
 151         return;
 152      }
 153      itsTracker.UnloadClusters();
 154      itsTracker.LoadClusters(itsTree);
 155      if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
 156        printf("exiting ITS tracker with code %d in event %d\n",trc,i);
 157        esdTree->Fill(); delete event;
 158        continue;
 159      }
 160
 161      //***** Propagate kITSin tracks to local x = 0 (through beam pipe)
 162      ToLocalX0(event);
 163
 164
 165      //***** Vertex from ESD tracks
 166      if(collcode==1) { // pp
 167        AliESDVertex *vertexTrks =
 168          (AliESDVertex*)vertexer->FindPrimaryVertex(event);
 169        event->SetPrimaryVertex(vertexTrks);
 170      }
 171
 172      esdTree->Fill();
 173      delete event;
 174
 175    }//<-----------------------------------The Loop over events ends here
 176    timer.Stop(); timer.Print();
 177
 178    // The AliESDs.root is born
 179    TFile *ef = TFile::Open("AliESDs.root","RECREATE");
 180    if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
 181
 182    //Write the TREE and close everything
 183    esdTree->Write();
 184    delete esdTree;
 185    ef->Close();
 186
 187    delete vertexer;
 188    delete initVertex;
 189    delete rl;
 190
 191    return;
 192 }
 193 //--------------------------------------------------------------------------
 194 Int_t ToLocalX0(AliESD *esd) {
 195
 196   Int_t ntracks = esd->GetNumberOfTracks();
 197   AliESDtrack *esdTrack = 0;
 198
 199   // loop on tracks
 200   for(Int_t tr = 0; tr<ntracks; tr++) {
 201     esdTrack = (AliESDtrack*)esd->GetTrack(tr);
 202     // ask for kITSin
 203     if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
 204     AliITStrackV2 *itsTrack = 0;
 205     try {
 206       itsTrack = new AliITStrackV2(*esdTrack);
 207       esdTrack = 0;
 208     }
 209     catch (const Char_t *msg) {
 210         Warning("FindPrimaryVertexForCurrentEvent",msg);
 211         continue;
 212     }
 213     // propagate track to beam pipe (0.8 mm of Be)
 214     if(itsTrack->GetX()>3.) itsTrack->PropagateTo(3.,0.0023,65.19);
 215     // propagate track to (0,0)
 216     itsTrack->PropagateTo(0.,0.,0.);
 217     itsTrack->UpdateESDtrack(AliESDtrack::kITSin);
 218     esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
 219     delete itsTrack;
 220   } // end loop on tracks
 221
 222   return 0;
 223 }